Question A

M = 1000
Nlist = c(1,10,100)

theDensities = lapply(1:length(Nlist), function(i) {
  n = Nlist[i]
simdatasets = matrix(
  rbinom(n*M, size=1, prob=0.2),
    #rpois(n*M, 1),
    #rexp(n*M),
    #rnorm(n*M),
    #rcauchy(n*M),
    nrow=n)
  str(simdatasets)
  
  simMeans = apply(simdatasets, 2, mean)
  str(simMeans)
  # summary(simMeans)
  # summary(colMeans(simdatasets))  ### same thing.
  density(simMeans, kernel="triangular", bw=0.02) 
  ### changing the bandwidth to make it more narrow.
}
)
##  int [1, 1:1000] 0 0 0 0 0 0 0 1 0 0 ...
##  num [1:1000] 0 0 0 0 0 0 0 1 0 0 ...
##  int [1:10, 1:1000] 1 0 0 0 0 0 1 0 1 1 ...
##  num [1:1000] 0.4 0.2 0.3 0.1 0.3 0.3 0.1 0.3 0.1 0.2 ...
##  int [1:100, 1:1000] 1 0 1 1 0 0 0 0 0 1 ...
##  num [1:1000] 0.26 0.18 0.15 0.22 0.19 0.15 0.24 0.19 0.2 0.21 ...
theXranges = sapply(theDensities, function(theDensity) range(theDensity$x))
theYranges = sapply(theDensities, function(theDensity) range(theDensity$y))

for(i in 1:length(Nlist)) {
  if(i==1)
    plot(theDensities[[i]],  type="b", col=2, pch=as.character(i)
          , xlim=c(min(theXranges),max(theXranges))
          , ylim=c(min(theYranges),max(theYranges))
      )
  else points(theDensities[[i]], type="b", col=i+1, pch=as.character(i))  
}

#legend()
legend(x="topright",
       legend=paste("Nlist =", Nlist),
       pch=as.character(1:length(Nlist)),
       lty=rep(1,length(Nlist)),
       col=1+1:length(Nlist)
)

M = 1000
Nlist = c(1,10,100)

theDensities = lapply(1:length(Nlist), function(i) {
  n = Nlist[i]
simdatasets = matrix(
  #rbinom(n*M, size=1, prob=0.2),
    rpois(n*M, 1),
    #rexp(n*M),
    #rnorm(n*M),
    #rcauchy(n*M),
    nrow=n)
  str(simdatasets)
  
  simMeans = apply(simdatasets, 2, mean)
  str(simMeans)
  # summary(simMeans)
  # summary(colMeans(simdatasets))  ### same thing.
  density(simMeans, kernel="triangular", bw=0.02) 
  ### changing the bandwidth to make it more narrow.
}
)
##  int [1, 1:1000] 1 2 2 0 1 3 2 2 1 0 ...
##  num [1:1000] 1 2 2 0 1 3 2 2 1 0 ...
##  int [1:10, 1:1000] 0 0 0 3 3 1 0 0 2 1 ...
##  num [1:1000] 1 0.7 0.5 0.6 0.9 1.1 0.9 0.8 0.6 0.8 ...
##  int [1:100, 1:1000] 2 2 0 0 3 1 0 2 1 0 ...
##  num [1:1000] 0.92 1.08 0.97 0.98 1.13 0.96 1.02 0.99 0.98 1 ...
theXranges = sapply(theDensities, function(theDensity) range(theDensity$x))
theYranges = sapply(theDensities, function(theDensity) range(theDensity$y))

for(i in 1:length(Nlist)) {
  if(i==1)
    plot(theDensities[[i]],  type="b", col=2, pch=as.character(i)
          , xlim=c(min(theXranges),max(theXranges))
          , ylim=c(min(theYranges),max(theYranges))
      )
  else points(theDensities[[i]], type="b", col=i+1, pch=as.character(i))  
}

#legend()
legend(x="topright",
       legend=paste("Nlist =", Nlist),
       pch=as.character(1:length(Nlist)),
       lty=rep(1,length(Nlist)),
       col=1+1:length(Nlist)
)

M = 1000
Nlist = c(1,10,100)

