# Example of LLN (Law of Large Number)
# t: sample size
# n: repetition time

# clean up workspace
rm(list=ls(all=TRUE))
# chiq-square(df=1), mu=1
fun_LLN <-function(t,n){
  x <- numeric(n)
  for (i in 1:n){
    x[i] <-mean(rchisq(t,df=1))
  }
  x
}

par(mfrow=c(2,2))
hist(fun_LLN(50,1000),xlim=range (0,2),freq=FALSE,main='T=50',xlab='Sample Mean')
hist(fun_LLN(100,1000),xlim=range (0,2),freq=FALSE,main='T=100',xlab='Sample Mean')
hist(fun_LLN(300,1000),xlim=range (0,2),freq=FALSE,main='T=300',xlab='Sample Mean')
hist(fun_LLN(1000,1000),xlim=range (0,2),freq=FALSE,main='T=1000',xlab='Sample Mean')


# clean up workspace
rm(list=ls(all=TRUE))
# student t-distribution(df=5), mu=0
fun_LLN <-function(t,n){
  x <- numeric(n)
  for (i in 1:n){
    x[i] <-mean(rt(t,df=5))
  }
  x
}

par(mfrow=c(2,2))
hist(fun_LLN(50,1000),xlim=range (-1,1),freq=FALSE,main='T=50',xlab='Sample Mean')
hist(fun_LLN(100,1000),xlim=range (-1,1),freq=FALSE,main='T=100',xlab='Sample Mean')
hist(fun_LLN(300,1000),xlim=range (-1,1),freq=FALSE,main='T=300',xlab='Sample Mean')
hist(fun_LLN(1000,1000),xlim=range (-1,1),freq=FALSE,main='T=1000',xlab='Sample Mean')


# clean up workspace
rm(list=ls(all=TRUE))
# student t-distribution(df=1), mu is not defined
fun_LLN <-function(t,n){
  x <- numeric(n)
  for (i in 1:n){
    x[i] <-mean(rt(t,df=1))
  }
  x
}

par(mfrow=c(2,2))
hist(fun_LLN(50,1000),xlim=range (-5,5),freq=FALSE,main='T=50',xlab='Sample Mean')
hist(fun_LLN(100,1000),xlim=range (-5,5),freq=FALSE,main='T=100',xlab='Sample Mean')
hist(fun_LLN(300,1000),xlim=range (-5,5),freq=FALSE,main='T=300',xlab='Sample Mean')
hist(fun_LLN(1000,1000),xlim=range (-5,5),freq=FALSE,main='T=1000',xlab='Sample Mean')

#######################################################
#######################################################
# Example of LLN (Central Limit Theorem)
# t: sample size
# n: repetition time

# chiq-square(df=1), mu=1, var=2
# t: sample size
# n: repetition time
# clean up workspace
rm(list=ls(all=TRUE))
fun_CLT <-function(t,n){
  x <- numeric(n)
  s.x <-numeric(n)
  k <- 1
  for (i in 1:n){
    x[i] <-mean(rchisq(t,df=1)) 
    s.x[i] <-sqrt(t)*(x[i]-k)/(sqrt(2*k))
  }
  s.x
}

par(mfrow=c(2,2))
xtemp <-fun_CLT(50,1000)
hist(xtemp,xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=50',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(100,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=100',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(300,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=300',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(1000,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=1000',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")



# student t-distribution(df=5), mu=0, var=5/3
# t: sample size
# n: repetition time
# clean up workspace
rm(list=ls(all=TRUE))
fun_CLT <-function(t,n){
  x <- numeric(n)
  s.x <-numeric(n)
  k <- 0
  var <-5/(5-2)
  for (i in 1:n){
    x[i] <-mean(rt(t,df=5)) 
    s.x[i] <-sqrt(t)*(x[i]-k)/(sqrt(var))
  }
  s.x
}

par(mfrow=c(2,2))
xtemp <-fun_CLT(50,1000)
hist(xtemp,xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=50',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(100,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=100',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(300,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=300',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(1000,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=1000',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")


# student t-distribution(df=2), mu=0, var is not defined
# t: sample size
# n: repetition time
# clean up workspace
rm(list=ls(all=TRUE))
fun_CLT <-function(t,n){
  x <- numeric(n)
  s.x <-numeric(n)
  sd_x <-numeric(n)
  mu <- 0
  for (i in 1:n){
    sample <-rt(t,df=2)
    x[i] <-mean(sample) 
    sd_x[i] <-sd(sample)
    s.x[i] <-sqrt(t)*(x[i]-mu)/(sqrt(sd_x[i]))
  }
  s.x
}

par(mfrow=c(2,2))
xtemp <-fun_CLT(50,1000)
hist(xtemp,xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=50',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(100,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=100',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(300,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=300',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")

hist(fun_CLT(1000,1000),xlim=range (-4,4),ylim=range(0,0.45),freq=FALSE,main='T=1000',xlab='Sample Mean')
lines(density(xtemp))
curve(dnorm(x, mean=0, sd=1), add=TRUE,col="red")