Monte Carlo Simulation Example (Manual and Package Based)

MCManualPlot.png

MonteCarloSimulationExamples.R

# Monte Carlo Simulation Examples

library(MonteCarlo)

#Monte Carlo Simulation using only native tools to estimate the value of pi 
# No package neccessary. See below for
#Monte Carlo Simulation Using Package

#NATIVE LIBRARIES

#Get the hypotenuse from assuming origin
getVectorLength <- function(x,y){
  return(sqrt(x^2 + y^2))
}

# r is the radius of the circle. 
inCircle <- function(x,y,r) {
  l <- getVectorLength(x,y)
  if(l >= r){
    return(F)
  } else return(T)
}

# Make sure functions are working
if(!inCircle(1,0,1) || !inCircle(0,0,1) || inCircle(2,0,1) || !inCircle(-1,0,1) || !inCircle(0,-1,1) || inCircle(1,1,1) ){
  print("Error. In Circle Function Not Working")
}

# Manual MC Simiulation 
n <- 1000000

#Generate a uniform distribution between -1 and 1
d <- cbind(runif(n,min=-1, max=1), runif(n,min=-1, max=1)) #there might be a better way to generate a random dist but this should work
labeledplot <- cbind(d, apply(d, 1, function(x){inCircle(x[1], x[2], 1)}))
g <- head(labeledplot, n = 5000) # show 5000 points in plot
plot(g[,1], g[,2], col = (g[,3] + 1), asp=1)

# get the probability
incir <- nrow(labeledplot[labeledplot[,3] == T,])
outcir <- nrow(labeledplot[labeledplot[,3] == F,])
piest <- (incir / ( incir + outcir)) * 4
print(paste("Pi Estimation is " , piest))
print(paste("Error is " , (1 - pi / piest)))

# USING PACKAGE 

# Why use the package? Because it handles the loop generation, parrellizes the computation 
# on a number of CPU units
# define parameter grid:

# 2 Ways to Do this. I've done both. 

# Way one initializes a vector with uniformely distibuted 
# random points producing a estimated distibution. The pacakge runs in simulation 
# and the estimations are averaged. 

# Second way simply provides result as T / F 
# Iterates over to form a list of T/F. Calculate the probablity by 
# counting the T and the F.

#need to rewrite the inCircle Function because the MonteCarlo Package
#this implementation generates a list and finds the probability. Mean probability can be calculated after. 
#Requires a list return 
inCirclePkgImplementation <- function(n) {
  xvec <- runif(n,min=-1, max=1) 
  yvec <- runif(n,min=-1, max=1)
  d <- cbind(xvec, yvec) 
  labeledplot <- cbind(d, apply(d, 1, function(x){inCircle(x[1], x[2], 1)}))
  incir <- nrow(labeledplot[labeledplot[,3] == T,])
  outcir <- nrow(labeledplot[labeledplot[,3] == F,])
  piest <- (incir / ( incir + outcir)) * 4
  return(list("decision"=piest))
}

# collect parameter grids in list:
param_list=list("n"=1)
MC_result<-MonteCarlo(func=inCirclePkgImplementation, nrep=1000, param_list=param_list)
print(paste("Pi Estimation is " , mean(MC_result$results$decision))) #Taking the average of the 1000 trails
print(paste("Error is " , (1 - pi / mean(MC_result$results$decision))))


#Requires a list return
inCirclePkgImplementation2 <- function(l) {
  res <- inCircle(runif(n,min=-1, max=1), runif(n,min=-1, max=1), 1)
  return(list("decision"=res))
}

param_list=list("l" = 1)
MC_result<-MonteCarlo(func=inCirclePkgImplementation2, nrep=100000, param_list=param_list)

incount <- sum(MC_result$results$decision[1, ] == 1)
outcunt <- sum(MC_result$results$decision[1, ] == 0)
piest <- (incount / ( incount + outcunt)) * 4

print(paste("Pi Estimation is " , piest)) #Taking the average of the 1000 trails
print(paste("Error is " , (1 - pi / piest)))