Download Monte Carlo Method: Estimation of Expected Values and Confidence Intervals and more Study notes Mathematical Statistics in PDF only on Docsity!
Statistics 512 Notes 9: The Monte Carlo Method
Continued
The Monte Carlo method:
Consider a function g X (^^ )of a random vector X^ where X
has density f^ (^ X^ ). Consider the expected value of g X (^^ ):
E g X [ ( )] g x f ( ) ( ) x dx
Suppose we take an iid random samples X^ 1 ,^ ^ ,^ Xm from
the density f^ (^ X^ ).
Then by the law of large numbers
1
[ ( )]
m
P i i
g X
E g X
m
The Monte Carlo method is to do a simulation to draw
1
m
X X from the density f ( X ) and estimate E g X [ ( )]
by
1
[ ( )]
m
i i
g X
E g X
m
In a simulation, we can make
m as large as we want.
Standard error of the estimate is
2
1
1
ˆ [ ( )]
m
m (^) i i
i i
E g X
g X
g X
m
S
m
By the Central Limit Theorem, an approximate 95%
confidence interval for E g X [^ (^^ )] is
ˆ [ ( )]
[ ( )] 1.
E g X
E g X S
Example: Monte Carlo estimation of
Define the unit square as a square centered at (0.5,0.5) with
sides of length 1 and the unit circle as the circle centered at
the origin with a radius of length 1. The ratio of the area of
the unit circle that lies in the first quadrant to the area of the
unit square is / 4.
Let U 1^ and U^ 2 be iid uniform (0,1) random variables. Let
1 2
g U ( , U )=1 if 1 2
( U , U )is in the unit circle and 0
otherwise. Then 1 2
[ ( , )]
E g U U
Monte Carlo method: Repeat the experiment of drawing
1 2
X ( U , U ),
1
U (^) and 2
U (^) iid uniform (0,1) random
variables, m times and estimate^ ^ by
1 2 1
n
i i i
g U U
m
. An approximate 95% confidence interval for^ ^ is
2
1 2 1
1 2 1
m
m (^) i i i
i i i
g U U
g U U
m
m
Because g U (^1^ ,^ U^ 2 )=0 or 1, (1) is equivalent to
m
In R, the command runif(n) draws n iid uniform (0,1)
random variables.
Here is a function for estimating pi:
$uci
[1] 3.
Back to Example 5.8.5:
The true size of the 0.05 nominal size t-test for random
samples of size 20 contaminated normal distribution A?
We want to estimate
1 20
E I t x [ { ( , , x ) 1.729}]
Monte Carlo method:
,1 , 1
1 20
[ { ( , , ) 1.729}]
m
i i i
I t x x
E I t x x
m
where (^) ,1 ,
i i
x x is a random sample of size 20 from the
contaminated normal distribution A.
[Here 1 20
X ( X , , X ) and f ( X )is the density of a
random sample of size 20 from the contaminated normal
distribution A and g X (^^ )^ ^ I t X { (^^ 1 ,^ ^ ,^ X 20 )^ 1.729}.]
How to draw a random observation from the contaminated
normal distribution A?
(1) Draw a Bernoulli random variable B with p=0.25;
(2) If B=0, draw a random observation from the
standard normal distribution. If B=1, draw a
random observation from the normal distribution
with mean 0 and standard deviation 25.
In R, the command rnorm(n,mean=0,sd=1) draws a random
sample of size n from the normal distribution with the
specified mean and SD. The command rbinom(n,size=1,p)
draws a random sample of size n from Bernoulli
distribution with probability of success p.
R function for obtaining Monte Carlo estimate
1 20
E I t x [ { ( , , x ) 1.729}]
empalphacn=function(nsims){
Obtains the empirical level of the test discussed in
Example 5.8.
nsims is the number of simulations
sigmac=25; # SD when observation is contaminated
probcont=.25; # Probability of contamination
alpha=.05; # Significance level for t-test
n=20; # Sample size
tc=qt(1-alpha,n-1); # Critical value for t-test
ic=0; # ic will count the number of times t-test is rejected
for(i in 1:nsims){
Bernoulli random variable which determines whether
each observation in sample is from standard normal or
normal with SD sigmac
b=rbinom(n,size=1,prob=probcont);
Sample observations from standard normal when b=0 and
normal with SD sigmac when b=
samp=rnorm(n,mean=0,sd=1+b*24);
Calculate t-statistics for testing mu=0 based on sample
tstat=mean(samp)/(var(samp)^.5/n^.5);
Theorem 5.8.1: Suppose the random variable U has a
uniform (0,1) distribution. Let F be the cdf of a random
variable that is strictly increasing on some interval I, where
F=0 to the left of I and F=1 to the right of I. Then the
random variable
1
X F ( U )
has cdf F, where
1
F (0)
=left
endpoint of I and
1
F (1)
=right endpoint of I.
Proof: A uniform distribution on (0,1) has the CDF
U
F u u (^) for u (0,1). Using the fact that the CDF F is a
strictly monotone increasing function on the interval I, then
on
1
1
[ ] [ ( ) ]
= [ ( ( )) ( )]
= [ ( )]
P X x P F U x
P F F U F x
P U F x
F x
Difficult to use this method when simulating random
variables whose inverse CDF cannot be obtained in closed
form.
Other methods for simulating a random variable:
(1) Accept-Reject Algorithm (Chapter 5.8.1)
(2) Markov chain Monte Carlo Methods (Chapter 11.4)
R commands for generating random variables
runif -- uniform random variables
rbinom -- binomial random variables
rnorm -- normal random variables
rt -- t random variables
rpois -- Poisson random variables
rexp -- exponential random variables
rgamma -- gamma random variables
rbeta -- beta random variables
rcauchy -- Cauchy random variables
rchisq -- chisquared random variables
rF -- F random variables
rgeom -- geometric random variables
rnbinom -- negative binomial random variables
Bootstrap Procedures
Bootstrap standard errors
1
n
X X iid with CDF F and variance
2
.
2
1
2 1
n
n
X X
Var Var X X
n n n
^
SD X ( )
n
We estimate SD X ( )by (^ )^
s
SE X
n
(^) where s (^) is the
sample standard deviation.
What about 1
n
SD Median X X? This SD depends in
a complicated way on the distribution F of the X’s. How to
approximate it?
Real World: F^ ^ X^ 1 ,^ ^ ,^ X^ n ^ Tn^ Median X (^^ 1 ,^ ^ ,^ Xn ).
- Draw
1
n
X X.
- Compute
1
n n
T g X X.
- Repeat steps 1 and 2 m times to get
,1 ,
n n m
T T
- Let
2
, (^1 1) ,
1 m 1 m
boot n i i r (^) n r
se T T
m m
The bootstrap involves two approximations:
not so small approx. error small approx. error
ˆ
n
F n (^) F n boot
SD T SD T se
R function for bootstrap estimate of SE(Median)
bootstrapmedianfunc=function(X,bootreps){
medianX=median(X);
vector that will store the bootstrapped medians
bootmedians=rep(0,bootreps);
for(i in 1:bootreps){
Draw a sample of size n from X with replacement and
calculate median of sample
Xstar=sample(X,size=length(X),replace=TRUE);
bootmedians[i]=median(Xstar);
seboot=var(bootmedians)^.5;
list(medianX=medianX,seboot=seboot);
Example: In a study of the natural variability of rainfall, the
rainfall of summer storms was measured by a network of
rain gauges in southern Illinois for the year 1960.
>rainfall=c(.02,.01,.05,.21,.003,.45,.001,.01,2.13,.07,.01,.
> median(rainfall)
[1] 0.
> bootstrapmedianfunc(rainfall,10000)
$medianX
[1] 0.
$seboot
