Demo monte carlo integration technique in R.

Monte Carlo integration.r

boundingRange <- function(x){
  max(x) - min(x)
}

runif_box <-function(n,x){
  b = boundingRange(x)
  u = runif(n)
  u*b+min(x)
}

xrange<-function(x1, x2){
  seq(x1,x2, length.out=1e4)
}

mc.integration <- function(f, x1, x2, n=1e4, plot.fig=F){
  if(x1>x2){
    t1=x2; t2=x1
    x1=t1; x2=t2
  }
  xv=xrange(x1,x2)
  x.random = runif_box(n, xv)
  y.random = runif_box(n, f(xv)*1.2 )
  under = (y.random < f(x.random) )&(y.random>0)
  over  = (y.random > f(x.random) )&(y.random<0) 
  prop = (sum(under)-sum(over) )/n
  totArea = (x2-x1)*boundingRange(f(xv))
  integ = (prop)*totArea
  
  if(plot.fig){
    plot(xv, f(xv)
         ,xlim=c(x1,x2)
         ,ylim=c(min(y.random),max(y.random))
         , type='n'
    )
    points(x.random, y.random , col=under|over)
    lines(xv, f(xv)
          ,xlim=c(x1,x2)
          ,ylim=c(min(y.random),max(y.random))
          , type='l'
          ,col="blue")
    
  }  
  
  list(integ, x.random, y.random, under|over, prop, totArea)  
}

################################

f<-function(x){sin(x)*x}
xv=xrange(-pi, pi)


results = mc.integration(f,0, 10*pi, 1e5, plot.fig=T)

results[[1]]
integrate(f,0, 10*pi)