# 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)
}

mcint.fcn <- function(nloop,my.n){
# Set parameters
my.a <- 0; my.b <- 1
I.approx.HM=c()
I.approx.MC=c()

for(i in 1:nloop){
# 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)

# Figure out points under the curve
under= my.rand.y < ex1.fcn(my.rand.x)

# Output proportion of points under the curve; final estimate of pi
my.m <- sum(under)
I.approx.HM[i] <- 4*(my.b - my.a)*(my.m/my.n)

I.approx.MC[i] <- ((my.b - my.a)/my.n)*sum(ex1.fcn(my.rand.x))
}

boxplot(I.approx.HM,I.approx.MC,xlab="")
axis(1,at=c(1,2),labels=c("Hit or Miss MC","Classic MC"))

}

mcint.fcn(100,1000)