Outlines

• Introduction
• Examples
  • Compare estimators for mean
  • Compare OLS and lasso estimators
  • Hypothesis testing, size
  • Hypothesis testing, power
  • Hypothesis testing, Confidence interval
  • Hypothesis testing, Confidence interval, finite sample sizes
• Further simulations

What is a simulation study?

Simulation:

• A numerical technique for conducting experiments on the computer.
• It involves random sampling from probability distributions.

Rationale for simulation study

• After we propose a statistical method, we want to validate the method using simulation so people can use it with confidence.
• Exact analytical derivations of properties are rarely possible.
• Large sample approximations to properties are often possible.
  • But what happens in the finite sample size case?
  • And what happens when assumptions are violated?

Common questions simulation study can answer

• Is an estimator consistent with infinite samples?
  • If so, Is an estimator biased in finite samples?
  • Is it still consistent under departures from assumptions?
• How does an estimator compare to its competing estimators on the basis of bias, mean square error, efficiency, etc.?
• Does a procedure for constructing a confidence interval for a parameter achieve the advertised nominal level of coverage?
• Does a hypothesis testing procedure attain the nominal (advertised) level of significance or size?
• How to estimate power from simulation study?

Monte Carlo simulation approximation

• An test statistic or estimator has a true sampling distribution.
• Ideally, we could want to know this true sampling distribution in order to address the issues on the previous slide.
• But derivation of the true sampling distribution is not tractable.
• Hence we approximate the sampling distribution of an estimator or a test statistic under a particular set of conditions.

How to approximate (Draw a picture to illustrate)

Typical Monte Carlo simulation involves the following procedure:

1. Generate \(B\) independent datasets \((1, 2, \ldots, B)\) under conditions of interest.
2. Compute the numerical value of an estimator/test statistic \(T\) from each dataset. Then we have \(T_1, T_2, \ldots, T_b, \ldots, T_B\).
3. If \(B\) is large enough, summary statistics across \(T_1, \ldots, T_B\) should be good approximations to the true properties of the estimator/test statistic.

Example:

Consider an estimator for a mean parameter \(\theta\):

• \(T_b\) is the value of \(T\) from the \(b^{th}\) data, \((b = 1, \ldots, B)\).
• The mean of all \(T_b\)’s is an estimate of the true mean.

We will have more concrete examples later.

Example 1, Compare estimator (normal distribution)

• Assume \(X \sim N(\mu, \sigma^2)\) with \(\mu = 1\) and \(\sigma^2 = 1\).
• We want to compare three estimators for the mean \(\mu\) of this distribution based on a random sample \(Y_1, Y_2, \ldots, Y_n\).
  • Sample Mean \(T^{mean}\)
  • Sample Median \(T^{median}\)
  • Sample 20% trimmed mean \(T^{mean20\%}\)

Note: 20% trimmed mean indicates average value after deleting first 10% quantile and last 10% quantile.

E.g. 3,2,4,7,1,6,5,8,9,10

• mean = mean(1:10)
• mean20% = mean(2:9)

Evaluation criteria

k: mean, median, mean20%

• \(T^{(k)}\): estimator of interest
• Bias: \(E({T}^{(k)}) - \mu\)
• variance: \(E(({T}^{(k)} - \mu)^2)\)
• Mean squared error: \(MSE = Bias^2 + vairance\)
  • MSE is a commonly used criteria for evaluating estimators.
• Relative efficiency
  • For any consistent estimator \(T^{(1)}\), \(T^{(2)}\), (i.e. \(E(T^{(1)}) = E(T^{(2)}) = \mu\))
  • \(RE = \frac{MSE(T^{(1)})}{MSE(T^{(2)})}\)
  • is the relative efficiency of estimator 2 to estimator 1.
  • \(RE < 1\) means estimator 1 is preferred.

\(T^{mean}\)

\[T^{mean} = \frac{1}{n}\sum_{i=1}^n Y_i\]

• Bias: \(E({T}^{(mean)}) - \mu = \frac{1}{n}\sum_{i=1}^n \mu - \mu = 0\)
• Variance: \(E(({T}^{(mean)} - \mu)^2) = \frac{\sigma^2}{n}\)
• MSE: \(MSE = Bias^2 + vairance = \frac{\sigma^2}{n}\)

Can be calculated analytically

\(T^{median}\)

\[T^{median} = Median_{i=1}^n Y_i\]