theDensities = lapply(1:length(Nlist), function(i) {
  n = Nlist[i]
simdatasets = matrix(
  #rbinom(n*M, size=1, prob=0.2),
    #rpois(n*M, 1),
   rexp(n*M),
    #rnorm(n*M),
    #rcauchy(n*M),
    nrow=n)
  str(simdatasets)
  
  simMeans = apply(simdatasets, 2, mean)
  str(simMeans)
  # summary(simMeans)
  # summary(colMeans(simdatasets))  ### same thing.
  density(simMeans, kernel="triangular", bw=0.02) 
  ### changing the bandwidth to make it more narrow.
}
)
##  num [1, 1:1000] 0.0599 0.2801 1.4005 1.1369 0.2974 ...
##  num [1:1000] 0.0599 0.2801 1.4005 1.1369 0.2974 ...
##  num [1:10, 1:1000] 1.38 0.258 0.668 0.433 0.392 ...
##  num [1:1000] 0.627 1.308 0.866 1.481 0.822 ...
##  num [1:100, 1:1000] 0.354 0.372 0.812 0.247 0.983 ...
##  num [1:1000] 0.97 1.013 0.857 0.863 1.16 ...
theXranges = sapply(theDensities, function(theDensity) range(theDensity$x))
theYranges = sapply(theDensities, function(theDensity) range(theDensity$y))

for(i in 1:length(Nlist)) {
  if(i==1)
    plot(theDensities[[i]],  type="b", col=2, pch=as.character(i)
          , xlim=c(min(theXranges),max(theXranges))
          , ylim=c(min(theYranges),max(theYranges))
      )
  else points(theDensities[[i]], type="b", col=i+1, pch=as.character(i))  
}

#legend()
legend(x="topright",
       legend=paste("Nlist =", Nlist),
       pch=as.character(1:length(Nlist)),
       lty=rep(1,length(Nlist)),
       col=1+1:length(Nlist)
)

M = 1000
Nlist = c(1,10,100)

theDensities = lapply(1:length(Nlist), function(i) {
  n = Nlist[i]
simdatasets = matrix(
  #rbinom(n*M, size=1, prob=0.2),
    #rpois(n*M, 1),
    #rexp(n*M),
   rnorm(n*M),
    #rcauchy(n*M),
    nrow=n)
  str(simdatasets)
  
  simMeans = apply(simdatasets, 2, mean)
  str(simMeans)
  # summary(simMeans)
  # summary(colMeans(simdatasets))  ### same thing.
  density(simMeans, kernel="triangular", bw=0.02) 
  ### changing the bandwidth to make it more narrow.
}
)
##  num [1, 1:1000] -0.8555 -0.7085 -1.0275 0.0759 0.2853 ...
##  num [1:1000] -0.8555 -0.7085 -1.0275 0.0759 0.2853 ...
##  num [1:10, 1:1000] 0.6025 0.0763 -0.0311 1.7558 -0.0246 ...
##  num [1:1000] 0.45929 0.00269 -0.37664 -0.27088 0.42689 ...
##  num [1:100, 1:1000] -2.02 -0.185 -0.389 0.764 -0.277 ...
##  num [1:1000] -0.0167 -0.0322 -0.0789 -0.0273 0.0526 ...
theXranges = sapply(theDensities, function(theDensity) range(theDensity$x))
theYranges = sapply(theDensities, function(theDensity) range(theDensity$y))

for(i in 1:length(Nlist)) {
  if(i==1)
    plot(theDensities[[i]],  type="b", col=2, pch=as.character(i)
          , xlim=c(min(theXranges),max(theXranges))
          , ylim=c(min(theYranges),max(theYranges))
      )
  else points(theDensities[[i]], type="b", col=i+1, pch=as.character(i))  
}

#legend()
legend(x="topright",
       legend=paste("Nlist =", Nlist),
       pch=as.character(1:length(Nlist)),
       lty=rep(1,length(Nlist)),
       col=1+1:length(Nlist)
)

M = 1000
Nlist = c(1,10,100)

