Code for Monte Carlo integration of a normal PDF

MCinte.R

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)