• Bias: \(E({T}^{(median)}) - \mu =0?\)
• Variance: \(E(({T}^{(median)} - \mu)^2) =?\)
• MSE: \(MSE = Bias^2 + vairance =?\)

Variance is hard to be calculated analytically

\(T^{mean20\%}\)

Can be calculated

\[T^{mean20\%} = \frac{1}{0.8n}\sum_{i=0.1n+1}^{0.9n} Y_i\]

• Bias: \(E({T}^{(mean20\%)}) - \mu =0?\)
• Variance: \(E(({T}^{(mean20\%)} - \mu)^2) =?\)
• MSE: \(MSE = Bias^2 + vairance =?\)

Variance is hard to be calculated analytically

Evaluation via simulation

Since sometimes it is hard to evaluate the estimator analytically, we can evaluate use simulation studies.

• algorithm:
  1. Sample \(n\) subject from \(f\) (i.e. \(N(\mu, \sigma^2)\))
  2. Calculate \({T}^{(k)}\)
  3. Repeat Step 1,2 B times to obtain \(T_1^{(k)}, \ldots, T_B^{(k)}\)
  4. Evaluate the property of \({T}^{(k)}\) using the Monte Carlo samples \(T_1^{(k)}, \ldots, T_B^{(k)}\).
• \(\bar{T}^{(k)} = \sum_{b=1}^B T_b^{(k)}\)
• \(\widehat{Bias} = \bar{T}^{(k)} - \mu\)
• \(\widehat{vairance} = \frac{1}{B - 1} \sum_{b=1}^B (T_b^{(k)} - \bar{T}^{(k)})\)
• Mean squared error: \(\widehat{MSE} = \widehat{Bias}^2 + \widehat{vairance}\)
• Relative efficiency
  • For any consistent estimator \(T^{(1)}\), \(T^{(2)}\), (i.e. \(E(T^{(1)}) = E(T^{(2)}) = \mu\))
  • \(\widehat{RE} = \frac{\widehat{MSE}(T^{(1)})}{\widehat{MSE}(T^{(2)})}\)
  • is the relative efficiency of estimator 2 to estimator 1.
  • \(\widehat{RE} < 1\) means estimator 1 is preferred.

Simulation experiment

mu <- 1; B <- 1000; n <- 20

T1 <- numeric(B); T2 <- numeric(B); T3 <- numeric(B)

for(b in 1:B){
  set.seed(b)
  x <- rnorm(n, mean = mu)
  
  T1[b] <- mean(x)
  T2[b] <- median(x)
  T3[b] <- mean(x, trim = 0.1)
}

bias_T1 <- abs(mean(T1) - mu); bias_T2 <- abs(mean(T2) - mu); bias_T3 <- abs(mean(T3) - mu)
bias <- c(bias_T1, bias_T2, bias_T3)

var_T1 <- var(T1); var_T2 <- var(T2); var_T3 <- var(T3)
vars <- c(var_T1, var_T2, var_T3)

MSE1 <- bias_T1^2 + var_T1; MSE2 <- bias_T2^2 + var_T2; MSE3 <- bias_T3^2 + var_T3
MSE <- c(MSE1, MSE2, MSE3)

RE <- MSE1/MSE

comparisonTable <- rbind(bias, vars, MSE, RE)
colnames(comparisonTable) <- c("mean", "median", "mean 20% trim")
comparisonTable

##              mean      median mean 20% trim
## bias 3.832736e-05 0.006673072   0.000688909
## vars 5.174830e-02 0.074680145   0.054550469
## MSE  5.174830e-02 0.074724675   0.054550944
## RE   1.000000e+00 0.692519633   0.948623449

Example 2, Compare estimator (exponential distribution)

• Assume \(X \sim Exp(\mu)\) with \(\mu = 1\). \[f(x) = \exp(-x), x\ge0\]
• We want to compare three estimators for the mean \(\mu\) of this distribution based on a random sample \(Y_1, Y_2, \ldots, Y_n\).
  • Sample Mean \(T^{mean}\)
  • Sample Median \(T^{median}\)
  • Sample 20% trimmed mean \(T^{mean20\%}\)

Simulation experiment

mu <- 1; B <- 1000; n <- 20

T1 <- numeric(B); T2 <- numeric(B); T3 <- numeric(B)

for(b in 1:B){
  set.seed(b)
  x <- rexp(n, rate = mu)
  
  T1[b] <- mean(x)
  T2[b] <- median(x)
  T3[b] <- mean(x, trim = 0.1)
}

