AA037014 Statistical Computing

8: Monte Carlo Methods in Inference

Dr. Diana Abdul Wahab
Introduction
· Monte Carlo methods may refer to any method in statistical inference or
numerical analysis where simulation is used.
- However, in this lecture only a subset of these methods are discussed.
- This lecture introduces some of the Monte Carlo methods for
statistical inference.
· Monte Carlo methods can be applied to estimate parameters of the
sampling distribution of a statistic, mean squared error (MSE), percentiles,
or other quantities of interest.
· Monte Carlo studies can be designed to assess the coverage probability for
confidence intervals, to find an empirical Type I error rate of a test
procedure, to estimate the power of a test, and to compare the
performance of different procedures for a given problem.

2/49
Introduction
· In statistical inference there is uncertainty in an estimate.
· The methods covered in this chapter use repeated sampling from a given
probability model, sometimes called parametric bootstrap, to investigate
this uncertainty.
· If we can simulate the stochastic process that generated our data,
repeatedly drawing samples under identical conditions, then ultimately we
hope to have a close replica of the process itself reflected in the samples.
· Other Monte Carlo methods, such as (nonparametric) bootstrap, are based
on resampling from an observed sample.
· Resampling methods are covered in the next lecture.

3/49
Content
1. Monte Carlo methods for estimation:
· standard error
· MSE
· confidence interval

2. Monte Carlo methods for hypothesis test:

· Type I error
· power of a test
· power comparison

3. Application

4/49
Monte Carlo Methods for Estimation
Suppose X1, . . . , Xn is a random sample from the distribution of X. An
estimator θ^ for a parameter θ is an n variate function

^ = θ
θ ^(X , . . . , X )
1 n

of the sample.

Functions of the estimator θ^ are therefore n-variate functions of the data,

also. Random variates from the sampling distribution of θ^ can be generated
by repeatedly drawing independent random samples x(j) and computing
(j) (j) (j)
^
θ ^(x
= θ
1
, . . . , xn ) for each sample.

5/49
Monte Carlo methods for estimation: standard
error
Suppose that X1 , X2 are iid from a standard normal distribution. Estimate
the mean difference E|X1 − X2|.

To obtain MC estimate of θ based on m replicates: 1. Generate random

samples of size 2 from the standard normal distribution. 2. Compute the
(j)
replicates of ^
θ . 3. Compute the mean of the replicates.

6/49
Monte Carlo methods for estimation: standard
error
For j=1000,

m <- 1000
g <- numeric(m)
for (i in 1:m) {
x <- rnorm(2)
g[i] <- abs(x[1] - x[2])
}
mean(g)

## [1] 1.106031

