public
Created

Code for Monte Carlo integration of a normal PDF

  • Download Gist
MCinte.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
x <- seq(-4,4,length=1000)
y <- dnorm(x,0,1)
plot(x,y,type='l',ylim=c(0,.5))
###Set the lower limits of interest
ll <- -1
ul <- 1
######Hit and miss MC integration#####
###First pick domain
x.dom <- seq(-4,4,length=1000)
y.dom <- seq(0,.5,length=1000)
prog <- txtProgressBar(min=0, max=n.sim, char="*", style=3)
n.sim <- 10000
x.r <- y.r <- vector()
under.curve <- rep(0,n.sim)
for(i in 1:n.sim){
###Draw random x and y
x.r[i] <- sample(x.dom,1)
y.r[i] <- sample(y.dom,1)
if(y.r[i] < dnorm(x.r[i],0,1)){
###Set further constraints on the x axis
if(x.r[i] >= ll && x.r[i] <= ul){under.curve[i]<-1}}
setTxtProgressBar(prog, i)
}
 
plot(-100,-100,ylim=range(y.dom),xlim=range(x.dom),col=2,lwd=2)
points(x.r,y.r,col=under.curve+1,pch=19)
lines(x,y,col=2,lwd=4)
 
f.h <- sum(under.curve)/n.sim
V.samp <- diff(range(x.dom))*diff(range(y.dom))
###Approximate integration value is area
area <- f.h*V.samp
###compare to the true value
pnorm(ul,0,1) - pnorm(ll,0,1)

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.