bias_T1 <- abs(mean(T1) - mu); bias_T2 <- abs(mean(T2) - mu); bias_T3 <- abs(mean(T3) - mu)
bias <- c(bias_T1, bias_T2, bias_T3)

var_T1 <- var(T1); var_T2 <- var(T2); var_T3 <- var(T3)
vars <- c(var_T1, var_T2, var_T3)

MSE1 <- bias_T1^2 + var_T1; MSE2 <- bias_T2^2 + var_T2; MSE3 <- bias_T3^2 + var_T3
MSE <- c(MSE1, MSE2, MSE3)

comparisonTable <- rbind(bias, vars, MSE)
colnames(comparisonTable) <- c("mean", "median", "mean 20% trim")
comparisonTable

##            mean     median mean 20% trim
## bias 0.01040018 0.26539760    0.13243138
## vars 0.05021121 0.05052217    0.04247023
## MSE  0.05031938 0.12095806    0.06000830

Example 3, Compare estimator (Cauchy distribution)

• Assume \(X \sim Cauchy(\mu,\gamma)\) with \(\mu = 0\) and \(\gamma = 1\). \[f(x; \mu, \gamma) = \frac{1}{\pi (1 + x^2)}\]

• We want to compare three estimators for the mean \(\mu\) of this distribution based on a random sample \(Y_1, Y_2, \ldots, Y_n\).
  • Sample Mean \(T^{mean}\)
  • Sample Median \(T^{median}\)
  • Sample 20% trimmed mean \(T^{mean20\%}\)

Simulation experiment

mu <- 0; B <- 1000; n <- 20

T1 <- numeric(B); T2 <- numeric(B); T3 <- numeric(B)

for(b in 1:B){
  set.seed(b)
  x <- rcauchy(n, location  = mu)
  
  T1[b] <- mean(x)
  T2[b] <- median(x)
  T3[b] <- mean(x, trim = 0.1)
}

bias_T1 <- abs(mean(T1) - mu); bias_T2 <- abs(mean(T2) - mu); bias_T3 <- abs(mean(T3) - mu)
bias <- c(bias_T1, bias_T2, bias_T3)

var_T1 <- var(T1); var_T2 <- var(T2); var_T3 <- var(T3)
vars <- c(var_T1, var_T2, var_T3)

MSE1 <- bias_T1^2 + var_T1; MSE2 <- bias_T2^2 + var_T2; MSE3 <- bias_T3^2 + var_T3
MSE <- c(MSE1, MSE2, MSE3)

comparisonTable <- rbind(bias, vars, MSE)
colnames(comparisonTable) <- c("mean", "median", "mean 20% trim")
comparisonTable

##             mean      median mean 20% trim
## bias   0.1187952 0.005811866    0.00296929
## vars 742.5037597 0.146730729    0.37387966
## MSE  742.5178721 0.146764506    0.37388847

UMVUE

• Uniformly minimum-variance unbiased estimator (UMVUE) is an unbiased estimator that has lower variance than any other unbiased estimator for all possible values of the parameter.

• Blue: the best linear unbiased estimator (BLUE) of the coefficients is given by the ordinary least squares (OLS) estimator,

\[\hat{\beta}^{ols} = \arg \min_\beta \|y - X\beta\|^2_2\]

• Modern statistics care about both bias and variance of an estimator instead of only unbiased estimator.

• Lasso estimator

\[\hat{\beta}^{lasso} = \arg \min_\beta \|y - X\beta\|^2_2 + \lambda \|\beta\|_1\]

Example 4, evaluate bias, variance and MSE of OLS, lasso

Model setup:

• \(\beta_0 = 1\), \(\beta_1 = 3\), \(\beta_2 = 5\), \(\beta_3 = \beta_4 = \beta_5 = 0\)
• \(X_1, \ldots, X_5 \sim N(0,1)\)
• \(Y_0 = \beta_0 + X_1 \beta_1 + X_2 \beta_2 +\ldots, X_5 \beta_5\)
• \(Y = Y_0 + N(0,3^2)\)
• \(n = 30\)

Goal:

• Wants to evaluate bias, variance and MSE of OLS, lasso estimators
• \(bias(\beta, \hat{\beta})^2 = \sum_{j=0}^5 bias(\beta_j, \hat{\beta}_j)^2\)
• \(var_T(\hat{\beta}) = \sum_{j=0}^5 var(\hat{\beta}_j)\)
• \(MSE(\hat{\beta}) = bias(\beta, \hat{\beta})^2 + var_T(\hat{\beta})\)

