## R code for Section 4.8: The Method of Monte Carlo ## I) Simulate random numbers from uniform distribution x = runif(1000, min=-1, max=1) hist(x, freq=F) curve(dunif(x, min=-1, max=1), lwd=2, lty=2, add=T) # Application: Estimate int_2^4 ln(x) dx 4*log(4)-2*log(2)-2 #truth: 2.158883 n = 100000 x = runif(n, min=2, max=4) gx = 2*log(x) est = mean(gx) # estimate of integral se = sd(gx)/sqrt(n) # estimate of standard error c(est-1.96*se, est+1.96*se) # 95% confidence interval ## II) 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 ## Alternative approach sample(x=c(1,2,3,4,5,6), size=100, replace=T) ## III) Use inverse cdf to generate random numbers ## Example: Pdf f(x) = 2*x, 01 # 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