Monte Carlo Methods

MC simulation and the boostrap

1 Some data (50 iid standard normal draws)

x = c(-0.027,-0.530,-1.582,-0.431,-0.248,0.478,0.697,1.484,-0.514,1.312,-0.584,0.854,0.473,
      -0.691,-0.651,-0.596,0.557,-0.950,-0.006,1.162,0.444,-1.156,-0.880,-1.073,-0.332,
      0.824,0.512,0.311,-1.027,-1.069,2.307,-0.126,1.903,-0.197,-0.581,0.170,0.236,0.481,
      1.246,0.341,0.664,-0.769,1.017,-0.353,-0.074,-0.515,-0.997,1.073,-0.110,2.526)

1.1 Estimation of a few population quantities

Assuming that the above data is a random sample of size \(n=50\) from a population with mean \(\mu\) and variance \(\sigma^2\). Assuming the data is Gaussian (normally distributed), we can show that the maximum likelihood estimators (MLEs) of \(\mu\) and \(\sigma^2\) are, respectively \[ {\widehat \mu}_{MLE} = \bar{x}_n = \frac{1}{n} \sum_{i=1}^n x_i \ \ \ \mbox{and} \ \ \ {\widehat \sigma}_{MLE}^2 = \frac{1}{n} \sum_{i=1}^n (x_i-\bar{x}_n)^2, \] such that \({\widehat \sigma}_{MLE}=\sqrt{{\widehat \sigma}_{MLE}^2}\). Finally , the mean absolute deviation is computed as \[ MAD = \frac{1}{n} \sum_{i=1}^n |x_i-\bar{x}_n|. \]

n = length(x)
mu.mle = mean(x)
sigma2.mle = mean((x-mu.mle)^2)
sigma.mle = sqrt(sigma2.mle)
mad = mean(abs(x-mu.mle))

 MLE for the population mean: 0.10006 
 MLE for the population variance: 0.8334424 
 MLE for the population standard deviation: 0.9129307 
 Mean absolute deviation (MAD): 0.7508248 

2 Sampling distributions

A few results are well known and easy to show:

• \({\widehat \mu}_{MLE}|\mu,\sigma^2 \sim N(\mu,\sigma^2/n)\), so the 95% confidence interval for \(\mu\) is given by \({\widehat \mu}_{MLE} \pm 1.95996 \sigma/\sqrt{n}\).

• \({\widehat \mu}_{MLE}|\mu,{\widehat \sigma}^2 \sim t_{n-1}(\mu,{\widehat \sigma}^2/n)\), so the 95% confidence interval for \(\mu\) is given by \({\widehat \mu}_{MLE} \pm 2.00958{\widehat \sigma}/\sqrt{n}\).

• \({\widehat \sigma}_{MLE}^2|\sigma^2 \sim \frac{\sigma^2}{n} \chi_{n-1}^2\). Similarly, a 95% confidence interval for \(\sigma^2\) is given by \[ [n{\widehat \sigma}_{MLE}^2/\chi^2_{0.975,n-1}; n{\widehat \sigma}_{MLE}^2/\chi^2_{0.025,n-1}] \]

• \({\widehat \sigma}_{MLE}|\sigma \sim \sigma \sqrt{\chi_{n-1}^2/n}\)    (sampling distribution of unknown format). A very rough approximation based on the Central Limit Theorem (CLT) is such that \[ {\widehat \sigma}_{MLE}|\sigma \sim N(\sigma((n-1)/n)^{1/2},\sigma^2/(2n)). \]

• The CLT also helps showing that \[ MAD | \sigma \overset{a}\sim N\left[(2\sigma^2/\pi)^{1/2};\frac{\sigma^2(1-2/\pi)}{n}\right] \equiv N\left[0.7979\sigma;0.3634\sigma^2/n\right]. \] Then, an approximate 95% confidence interval for \(\sigma\) based on the MAD is given by \[ [MAD/(0.7979-0.7122509/\sqrt{n});MAD/(0.7979+0.7122509/\sqrt{n})] \]

2.1 Confidence intervals

IC.mu = c(mu.mle+qt(c(0.025,0.975),n-1)*sigma.mle/sqrt(n))
IC.sigma2 = c(n*sigma2.mle/qchisq(0.025,n-1),n*sigma2.mle/qchisq(0.975,n-1))
IC.sigma = c(sigma.mle/(sqrt((n-1)/n)+1.96/sqrt(2*n)),sigma.mle/(sqrt((n-1)/n)-1.96/sqrt(2*n)))
IC.sigma1 = c(mad/(0.7979+0.7122509/sqrt(n)),mad/(0.7979-0.7122509/sqrt(n)))