For simulation of OLS (ordinary least square)

• Set \(B = 1000\)
• For \(b (1 \le b \le B)\)
  1. simulate \(Y_b = Y_0 + N(0,3^2)\)
  2. Calculate \(\hat{\beta}^{ols}_b\)
  3. Calculate \(\hat{\bar{\beta}}^{ols} = \frac{1}{B}\sum_{b=1}^B \hat{\beta}^{ols}_b\)
  4. Calculate bias \(bias(\beta, \hat{\bar{\beta}}^{ols})^2 = \sum_{j=0}^5 bias(\beta_j, \hat{\bar{\beta}}^{ols}_j)^2\)
  5. Calculate variance: \(var(\hat{\beta}) = \frac{1}{B-1} (\hat{\beta}^{ols}_b - \hat{\bar{\beta}}^{ols})^2\) , \(var_T(\hat{\beta}) = \sum_{j=0}^5 var(\hat{\beta}_j)\)
  6. \(MSE(\hat{\beta}) = bias(\beta, \hat{\beta})^2 + var_T(\hat{\beta})\)

For simulation of lasso

• Set \(B = 1000\)
• For \(b (1 \le b \le B)\), given a specific tuning parameter,
  1. simulate \(Y_b = Y_0 + N(0,3^2)\)
  2. Calculate \(\hat{\beta}^{lasso}_b\)
  3. Calculate \(\hat{\bar{\beta}}^{lasso} = \frac{1}{B}\sum_{b=1}^B \hat{\beta}^{lasso}_b\)
  4. Calculate bias \(bias(\beta, \hat{\bar{\beta}}^{lasso})^2 = \sum_{j=0}^5 bias(\beta_j, \hat{\bar{\beta}}^{lasso}_j)^2\)
  5. Calculate variance: \(var(\hat{\beta}) = \frac{1}{B-1} (\hat{\beta}^{lasso}_b - \hat{\bar{\beta}}^{lasso})^2\) , \(var_T(\hat{\beta}) = \sum_{j=0}^5 var(\hat{\beta}_j)\)
  6. \(MSE(\hat{\beta}) = bias(\beta, \hat{\beta})^2 + var_T(\hat{\beta})\)
• Repeat this procedure for all tuning parameters

Simulation experiment

library(lars)

## Loaded lars 1.2

beta0 <- 1; beta1 <- 3; beta2 <- 5
beta3 <- beta4 <- beta5 <- 0
beta <- c(beta0,beta1,beta2,beta3,beta4,beta5)
p <- 5; n <- 30
set.seed(32611)
X0 <- matrix(rnorm(p*n),nrow=n)
X <- cbind(1, X0) ## X fixed

Y0 <- X %*% as.matrix(beta, ncol=1)

B <- 1000
beta_ols <- replicate(n = B, expr = list())
beta_lasso <- replicate(n = B, expr = list())

for(b in 1:B){
  set.seed(b)

  error <- rnorm(length(Y0), sd = 3)
  Y <- Y0 + error
  
  abeta_ols <- lm(Y~X0)$coefficients
  beta_ols[[b]] <- abeta_ols
  
  alars <- lars(x = X, y = Y, intercept=F)
  beta_lasso[[b]] <- alars
}

Simulation experiment (continue 2)

beta_ols_hat <- Reduce("+",beta_ols)/length(beta_ols)
beta_ols_bias <- sqrt(sum((beta_ols_hat - beta)^2))
beta_ols_variance <-  sum(Reduce("+", lapply(beta_ols,function(x) (x - beta_ols_hat)^2))/(length(beta_ols) - 1))
beta_ols_MSE <- beta_ols_bias + beta_ols_variance

as <- seq(0,10,0.1)
lassoCoef <- lapply(beta_lasso, function(alars){
  alassoCoef <- t(coef(alars, s=as, mode="lambda"))
  alassoCoef
} )

beta_lasso_hat <- Reduce("+",lassoCoef)/length(lassoCoef)
beta_lasso_bias <- sqrt(colSums((beta_lasso_hat - beta)^2))
beta_lasso_variance <-  colSums(Reduce("+", lapply(lassoCoef,function(x) (x - beta_lasso_hat)^2))/(length(beta_ols) - 1))
beta_lasso_MSE <- beta_lasso_bias + beta_lasso_variance

