## R code: The Method of Monte Carlo ##

## Simulate the outcomes of tossing a fair die
#  Obtain a random sequence of {1,2,3,4,5,6}
#  n is the sample size

dietossing <- function(n) 
{
  ceiling(runif(n)*6)
}
dietossing(100)          # run the function with n=100

## Estimation of Pi ##
## Pi=3.1415926...

# Obtain the estimate of pi and its standard error
# n is the number of simulations

pi.est <- function(n) 
{
  u1 <- runif(n)               # uniform random sample of size n
  u2 <- runif(n)
  cnt <- rep(0,n)              # repeat 0 for n times
  chk <- u1^2 + u2^2 - 1       # vector operation
  cnt[chk < 0] <- 1            # if chk[i]<0, then cnt[i]=1
  est <- mean(cnt)             # average
  se <- 4*sqrt(est*(1-est)/n)
  est <- 4*est
  list(estimate=est, standard_error=se)   # ouput
}
pi.est(1000)                          # n=100


## Monte Carlo Integration ##

# Obtain the estimate of pi and its standard error
# n is the number of simulations

pi.est2 <- function(n) 
{
  samp <- 4*sqrt(1-runif(n)^2)
  est <- mean(samp)
  se <- sqrt(var(samp)/n)
  list(est=est,se=se)
}
pi.est2(1000)


## Generate a random sample from expoential distribution ##

n <- 200
Usamp <- runif(n)          # uniform random sample of size n  
Y <- -log(1 - Usamp)       # inverse of cdf of expoential function
hist(Y, freq=F)

densityf <- function(x) { exp(-x) }   # exponential density function
x <- seq(0, 10, by=0.05)
points(x, densityf(x), type="l")  

# http://www.stat.wmich.edu/mckean/HMC/Rcode/AppendixB/ #