#################################################################
### Simulation study: (Empirical) distribution of sample mean ###
#################################################################

# 'True' mean and standard deviation
mu <- 178
sigma <- 12

# Sample size
n <- 10

# Simulate normal distributed X_i for n = 10
x <- rnorm(n = n, mean = mu, sd = sigma)
round(x)

# Empirical density
hist(x, prob = TRUE, col = 'blue')

# Compute sample mean
mean(x)

# Repeat the simulated sampling many times (1000 samples)
mat <- replicate(1000, rnorm(n = n, mean = mu, sd = sigma))

# Compute the sample mean for each sample
xbar <- apply(mat, 2, mean)
xbar

# See the distribution of the sample means
hist(xbar, prob = TRUE, col = 'blue', nclass = 10, xlab = "sample mean",
     main = "Empirical distribution of the sample mean")

# Empirical mean of sample means
mean(xbar)


#####################################################################
### Density function for t-distribution with 9 degrees of freedom ###
### plotted together with that of a standard normal distribution  ###
#####################################################################

x <- seq(-4, 4, by = 0.01)
plot(x, dt(x, df = 9), type = "l", col = "red", ylab = "Density(x)")
lines(x, dnorm(x), type = "l")
text(2.5, 0.3,"Black: N(0,1)")
text(3, 0.1,"Red: t(9)", col = "red")

#############################################################
### Example: Heights                                      ###
### Computing a 99% CI for mu using the 't.test' function ###
#############################################################

# Data
x <- c(168, 161, 167, 179, 184, 166, 198, 187, 191, 179)

# 99% CI for mu
t.test(x, conf.level = 0.99)

################################################################################
### CLT Example: Mean of independent and U(0,1)-distributed random variables ###
################################################################################

n <- 1 # Sample size
k <- 1000 # No. of samples (i.e. no. of means to be computed)

# Simulations from U(0,1)-distribution (k = 1000 samples, each of size n = 1)
u <- matrix(runif(k*n), ncol = n)

# Empirical density of means
hist(apply(u, 1, mean), col = "blue", main = "n = 1", xlab = "Means", prob = TRUE)

n <- 2 # Sample size
k <- 1000 # No. of samples (i.e. no. of means to be computed)

# Simulations from U(0,1)-distribution (k = 1000 samples, each of size n = 2)
u <- matrix(runif(k*n), ncol = n)

# Empirical density of means
hist(apply(u, 1, mean), col = "blue", main = "n = 2", xlab = "Means", xlim = c(0,1), prob = TRUE)

n <- 6 # Sample size
k <- 1000 # No. of samples (i.e. no. of means to be computed)

# Simulations from U(0,1)-distribution (k = 1000 samples, each of size n = 6)
u <- matrix(runif(k*n), ncol = n)

# Empirical density of means
hist(apply(u, 1, mean), col = "blue", main = "n = 6", xlab = "Means", xlim = c(0,1), prob = TRUE)

n <- 30 # Sample size
k <- 1000 # No. of samples (i.e. no. of means to be computed)

# Simulations from U(0,1)-distribution (k = 1000 samples, each of size n = 30)
u <- matrix(runif(k*n), ncol = n)

# Empirical density of means
hist(apply(u, 1, mean), col = "blue", main = "n = 30", xlab = "Means", xlim = c(0,1), prob = TRUE)


#####################################################################
### Simulation study: (Empirical) distribution of sample variance ###
#####################################################################

# 'True' mean and standard deviation
mu <- 178
sigma <- 12

# Sample size
n <- 10

# Simulate normal distributed X_i for n = 10
x <- rnorm(n = n, mean = mu, sd = sigma)
round(x)

# Empirical density
hist(x, prob = TRUE, col = 'blue')

# Compute sample mean
mean(x)

# Repeat the simulated sampling many times (1000 samples)
mat <- replicate(1000, rnorm(n = n, mean = mu, sd = sigma))

# Compute the sample variance for each sample
xbar <- apply(mat, 2, var)
xbar

# See the empirical distribution of the sample variances
hist((n-1)*xbar/sigma^2, prob = TRUE, col = 'blue', nclass = 30,
     main = "Empirical pdf of sample variance", xlab = "sample variance")

# Compare with ChiSquare pdf
xp <- seq(0,40,0.1)
lines(xp, dchisq(xp,df=n-1), col = "red", lw = 3)
text(25, 0.06,"ChiSquare(df=9)", col = "red")