Simulation experiment (continue 3)

par(mfrow=c(2,2))
plot(x = as, y = beta_lasso_bias, type = "l", xlab="lambda", main="bias")
abline(h = beta_ols_bias, lty=2)
legend("topleft", legend = c("lasso", "OLS"), col = 1, lty = c(1,2))

plot(x = as, y = beta_lasso_variance, col=4, type = "l", xlab="lambda", main="variance")
abline(h = beta_ols_variance, lty=2, col=4)
legend("bottomleft", legend = c("lasso", "OLS"), col = 4, lty = c(1,2))

plot(x = as, y = beta_lasso_MSE, col=2, type = "l", xlab="lambda", main="MSE")
abline(h = beta_ols_MSE, lty=2, col=2)
legend("topleft", legend = c("lasso", "OLS"), col = 2, lty = c(1,2))


plot(x = as, y = beta_lasso_MSE, ylim = c(min(beta_lasso_bias), max(beta_lasso_MSE)), 
     type="n", xlab="lambda", main="MSE = bias^2 + variance")
lines(x = as, y = beta_lasso_MSE, type = "l", col=2)
lines(x = as, y = beta_lasso_variance, col=4)
lines(x = as, y = beta_lasso_bias, col=1)
legend("bottomright", legend = c("bias", "variance", "MSE"), col = c(1,4,2), lty = c(1))

Hypothesis testing (one-sided)

Hypothesis testing (two-sided)

Does a hypothesis testing procedure attain the nominal (advertised) level of significance or size?

Consider the following tests:

• Parametric: t.test()
• Non-Parametric: wilcox.test()

Simulation for sizes of hypothesis tests

• Test \(H_0: \mu = \mu_0\) against \(H_A: \mu \neq \mu_0\)
• To evaluate whether size/level of test achieves advertised \(\alpha\)
  • Generate data under \(H_0: \mu = \mu_0\) and calculate the proportion of rejections of \(H_0\)
• Approximation the true probability of rejecting \(H_0\) when it is true
  • Proportion should be approximately equal to \(\alpha\)

Consider using the following simulation setting

• n = 20
• B = 1000
• alpha = 0.05
• two.sided test

Example 5, t.test() size

alpha <- 0.05; mu=0; n <- 20; B <- 1000

counts <- 0
for(b in 1:B){
  set.seed(b)
  
  x <- rnorm(n, mean=mu)
  atest <- t.test(x)
  
  if(atest$p.value < alpha){
    counts <- counts + 1
  }
}

counts/B

## [1] 0.05

Example 6, wilcox.test() size

alpha <- 0.05; mu=0; n <- 20; B <- 1000

counts <- 0
for(b in 1:B){
  set.seed(b)
  
  x <- rnorm(n, mean=mu)
  atest <- wilcox.test(x)
  
  if(atest$p.value < alpha){
    counts <- counts + 1
  }
}

counts/B

## [1] 0.047

Example 7, t.test() only use first 10 samples size

alpha <- 0.05; mu=0; n <- 20; B <- 1000

counts <- 0
for(b in 1:B){
  set.seed(b)
  
  x <- rnorm(n, mean=mu)[1:10]
  atest <- t.test(x)
  
  if(atest$p.value < alpha){
    counts <- counts + 1
  }
}

counts/B

## [1] 0.056

Compare these test procedures

• If several tests has the same size, how to compare which test is better? - Compare power
• To evaluate power
  • Generate data under some alternatives, \(\mu \neq \mu_0\) (i.e. \(\mu = \mu_1\)) and calculate the proportion of rejections of \(H_0\)
• Approximate the true probability of rejecting \(H_0\) when the alternative is true.

Example 8, t.test() power

alpha <- 0.05; mu=0; mu1 <- 1; n <- 20; B <- 1000

counts <- 0
for(b in 1:B){
  set.seed(b)
  
  x <- rnorm(n, mean=mu1)
  atest <- t.test(x)
  
  if(atest$p.value < alpha){
    counts <- counts + 1
  }
}
counts / B

## [1] 0.985

Example 8, Validate the result using R function

library(pwr)
alpha <- 0.05; mu=0; mu1 <- 1; n <- 20; B <- 1000
pwr.t.test(n = n, d = mu1, sig.level = alpha, type = "one.sample", alternative = "two.sided")