tab = rbind(IC.mu,IC.sigma2,IC.sigma,IC.sigma1)
tab

##                 [,1]      [,2]
## IC.mu     -0.1593920 0.3595120
## IC.sigma2  1.3206221 0.5934305
## IC.sigma   0.7697888 1.1498599
## IC.sigma1  0.8355240 1.0769570

3 Monte Carlo simulation exercises

set.seed(20102025)
alpha = 0.05      # Level of the confidence interval
mu     = 0           # Normal(mu,sigma^2)
sigma = 1
R        = 1000     # Replications

par(mfrow=c(3,4))
for (n in c(20,50,100)){
  estimators = matrix(0,R,4)
  for (r in 1:R){
    x = rnorm(n,mu,sigma)
    mu.mle = mean(x)
    sigma2.mle = mean((x-mu.mle)^2)
    sigma.mle = sqrt(sigma2.mle)
    mad = mean(abs(x-mu.mle))
    estimators[r,] = c(mu.mle,sigma2.mle,sigma.mle,mad)
  }
  hist(estimators[,1],xlab="mu",main="mu (MLE)",prob=TRUE)
  abline(v=mu,col=2,lwd=2)
  legend("topright",legend=paste("n = ",n,sep=""),bty="n")
  hist(estimators[,2],xlab="sigma2",main="sigma2 (MLE)",prob=TRUE)
  abline(v=sigma^2,col=2,lwd=2)
  hist(estimators[,3],xlab="sigma",main="sigma (MLE)",prob=TRUE)
  abline(v=sigma,col=2,lwd=2)
  hist(estimators[,4],xlab="sigma",main="sigma (MAD)",prob=TRUE)
  abline(v=sqrt(2/pi)*sigma,col=2,lwd=2)
}

4 Approximating sampling distributions via bootstrap

x = c(-0.027,-0.530,-1.582,-0.431,-0.248,0.478,0.697,1.484,-0.514,
          1.312,-0.584,0.854,0.473,-0.691,-0.651,-0.596,0.557,-0.950,
          -0.006,1.162,0.444,-1.156,-0.880,-1.073,-0.332,0.824,0.512,
          0.311,-1.027,-1.069,2.307,-0.126,1.903,-0.197,-0.581,0.170,
          0.236,0.481,1.246,0.341,0.664,-0.769,1.017,-0.353,-0.074,
          -0.515,-0.997,1.073,-0.110,2.526)
n = length(x)
mu0.mle = mean(x)
sigma20.mle = mean((x-mu0.mle)^2)
sigma0.mle = sqrt(sigma20.mle)
mad0 = mean(abs(x-mu0.mle))

set.seed(3141593)
R = 10000  # Bootstrap replications
estimators = matrix(0,R,2)
for (r in 1:R){
  x1 = sample(x,size=n,replace=TRUE)
  mu.mle = mean(x1)
  sigma2.mle = mean((x1-mu.mle)^2)
  estimators[r,] = c(mu.mle,sigma2.mle)
}

quantiles = apply(estimators,2,quantile,c(0.025,0.975))

par(mfrow=c(1,2))
hist(estimators[,1],xlab="mu",main="mu (MLE)",prob=TRUE)
abline(v=mu0.mle,col=2,lwd=4)
abline(v=mean(estimators[,1]),lwd=4)
segments(IC.mu[1],0,IC.mu[2],0,col=2,lwd=4)
segments(quantiles[1,1],0.1,quantiles[2,1],0.1,lwd=4)
legend("topright",legend=c("Exact","Bootstrap"),col=2:1,lwd=4,bty="n")

hist(estimators[,2],xlab="sigma2",main="sigma2 (MLE)",prob=TRUE)
abline(v=sigma20.mle,col=2,lwd=4)
abline(v=mean(estimators[,2]),lwd=4)
segments(IC.sigma2[1],0,IC.sigma2[2],0,col=2,lwd=4)
segments(quantiles[1,2],0.1,quantiles[2,2],0.1,lwd=4)

Notice that the exact confidence intervals are relatively close to the bootstrap ones, but they do not match. The intervals are almost identical for the estimator of \(\mu\). For \(\sigma^2\), the exact sampling distribution is a \(\chi^2\), which skewed, while the bootstrap approximation resembles a Gaussian distribution. The differences decrease when the sample size \(n\) increases and the bootstrap will eventually converse to the sampling distribution of the estimator.