## R code for Section 5.8: The Method of Monte Carlo

## simulate the outcomes of tossing a fair die
# function to obtain a random sequence of {1,2,3,4,5,6}
# n is the length
dietossing <- function(n) {
  ceiling(runif(n)*6)
}
dietossing(100)          # run the function with n=100


## Example 5.8.1 (Estimation of Pi)
## Pi=3.141592653589793...
# function to 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)
}
piest(100)               # run the function with n=100

## Example 5.8.3 (Monte Carlo Integration)
# function to 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)
}

## Exercise 5.8.17 (Accept-Reject Algorithm)
# function to generate random sample from f(x)=beta*x^{beta-1}, 0<x<1, beta>1
# beta is the parameter; n is the sample size
samplef <- function(beta, n) {
  samp <- rep(0,n)
  num <- 1
  while(num <= n) {
    y <- runif(1)
    u <- runif(1)
    if(u <= y^(beta-1)) {
      samp[num] <- y
      num <- num + 1
    }
  }
  samp
}
b <- 3              # beta=3
samplef(b,100)

# test the function
Y <- samplef(b,1000)
densityf <- function(x) { b*x^(b-1); }
x <- seq(0, 1, by=0.01)
hist(Y, freq=F)
points(x, densityf(x), type="l")  # add density function to the histogram