## 
##      One-sample t test power calculation 
## 
##               n = 20
##               d = 1
##       sig.level = 0.05
##           power = 0.9885913
##     alternative = two.sided

Example 9, wilcox.test() power

alpha <- 0.05; mu=0; mu1 <- 1; n <- 20; B <- 1000

counts <- 0
for(b in 1:B){
  set.seed(b)
  
  x <- rnorm(n, mean=mu1)
  atest <- wilcox.test(x)
  
  if(atest$p.value < alpha){
    counts <- counts + 1
  }
}
counts / B

## [1] 0.979

Example 10, t.test() only use first 10 samples power

alpha <- 0.05; mu=0; mu1 <- 1; n <- 20; B <- 1000

counts <- 0
for(b in 1:B){
  set.seed(b)
  
  x <- rnorm(n, mean=mu1)[1:10]
  atest <- t.test(x)
  
  if(atest$p.value < alpha){
    counts <- counts + 1
  }
}
counts / B

## [1] 0.805

Example 10, , Validate the result using R function

alpha <- 0.05; mu=0; mu1 <- 1; n <- 20; B <- 1000
pwr.t.test(n = 10, d = mu1, sig.level = alpha, type = "one.sample", alternative = "two.sided")

## 
##      One-sample t test power calculation 
## 
##               n = 10
##               d = 1
##       sig.level = 0.05
##           power = 0.8030969
##     alternative = two.sided

Confidence interval

A \(1 - \alpha\) confidence interval for parameter \(\mu\) is an interval \(C_n = (a,b)\) where \(a = a(X_1, \ldots, X_n)\) and \(b = b(X_1, \ldots, X_n)\) are functions of the data such that \[P(\mu \in C_n) \ge 1 - \alpha\] In other words, \((a,b)\) traps \(\mu\) with probability \(1 - \alpha\). We call \(1 - \alpha\) the coverage of the confidence interval.

Interpretation:

• A confidence interval is not a probability statement about \(\mu\) since \(\mu\) is a fixed quantity, not a random variable.
• if I repeat the experiment over and over, the interval will contain the parameter \(\mu\) 95 percent of the time.

Example 11, Confidence interval (\(\mu\) can be random)

• \(\mu = 1\)
• \(n=10, i = 1, \ldots, n\)
• \(X_i \sim N(\mu, 1)\)
• \(95\%\) confidence interval: \((\bar{X} - qt(0.975, n-1) \times \frac{\hat{\sigma}}{\sqrt{n}}, \bar{X} + qt(0.975, n-1) \times \frac{\hat{\sigma}}{\sqrt{n}})\)

mu <- 1
B <- 1000
n <- 10

x <- rnorm(n, mean=mu)
c(mean(x) - qt(0.975, n-1) * sd(x) / sqrt(n), mean(x) + qt(0.975, n-1) * sd(x) / sqrt(n))

## [1] 0.400803 2.219041

t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 3.2595, df = 9, p-value = 0.009847
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.400803 2.219041
## sample estimates:
## mean of x 
##  1.309922

counts <- 0
for(b in 1:B){
  set.seed(b)

  x <- rnorm(n, mean=mu)
  atest <- t.test(x)
  
  if(atest$conf.int[1] < mu & mu < atest$conf.int[2]){
    counts <- counts + 1
  }
}

counts/B

## [1] 0.944

Example 12, Confidence interval (using normal approximation, n=10)

Is the Confidence interval accurate in finite samples?

• \(\mu = 1\)
• \(n=10, i = 1, \ldots, N\)
• \(X_i \sim N(\mu, 1)\)
• \(95\%\) confidence interval: \((\bar{X} - 1.96 \times \frac{\hat{\sigma}}{\sqrt{n}}, \bar{X} + 1.96 \times \frac{\hat{\sigma}}{\sqrt{n}})\)

set.seed(32611)
B <- 1000 ## number of test
counts <- 0 ## number of accepted tests
n <- 10
mu <- 1
alpha <- 0.05
for(b in 1:B){
  x <- rnorm(n, mean=mu, sd = 1)
  lb <- mean(x) - qnorm(alpha/2, lower.tail = F) * sd(x)/sqrt(n)
  ub <- mean(x) + qnorm(alpha/2, lower.tail = F) * sd(x)/sqrt(n)
  if(lb < mu & ub > mu){
    counts <- counts + 1
  }
}
print(counts/B)

## [1] 0.909

Example 12, Confidence interval (using normal approximation, n=100)

