# Monte Carlo Methods

# the hit or miss method


g <- function(x) -2*x * log(x)/log(x+2)

a <- 0
b <- 1
c <- 1

x <- seq(0.001,1,by=0.001)
gx <- g(x)
plot(gx~x,type = "l",ylim = c(0,c))

N <- 1000

X <- runif(N,a,b)
Y <- runif(N,0,c)
points(X,Y,col = ifelse(Y <= g(X),"blue","red"))

I <- c*(b-a)*sum(Y <= g(X))/N
I


# Wolfram alpha gives 0.5685

# The classical Monte Carlo approach to approximating an integral

N <- 10000
a <- 0
b <- 1
U <- runif(N,a,b)

I <- (b-a) * mean(g(U))
I


# Monte Carlo optimization -- minimizing a function

g <- function(x) x**2 - 3*sin(pi*x) + 2 * cos(3*pi*x)

x <- seq(-3,3,length=201)
gx <- g(x)
plot(gx~x,type="l")

a <- -3
b <- 3
N <- 500
U <- runif(N,a,b)

abline(v=U,lwd=.5, col = rgb(.5,0,0,.1))

j <- which.min(g(U)) # give the index of the minimum entry
xhat <- U[j]

abline(v = xhat,col="chartreuse", lwd= 2)
xhat



## What about a bivariate function??

h <- function(x,y){
  
  a <- (x*sin(20*y) + y*sin(20*x))^2*cosh(sin(10*x)*x)
  b <- (x*cos(10*y) - y*cos(10*x))^2*cosh(cos(20*y)*y)
  c <- sqrt(x^2 + y^2)/10
  val <- a + b + c
  return(val)  
  
}

x <- seq(-1,1,length=101)
y <- x
z <- outer(x,y,FUN=h)

par(mar=c(1.1,1.1,1.1,1.1))
persp(x,y,z,      
      theta = 24, 
      phi = 30,
      cex.axis = 0.8,
      expand = 0.5,
      border = rgb(.545,0,0),
      box = F,
      lwd = .8)

# try to find the minimizer using the uniform sampling method

N <- 10000
# range of X
a1 <- -1
b1 <- 1
X <- runif(N,a1,b1) # generate X coordinates

a2 <- -1
b2 <- 1
Y <- runif(N,a2,b2) # generate Y coordinates

hXY <- numeric(N)
for(k in 1:N){ # evaluate the function at all of these coordinates

  hXY[k] <- h(X[k],Y[k])  
  
}

j <- which.min(hXY)

minimizer <- c(X[j],Y[j])
minimizer


par(mar=c(.1,.1,.1,.1))
contour(x,y,z,
        nlevels = 100,
        col="black", 
        lwd = .5,
        xaxt = "n",
        bty = "n",
        yaxt = "n")
points(X,Y,col = rgb(0,0,.545,.1))
points(X[j],Y[j],col="red",pch = 19)
abline(h = Y[j],col="red")
abline(v = X[j],col="red")