### This is the script file for simulating the LLN, CLT, and
### independence of the sample mean and sample variance.

#################################################################################
# cummean function: calculates the cummulative mean of a vactor.

cummean <- function(x){
	n <- length(x)
	y <- numeric(n)
	z <- c(1:n)
	y <- cumsum(x)
	y <- y/z
	return(y)
}

# LLN

n <- 10000

z <- rnorm(n)

hist(z, main= 'Standard Normal Random Deviates')

x <- seq(1,n,1)

y <- cummean(z)

plot(x,y,type= 'l',main= 'Convergence Plot')

epsilon <- 0.1

Count <- 0

for(i in 1:n){
	if(abs(y[i]) > epsilon) Count <- Count + 1
}

Count/n 	# approximate probability that abs(y[i]-mu) > epsilon

#################################################################################
# CLT

# Normal

n <- 30		# sample size
k <- 1000		# number of samples

mu <- 5
sigma <- 2
SEM <- sigma/sqrt(n)

x <- matrix(rnorm(n*k,mu,sigma),n,k)		# This gives a matrix with the samples 
							# down the columns.

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

x.down <- mu - 4*SEM
x.up <- mu + 4*SEM

y.up <- 1.5

hist(x.mean, prob = T, xlim = c(x.down,x.up), ylim = c(0,y.up), main= 'Sampling 
	distribution of the sample mean, Normal case')

par(new = T)

x <- seq(x.down,x.up,0.01)

y <- dnorm(x,mu,SEM)

plot(x, y, type= 'l', xlim = c(x.down,x.up), ylim = c(0,y.up))

# Exponential

n <- 30		# sample size
k <- 1000		# number of samples

theta <- 3

lambda <- 1/theta

mu <- theta
sigma <- theta
SEM <- sigma/sqrt(n)


x <- matrix(rexp(n*k,lambda),n,k)		# This gives a matrix with the 
							# samples down the columns.

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

x.down <- 0
x.up <- 6

y.up <- 1

hist(x.mean, prob = T, xlim= c(x.down,x.up), ylim = c(0,y.up), main = 'Sampling 
	distribution of the sample mean, Normal case')

par(new = T)

x <- seq(x.down,x.up,0.01)

y <- dnorm(x,mu,SEM)

plot(x, y, type = 'l', xlim = c(x.down,x.up), ylim = c(0,y.up))

#################################################################################
# Independence of x.bar and s.

# Normal

n <- 30		# sample size
k <- 1000		# number of samples

x <- matrix(rnorm(n*k,mu,sigma),n,k)		# This gives a matrix with the samples 
							# down the columns.

x.mean <- apply(x,2,mean)
x.sd <- sqrt(apply(x,2,var))

plot(x.mean,x.sd)

cor(x.mean,x.sd)

# Exponential

n <- 30		# sample size
k <- 1000		# number of samples

theta <- 3

lambda <- 1/theta

mu <- theta
sigma <- theta
SEM <- sigma/sqrt(n)

x <- matrix(rexp(n*k,lambda),n,k)		# This gives a matrix with the samples 
							# down the columns.

x.mean <- apply(x,2,mean)
x.sd <- sqrt(apply(x,2,var))

plot(x.mean,x.sd)

cor(x.mean,x.sd)