• behavior in large sample

• \(\mu = 1\)
• \(n=100, i = 1, \ldots, N\)
• \(X_i \sim N(\mu, 1)\)

• \(95\%\) confidence interval: \((\bar{X} - 1.96 \times \frac{\hat{\sigma}}{\sqrt{n}}, \bar{X} + 1.96 \times \frac{\hat{\sigma}}{\sqrt{n}})\)

set.seed(32611)
B <- 1000 ## number of test
counts <- 0 ## number of accepted tests
n <- 100
mu <- 1
alpha <- 0.05
for(b in 1:B){
  x <- rnorm(n, mean=mu, sd = 1)
  lb <- mean(x) - qnorm(alpha/2, lower.tail = F) * sd(x)/sqrt(n)
  ub <- mean(x) + qnorm(alpha/2, lower.tail = F) * sd(x)/sqrt(n)
  if(lb < mu & ub > mu){
    counts <- counts + 1
  }
}
print(counts/B)

## [1] 0.945

Principle for simulation study

1. A Monte Carlo simulation is just like any other experiment
   • Factors that are of interest to vary in the experiment: sample size n, distribution of the data, magnitude of variation…
   • Each combination of factors is a separate simulation, so that many factors can lead to very large number of combinations and thus number of simulations (time consuming)
   • Results must be recorded and saved in a systematic way
   • Don’t only choose factors favorable to a method you have developed!
   • Sample size B (number of data sets) must deliver acceptable precision.
2. Save everything!
   • Save the individual estimates in a file and then analyze
   • as opposed to computing these summaries and saving only them
   • Critical if the simulation takes a long time to run!
3. Keep S small at first
   • Test and refine code until you are sure everything is working correctly before carrying out final production runs
   • Get an idea of how long it takes to process one data set
4. Set a different seed for each run and keep records
   • Ensure simulation runs are independent
   • Runs may be replicated if necessary
5. Document your code

Further simulation topics

• model assumption is violated
• central limit theorem
• False discovery rate
• machine learning
• validate a distribution as a mixture of other distributions.
• compare ward test, score test and likelihood ratio test

Reference

• https://www.stat.nus.edu.sg/~stazjt/teaching/ST2137/lecture/lec%2011.pdf
• http://www4.stat.ncsu.edu/~davidian/st810a/simulation_handout.pdf

knitr::purl("simulation.rmd", output = "simulation.R ", documentation = 2)

## 
## 
## processing file: simulation.rmd

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |..                                                               |   3%
  |                                                                       
  |....                                                             |   5%
  |                                                                       
  |.....                                                            |   8%
  |                                                                       
  |.......                                                          |  11%
  |                                                                       
  |.........                                                        |  14%
  |                                                                       
  |...........                                                      |  16%
  |                                                                       
  |............                                                     |  19%
  |                                                                       
  |..............                                                   |  22%
  |                                                                       
  |................                                                 |  24%
  |                                                                       
  |..................                                               |  27%
  |                                                                       
  |...................                                              |  30%
  |                                                                       
  |.....................                                            |  32%
  |                                                                       
  |.......................                                          |  35%
  |                                                                       
  |.........................                                        |  38%
  |                                                                       
  |..........................                                       |  41%
  |                                                                       
  |............................                                     |  43%
  |                                                                       
  |..............................                                   |  46%
  |                                                                       
  |................................                                 |  49%
  |                                                                       
  |.................................                                |  51%
  |                                                                       
  |...................................                              |  54%
  |                                                                       
  |.....................................                            |  57%
  |                                                                       
  |.......................................                          |  59%
  |                                                                       
  |........................................                         |  62%
  |                                                                       
  |..........................................                       |  65%
  |                                                                       
  |............................................                     |  68%
  |                                                                       
  |..............................................                   |  70%
  |                                                                       
  |...............................................                  |  73%
  |                                                                       
  |.................................................                |  76%
  |                                                                       
  |...................................................              |  78%
  |                                                                       
  |.....................................................            |  81%
  |                                                                       
  |......................................................           |  84%
  |                                                                       
  |........................................................         |  86%
  |                                                                       
  |..........................................................       |  89%
  |                                                                       
  |............................................................     |  92%
  |                                                                       
  |.............................................................    |  95%
  |                                                                       
  |...............................................................  |  97%
  |                                                                       
  |.................................................................| 100%

## output file: simulation.R

## [1] "simulation.R "