7/49
Monte Carlo methods for estimation: standard
error
Estimating the standard error of the mean
−−−−− −
^) = √V ar( ^)
se(θ θ

sqrt(sum((g - mean(g))^2)) / m

## [1] 0.02641957

8/49
Monte Carlo estimation of MSE
^ − θ) 2 ]
^) = E[(θ
M SE(θ

If m random samples are generated from the distribution of X, then a Monte

Carlo estimate of the MSE is
m
1 (j)
2
M^
SE = ^
∑ (θ − θ)
m
j=1

9/49
Example: estimating the MSE of a trimmed
mean
· A trimmed mean is sometimes applied to estimate the center of a
continuous symmetric distribution that is not necessarily normal.
· In this example, we compute an estimate of the MSE of a trimmed mean at
k-th level:

n−k
1
X̄[−k] = ∑ X(i)
n − 2k
i=k+1

Goal: obtain a Monte Carlo estimate of the M SE(X̄[−1] ) of the first level
trimmed mean assuming that the sampled distribution is standard normal.

10/49
Example: estimating the MSE of a trimmed
mean
Steps:
(j) (j)
1. Generate x1 , . . . , xn idd from distribution of X.

2. Sort x1 in increasing order.

(j) (j)
, . . . , xn

3. Compute T (j) .
1 (j)
= ∑x
n−2 (i)

4. Compute M SE(T
¯
).

11/49
Example: estimating the MSE of a trimmed
mean
Exercise: compute the MSE of a trimmed mean where n=20 with m=1000
replications.

n <- 20
m <- 1000
tmean <- numeric(m)
for (i in 1:m) {
x <- sort(rnorm(n))
tmean[i] <- sum(x[2:(n-1)]) / (n-2)
}
mse <- mean(tmean^2)

12/49
Example: estimating the MSE of a trimmed
mean

mse

## [1] 0.05236238

The standard error is

sqrt(sum((tmean - mean(tmean))^2)) / m

## [1] 0.007235998

The estimate of MSE for the trimmed mean in this run is approximately
0.0523624 with se 0.007236. For comparison, the MSE of the sample mean X
is Var(X)/n, which is 1/20 = 0.05 in this example.

13/49
Example: estimating the MSE of a median
Another example on the median:

n <- 20
m <- 1000
tmean <- numeric(m)
for (i in 1:m) {
x <- sort(rnorm(n))
tmean[i] <- median(x)
}
mse <- mean(tmean^2)

14/49
Example: estimating the MSE of a median

mse

## [1] 0.07557061

sqrt(sum((tmean - mean(tmean))^2)) / m #se

## [1] 0.008688686

The estimate of MSE for the sample median is approximately 0.0755706 and
se(MSE) = 0.0086887.

15/49
Estimating a confidence level
· One type of problem that arises frequently in statistical applications is the
need to evaluate the cdf of the sampling distribution of a statistic, when
the density function of the statistic is unknown or intractable.
· For example, many commonly used estimation procedures are derived
under the assumption that the sampled population is normally
distributed.
· In practice, it is often the case that the population is non-normal and in
such cases, the true distribution of the estimator may be unknown or
intractable.
· The following examples illustrate a Monte Carlo method to assess the
confidence level in an estimation procedure.

16/49
Estimating a confidence level
· If (U,V) is a confidence interval estimate for an unknown parameter θ, then
U and V are statistics with distributions that depend on the distribution FX
of the sampled population X.
· The confidence level is the probability that the interval (U,V) covers the true
value of the parameter θ.
· Evaluating the confidence level is therefore an integration problem.

17/49
Estimating a confidence level
· Note that the sample-mean Monte Carlo approaches to evaluating an
integral ∫ g(x)dx do not require that the function g(x) is specified.
- It is only necessary that the sample from the distribution g(X) can be
generated.
· It is often the case in statistical applications, that g(x) is in fact not
specified, but the variable g(X) can be easily generated.

18/49
Estimating a confidence level
· We consider the confidence interval estimation procedure for variance.
· This procedure is sensitive to mild departures from normality.
· We use Monte Carlo methods to estimate the true confidence level when
the normal theory confidence interval for variance is applied to non-
normal data.
· The classical procedure based on the assumption of normality is outlined
first.

19/49
Estimating a confidence level
If X1 , . . . , Xn is a random sample from a Normal(μ, σ 2 ) distribution, n>=2,
and S 2 is the sample variance, then
2
(n − 1)S 2
V = ∼ χ (n − 1)
2
σ

A one side 100(1−α)% confidence interval is given by (0, (n−1)S2/χ2 α ), where

χ α is the α-quantile of the χ (n − 1) distribution. If the sampled population
2 2

is normal with variance σ 2 , then the probability that the confidence interval
contains σ 2 is 1 − α.

20/49
Estimating a confidence level
The calculation of the 95% upper confidence limit (UCL) for a random sample
size n = 20 from a Normal(0, σ 2 = 4) distribution is shown below.

n <- 20
alpha <- .05
x <- rnorm(n, mean=0, sd=2)
UCL <- (n-1) * var(x) / qchisq(alpha, df=n-1)
UCL

## [1] 6.579472

21/49
Estimating a confidence level
Several run tests…

x <- rnorm(100, mean=0, sd=2)

(UCL <- (n-1) * var(x) / qchisq(alpha, df=n-1))

## [1] 6.696469

x <- rnorm(1000, mean=0, sd=2)

(UCL <- (n-1) * var(x) / qchisq(alpha, df=n-1))

## [1] 7.66146

x <- rnorm(10000, mean=0, sd=2)

(UCL <- (n-1) * var(x) / qchisq(alpha, df=n-1))

## [1] 7.45693

22/49
Estimating a confidence level
Approximately 95% of the intervals based should contain σ 2 .

· Empirical confidence level is an estimate of the confidence level obtained

by simulation.
· For the simulation experiment, repeat the steps above a large number of
times, and compute the proportion of intervals that contain the target
parameter.

23/49
Monte Carlo experiment to estimate a
confidence level
Steps:

1. Generate j random samples.

2. Compute the confidence interval Cj
3. Compute the proportion of intervals that contain the target parameter
I (θ ∈ Cj )

4. Compute the empirical confidence level

24/49
Monte Carlo experiment to estimate a
confidence level
For μ = 0, σ = 2, n = 20, m = 1000, α = 0.05 :

n <- 20
alpha <- .05
UCL <- replicate(1000, expr = {
x <- rnorm(n, mean = 0, sd = 2)
(n-1) * var(x) / qchisq(alpha, df = n-1)
} )

25/49
Monte Carlo experiment to estimate a
confidence level

# count the number of intervals that contain sigma^2=4

sum(UCL > 4)

## [1] 949

#or compute the mean to get the confidence level

mean(UCL > 4)

## [1] 0.949

26/49
Monte Carlo experiment to estimate a
confidence level
The result is that 0.949 intervals satisfied (UCL > 4), so the empirical
confidence level is 94.9% in this experiment. The result will vary but should
be close to the theoretical value, 95%. The standard error of the estimate is
(0.95(1 − 0.95)/1000)1/2 = 0.00689.

27/49
Monte Carlo experiment to estimate a
confidence level
The interval estimation procedure for estimating variance is sensitive to
departures from normality, so the true confidence level may be different
than the stated confidence level when data are non-normal.

28/49
Monte Carlo Methods for Hypothesis Tests
Suppose that we wish to test a hypothesis concerning a parameter θ that lies
in a parameter space Θ. The hypotheses of interest are

H0 : θ ∈ Θ0

vs

H1 : θ ∈ Θ1

where Θ0 and Θ1 partition the parameter space Θ.

29/49
Monte Carlo Methods for Hypothesis Tests
Two types of error can occur in statistical hypothesis testing.

1. A Type I error occurs if the null hypothesis is rejected when in fact the null
hypothesis is true.
2. A Type II error occurs if the null hypothesis is not rejected when in fact the
null hypothesis is false.

30/49
Monte Carlo Methods for Hypothesis Tests
· The significance level of a test is denoted by α, and α is an upper bound
on the probability of Type I error.
· The probability of rejecting the null hypothesis depends on the true value
of θ.

31/49
Monte Carlo Methods for Hypothesis Tests
· The probability of Type I error is the conditional probability that the null
hypothesis is rejected given that H0 is true.
· Thus, if the test procedure is replicated a large number of times under the
conditions of the null hypothesis, the observed Type I error rate should be
at most (approximately) α.
· If T is the test statistic and T* is the observed value of the test statistic,
then T* is significant if the test decision based on T* is to reject H0.
· The significance probability or p-value is the smallest possible value of α
such that the observed test statistic would be significant.

32/49
Empirical Type I error rate
· An empirical Type I error rate can be computed by a Monte Carlo
experiment.
· The test procedure is replicated a large number of times under the
conditions of the null hypothesis.
· The empirical Type I error rate for the Monte Carlo experiment is the
sample proportion of significant test statistics among the replicates.

33/49
Monte Carlo experiment to assess Type I error
rate:
Steps:

1. Generate j random samples

2. Compute the test statistics Tj
3. Record the test decision Ij = 1 if H0 is rejected, otherwise 0.
4. Compute the proportion of significant tests.

This proportion is the observed Type I error rate.

34/49
Monte Carlo experiment to assess Type I error
rate:
Suppose that X1 , . . . , X20 is a random sample from a N (μ, σ 2 )
distribution. Test H0 : μ = 500 H1 : μ > 500 at α = 0.05. Under the null
hypothesis

X̄ − 500
T∗ = −− ∼ t(19)
S/√20

Large values of T* support the alternative hypothesis.

35/49
Monte Carlo experiment to assess Type I error
rate:
Use a Monte Carlo method to compute an empirical probability of Type I
error when σ = 100, and check that it is approximately equal to α = 0.05.

36/49
Monte Carlo experiment to assess Type I error
rate:

n <- 20
alpha <- .05
mu0 <- 500
sigma <- 100
m <- 10000
#number of replicates
p <- numeric(m)
#storage for p-values
for (j in 1:m) {
x <- rnorm(n, mu0, sigma)
ttest <- t.test(x, alternative = "greater", mu = mu0)
p[j] <- ttest$p.value
}

37/49
Monte Carlo experiment to assess Type I error
rate:

p.hat <- mean(p < alpha)

se.hat <- sqrt(p.hat * (1 - p.hat) / m)
print(c(p.hat, se.hat))

## [1] 0.050400000 0.002187689

38/49
Monte Carlo experiment to assess Type I error
rate:
· The observed Type I error rate in this simulation is 0.0504, and the
standard error of the estimate is approximately sqrt(0.05 x 0.95/m) =
0.0022.
· Estimates of Type I error probability will vary, but should be close to the
nominal rate α = 0.05 because all samples were generated under the null
hypothesis from the assumed model for a t-test (normal distribution).
· In this experiment the empirical Type I error rate differs from α = 0.05 by
less than one standard error.

39/49
Monte Carlo experiment to assess Type I error
rate:
· Theoretically, the probability of rejecting the null hypothesis when μ = 500
is exactly α = 0.05 in this example.
· The simulation really only investigates empirically whether the method of
computing the p-value in t.test (a numerical algorithm) is consistent with
the theoretical value α = 0.05.

40/49
Power of a test
· When we test two hypotheses, H0 and H1, a Type II error happens when
H1 is actually true, but we mistakenly fail to reject H0.
· The power of a test is a measure of how likely we are to correctly reject H0
when it’s false, and it’s represented by a function called the power
function.
· The power function tells us the probability of rejecting H0 when the true
value of the parameter is a specific value, denoted as θ.
· For a given θ1 , the probability of making a Type II error (failing to reject H0
when it’s false) for that value is 1 minus the power of the test at θ1 , or in
other words, 1 − π(θ1 ).

41/49
Power of a test
· Ideally, we want a test that has a low probability of error.
· The probability of making a Type I error (rejecting H0 when it’s actually
true) is controlled by something called the significance level, denoted as α.
· However, we are more interested in minimizing Type II errors, which
means having high power under the alternative hypothesis (H1).
· When comparing different test procedures for the same hypotheses at the
same significance level, we focus on comparing their power.
· A test with higher power is preferable because it has a lower probability of
making a Type II error.

42/49
Power of a test
Use simulation to estimate power and plot an empirical power curve for the
t-test.

· To plot the curve, we need the empirical power for a sequence of

alternatives θ along the horizontal axis.
· Each point corresponds to a Monte Carlo experiment.
· The outer for loop varies the points θ (mu) and the inner replicate loop
estimates the power at the current θ.

43/49
Power of a test

n <- 20
m <- 1000
mu0 <- 500
sigma <- 100
mu <- c(seq(450, 650, 10))
#alternatives
M <- length(mu)
power <- numeric(M)

44/49
Power of a test

for (i in 1:M) {
mu1 <- mu[i]
pvalues <- replicate(m, expr = {
#simulate under alternative mu1
x <- rnorm(n, mean = mu1, sd = sigma)
ttest <- t.test(x,
alternative = "greater", mu = mu0)
ttest$p.value
} )
power[i] <- mean(pvalues <= .05)
}

45/49
Power of a test
The estimated power π
^(θ) values are now stored in the vector power.

Next, plot the empirical power curve, adding vertical error bars at
π
^(θ) ± se^ (π^(θ)) using the errbar function in the Hmisc package.

46/49
Power of a test

library(Hmisc) #for errbar

plot(mu, power)
abline(v = mu0, lty = 1)
abline(h = .05, lty = 1)
#add standard errors
se <- sqrt(power * (1-power) / m)
errbar(mu, power, yplus = power+se, yminus = power-se,
xlab = bquote(theta))
lines(mu, power, lty=3)

47/49
Power of a test

48/49
Power of a test
Note that the empirical power π^(θ) is small when θ is close to θ0 = 500, and

increasing as θ moves farther away from θ0 , approaching 1 as θ -> ∞.

49/49