### What is the difference between 

### s^2 = sum(x_i - x.bar)^2 / (n - 1)	(unbiased)
###
### t^2 = sum(x_i - x.bar)^2 / n			(biased)

### when estimating the variance of the exponetial?

# population parameters

mu <- 10
sigma <- 5
sigma2 <- sigma^2

# 100 samples of size n = 30

x <- matrix(rnorm(3000,mu,sigma),nrow=30)

# compute the s^2 for each sample

x.s2 <- apply(x,2,var)

# plot the sampling distribution

par(mfrow=c(2,2))

hist(x.s2,probability = T,main="sampling distribution of s2")

plot(density(x.s2),main="sampling distribution of x.s2",type="l")

# compute the mean, bias, se, rmse

mean(x.s2)

bias.s2 <- mean(x.s2) - sigma2
bias.s2

se.s2 <- stdev(x.s2)
se.s2

RMSE.s2 <- sqrt(bias.s2^2 + se.s2^2)
RMSE.s2

# function to compute t2

t2 <- function(x){
	t2 <- (n-1)/n*var(x)
}

# compute the t^2 for each sample

x.t2 <- apply(x,2,t2)

# plot the sampling distribution

hist(x.t2,probability = T,main="sampling distribution of t2")

plot(density(x.t2),main="sampling distribution of x.t2",type="l")

# compute the mean, bias, se, rmse

mean(x.t2)

bias.t2 <- mean(x.t2) - sigma2
bias.t2

se.t2 <- stdev(x.t2)
se.t2

sqrt(bias.t2^2 + se.t2^2)