Exercise 5.1

You can create a function to evaluate the MC estimate of the integral \(\int_{0}^{\pi/3} sin(t) dt\) and compare your value with the exact integral value, i.e., 0.5. Also, you should be able to demonstrate that by taking larger samples we reduce the error of our estimation and thus have a better estimate to the true integral value.

ex1 <- function(n){
u <- runif(n, min = 0, max = pi/3)
x <- mean(sin(u))
sol <- x*(pi/3)
return(sol)
}

#Compute Estimate
ex1(100000)
## [1] 0.499892

Exercise 5.2

Refer to Example 5.3. Compute a Monte Carlo estimate of the standard normal cdf, by generating from the Unifrom(0, x) distribution. Compare your estimates with the normal cdf function pnorm. Compute an estimate of the variance of your Monte Carlo estimate of \(\Phi(2)\), and a 95% confidence interval for \(\Phi(2)\).

mcZtable<-function(n){
  x<-seq(-2.5,2.5,.01)   #Z-scores
  u<-runif(n)
  Zscores<-data.frame("X"=x)
  Zscores$NormCDF<-apply(Zscores,1,function(r){
    if (r<0) {
  calc1<-mean(-r['X']*exp(-1*r['X']^2*u^2/2))
  theta.hat<-.5-calc1/sqrt(2*pi)
} else {
  calc1<-mean(r['X']*exp(-r['X']^2*u^2/2))
  theta.hat<-calc1/sqrt(2*pi)+.5
}
    return(theta.hat)
  })
  Zscores$PNorm<-apply(Zscores['X'],1,pnorm)
  View(Zscores)
}

We estimate the pdf of Z=2 ( i.e., object name “estimate”" in R,is our point estimate). We find the variance using the formula for the estimate of variance of MC estimate of integra;/( i.e., object name “v,” in R) when n= 1000. And we find a 95% Confidence Interval for Z= 2 ( i.e., object names LB and UB in R).

Also, as n → ∞ our MC estimator follows a normal distribution we can also create a large sample 95% confidence interval, which takes the form, point estimate +-Margin of error.

EstVar<-function(x,n){
  u<-runif(n)
  if (x<0) {
  samples<-.5-(-x*exp(-x^2*u^2/2))/sqrt(2*pi)
  estimate<-mean(samples)
} else {
  samples<-x*exp(-x^2*u^2/2)/sqrt(2*pi)+.5
  estimate<-mean(samples)
}
  v<-var(samples)/n
  #confint(samples,level=0.95)
  LB<-estimate-1.96*sqrt(v)
  UB<-estimate+1.96*sqrt(v)
  cat("We estimate the pdf of Z=",x,"to be",estimate,"with variance",v,"when n=",n,".\n And a 95% Confidence Interval for Z=",x," is (",LB,",",UB,").")
}

EstVar(2,1000)
## We estimate the pdf of Z= 2 to be 0.9704412 with variance 5.453455e-05 when n= 1000 .
##  And a 95% Confidence Interval for Z= 2  is ( 0.9559671 , 0.9849153 ).

5.3

Compute a Monte Carlo estimate \(\hat{\theta}\) of \(θ = \int_{0}^{0.5}e^{-x}dx\) by sampling from Uniform(0,.5), and estimate the variance of \(\hat{\theta}\) . Find another Monte Carlo estimate \(θ^{*}\) by sampling from the exponential distribution. Which of the variances is smaller, and why?

First we calculate \(\hat{\theta}\) by sampling from Uniform(0,.5)

#Function to generate thetahat and variance, n is number of iterations
thetahat <- function(n){
u <- runif(n, min = 0, max = .5)
that <- exp(-u)*.5
#varthat <- (1/n)*(mean((that - mean(that))^2))
varthat <- var(that)
#Print Output
print(paste0("Theta hat is ", mean(that), ", and the variance is ", varthat, "."))
}
thetahat(10000)
## [1] "Theta hat is 0.393211646608344, and the variance is 0.00322852126561783."

Now we calculate \(\theta^{*}\) by sampling from Exp(1).

#Function to generate thetastar, n is the number of iterations
thetastar <- function(n){
x <- rexp(n,1)
g <- (exp(-x)*(x>0)*(x<.5))/(exp(-x))
thetastar <- mean(g)
vartstar <- var(g)
#Print Results
print(paste0("Theta star is ", thetastar, ", and the variance is ", vartstar, "."))
}
thetastar(10000)
## [1] "Theta star is 0.3952, and the variance is 0.239040864086409."

For the same n, the variance of \({\hat{\theta}}\) is significantly smaller than the variance of \({\theta^{*}}\). Can you think of a reason why?