# Hit-or-miss Monte Carlo integration code # Define function ex1.fcn <- function(x){ h <- 4/(1 + x^2) # return is a more formal alternative than simply listing h return(h) } # Plot function and include a shaded polygon x=seq(0,1,.05) hx=ex1.fcn(x) xx=c(x,1,0) hxx=c(hx,0,0) plot(xx,hxx,type="n",ylim=c(0,4),xlab="x",ylab="h(x) ",main="h(x)") polygon(xx,hxx,col="pink",border="red") # Set parameters my.n <- 10000 my.a <- 0; my.b <- 1 # Generate random numbers uniformly inside your box my.rand.x <- runif(my.n, min=my.a, max=my.b) my.rand.y <- runif(my.n, min=0, max=4) # Here is the entire sample points(my.rand.x,my.rand.y,pch=19,cex=.2) # Start over again for a more interesting graph dev.off() x=seq(0,1,.05) plot(x,4/(1+x^2),type="l",col="red",ylab=expression(h(x)),ylim=c(0,4)) # Use an outline rather than a polygon this time lines(c(0,0),c(0,4),col="red") lines(c(1,1),c(0,2),col="red") lines(c(0,1),c(0,0),col="red") # Figure out points under the curve under= my.rand.y < ex1.fcn(my.rand.x) # Plot points both under and over the curve #These steps could be written in one command points(my.rand.x[under],my.rand.y[under],col="red",cex=.2) points(my.rand.x[!under],my.rand.y[!under],col="black",cex=.2) lines(c(0,1),c(4,4),col="black") lines(c(1,1),c(2,4),col="black") title("Hit or Miss simulation") legend(.1,.5,pch=c(19,19),c("Below h(x)","Above h(x)"),pt.cex=.2,col=c("red","black"),xjust=0,yjust=0,cex=1,bg="white") # An alternative plotting approach hm.df=data.frame(my.rand.x,my.rand.y,under) with(hm.df,plot(my.rand.x,my.rand.y,cex=.2,pch=as.numeric(under+1),col=as.numeric(under+1),ylab="h(x)",xlab="x")) title("Hit or Miss simulation") legend(.1,.5,pch=c(1,2),c("Below h(x)","Above h(x)"),pt.cex=.2,col=c(1,2),xjust=0,yjust=0,cex=1,bg="white") lines(x,4/(1+x^2),col="red") lines(c(0,0),c(0,4),col="red") lines(c(1,1),c(0,2),col="red") lines(c(0,1),c(0,0),col="red") lines(c(0,1),c(4,4),col="black") lines(c(1,1),c(2,4),col="black") # Output proportion of points under the curve; final estimate of pi my.m <- sum(my.rand.y < ex1.fcn(my.rand.x)) my.prop=my.m/my.n I.approx.HM <- 4*(my.b - my.a)*(my.m/my.n) print(c(I.approx.HM,my.prop)) # R code to integrate h() using classical Monte Carlo method # Example 1—we will reuse the random sample and parameters from the Hit-or-Miss simulation I.approx.MC <- ((my.b - my.a)/my.n)*sum(ex1.fcn(my.rand.x)) print(c(I.approx.HM,I.approx.MC))