theDensities = lapply(1:length(Nlist), function(i) {
  n = Nlist[i]
simdatasets = matrix(
  #rbinom(n*M, size=1, prob=0.2),
    #rpois(n*M, 1),
    #rexp(n*M),
    #rnorm(n*M),
    rcauchy(n*M),
    nrow=n)
  str(simdatasets)
  
  simMeans = apply(simdatasets, 2, mean)
  str(simMeans)
  # summary(simMeans)
  # summary(colMeans(simdatasets))  ### same thing.
  density(simMeans, kernel="triangular", bw=0.02) 
  ### changing the bandwidth to make it more narrow.
}
)
##  num [1, 1:1000] -2.435 -0.287 -0.162 -1.328 0.402 ...
##  num [1:1000] -2.435 -0.287 -0.162 -1.328 0.402 ...
##  num [1:10, 1:1000] -7.799 9.902 3.177 -2.214 -0.365 ...
##  num [1:1000] -0.0152 -2.6944 0.4839 -19.754 -0.328 ...
##  num [1:100, 1:1000] -2.152 -0.112 0.595 16.616 28.079 ...
##  num [1:1000] -2.353 1.88 -0.591 0.422 1.251 ...
theXranges = sapply(theDensities, function(theDensity) range(theDensity$x))
theYranges = sapply(theDensities, function(theDensity) range(theDensity$y))

for(i in 1:length(Nlist)) {
  if(i==1)
    plot(theDensities[[i]],  type="b", col=2, pch=as.character(i)
          , xlim=c(min(theXranges),max(theXranges))
          , ylim=c(min(theYranges),max(theYranges))
      )
  else points(theDensities[[i]], type="b", col=i+1, pch=as.character(i))  
}

#legend()
legend(x="topright",
       legend=paste("Nlist =", Nlist),
       pch=as.character(1:length(Nlist)),
       lty=rep(1,length(Nlist)),
       col=1+1:length(Nlist)
)

##  Discussion  ##

Discussion

This experiment is to observe the effect of sample size on each distribution of the sample mean. The distribution of the sample mean is a probability distribution for all possible values of a sample mean, computed from a sample of size n. In this question, I generate M random samples with various sample size, and plot the density of the mean of each random sample for a given distribution. Interestingly, in all of the random simulation data except Cauchy distribution, the distribution of the smaple mean shows gaussian curve as the sample size is large. This is the application of Central Limit Theorum; when sample size grows, the distribution of sample mean approaches normal distribution (gausian curve).

In the first random binomial simulation data, the 1000 mean was estimated close to 0.2 as the sample size (Nlist) is larger. The distribution of the sample mean was dispersed as the sample size is small.

The second random Poisson simulation data show discreted distribution of the sample mean when the sample size is small. It is because Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of space.

The third random exponential simulation data show the sample mean converges to 1 as the sample size is large. When, sample size is small (Nlist=10), I observed the sample mean is in the range of 0~3 approximately.

The forth random normal simulation data should the sample mean that converges to 0. The three plots of the sample mean disbribution show gaussian curve.

The fifth random cauchy simulation data show very sharp convergence of sample mean at 0 when the sample size is 100. The mean and variance of Cauchy distribution are undefined as Cauchy distribution has infinite variance. Therefore, the distribution of the sample mean even when the sample size is large (Nlist=100) is very dispersed and don’t follow normal distribution. But, interestingly, the sample mean density for mean=0 is very high regardless of the sample size.

##  Question B  ##

Question B

The likelihood function of θ as a function of the data with random variables x1,…,xn is defined as lik(θ) = f(x1,…,xn | θ). If the Xi are i.i.d., the joint density is the product of the marginal densities: lik(θ) = ∏f(xi|θ). Typically, it is more convenient to work with the (natural) logarithm of the likelihood, so-called log likelihood. In the i.i.d. case, this is given by l(θ) = ∑log(f(xi|θ). Then, this function becomes a sum of a series of log probabilities. Since the data are random, log likelihood, l(θ), is random varialbles too; it’s a random variable made up of a sum of other random variables.

The Central Limit Theorum says that the sum of random variables from a large sample size will have approximately normal distribution. Therefore, the log-likelihood function, l(θ), will have approximately normal distribution under the condition that sample size is large enough, because log-likelihood function is a sum of random variables. The Central Limit Theorem gives a first order approximation to the distribution of the log-likelihood function under the regularity conditions. According to Central Limit Theorem, we can establish the asymptotic normality of the gradient of the log- likelihood, so as to be able to characterize the likelihood itself.

Question C

5% of the normal distribution lies outsides an interval approximately +2s -2s about the mean; to what extend is this true? Where are the limits corresponding to 1%, 0.5%, and 0.1%?

I tried to find answers, but I couldn’t find clearn answer yet.