## R code for Section 5.8: 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

## Example 5.8.1 (Estimation of Pi), code in Appendix B ##
## Pi=3.141592653589793...

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

piest <- 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=se)   # ouput
}
piest(100)                          # n=100


## Example 5.8.3 (Monte Carlo Integration) ##

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

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



## Source: ##
# http://www.stat.wmich.edu/mckean/HMC/Rcode/AppendixB/ #
# http://www.math.uic.edu/~jyang06/stat411/rcode/5_8.R  #