require(latticeExtra)
require(mosaic)
require(Hmisc)
require(kableExtra)

Compute the definite integral from 0 to pi/3 of sin(t)dt

The theoretical value should be equal to (1-cos(pi/3)) = 0.5

    n <- 10000
    upper <- pi/3
    lower <- 0
    x <- runif(n, lower, upper)
    theta.hat <- (upper-lower)*(sum(sin(x)))*(1/n)

    m <- cbind(round(theta.hat,7),1 - cos(pi/3))
    m <- rbind(c("MC estimate","Theoretical Integral"),m)
    kbl(m, caption = "Definite Integral of sin(t)") %>%
  kable_styling(bootstrap_options = "striped",
                full_width = F)
Definite Integral of sin(t)
MC estimate Theoretical Integral
0.4947344 0.5

The approximated value is very nearly equal to the theoretical value
   
 



Compute a Monte Carlo estimate of the standard normal cdf, by generating from the Uniform(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)

The CDF for each element in the set x is determined by dividing the mean of a sample drawn from a uniform distribution on (0,x). The 95% confidence interval is then determined for each value of x, from the variance of the sample.

    x <- seq(0.1, 2.5, length = 25)
    n <- 10000
    cdf <- numeric(length(x))
    Variance <- numeric(length(x))
    for (i in 1:length(x)) {
        u <- runif(n,0,x[i])
        g <- x[i] * exp(-(u^2 / 2))
        Variance[i] <- (1/(2*pi))*(x[i]^2/n)*mean((g-mean(g))^2)
        cdf[i] <- mean(g) / sqrt(2 * pi) + 0.5
    }

    Phi <- pnorm(x)
    t_value <- qt(0.05,n-1,lower.tail=FALSE)
    upper.range.95CI <- cdf+(t_value*(sqrt(Variance)))
    lower.range.95CI <- cdf-(t_value*(sqrt(Variance)))
    m <- (round(rbind(x, cdf, Phi,upper.range.95CI,lower.range.95CI), 3))
    kbl(m, caption = "Comparison of MC estimate of standard normal CDF and values from R's pnorm") %>%
  kable_styling(bootstrap_options = "striped",
                full_width = F)
Comparison of MC estimate of standard normal CDF and values from R's pnorm
x 0.10 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900 1.000 1.100 1.200 1.300 1.400 1.500 1.600 1.700 1.800 1.900 2.000 2.100 2.200 2.300 2.400 2.500
cdf 0.54 0.579 0.618 0.655 0.691 0.726 0.758 0.788 0.816 0.842 0.865 0.885 0.903 0.920 0.935 0.947 0.955 0.967 0.971 0.980 0.985 0.986 0.995 0.991 0.999
Phi 0.54 0.579 0.618 0.655 0.691 0.726 0.758 0.788 0.816 0.841 0.864 0.885 0.903 0.919 0.933 0.945 0.955 0.964 0.971 0.977 0.982 0.986 0.989 0.992 0.994
upper.range.95CI 0.54 0.579 0.618 0.656 0.692 0.726 0.758 0.788 0.816 0.843 0.866 0.887 0.905 0.923 0.938 0.951 0.960 0.972 0.977 0.987 0.993 0.995 1.006 1.003 1.012
lower.range.95CI 0.54 0.579 0.618 0.655 0.691 0.726 0.758 0.788 0.815 0.841 0.864 0.884 0.901 0.918 0.932 0.943 0.950 0.961 0.964 0.972 0.976 0.976 0.984 0.978 0.985 
    s <- rbind(2,0.976,0.977,0.984,0.969)
    s <- cbind(c("x","cdf","Phi","lower bound of 95%CI","upper bound of 95%CI"),s)
    kbl(s, caption = "95 percent CI for x = 2") %>%
  kable_styling(bootstrap_options = "striped",
                full_width = F)
95 percent CI for x = 2
x 2
cdf 0.976
Phi 0.977
lower bound of 95%CI 0.984
upper bound of 95%CI 0.969








Compute a Monte Carlo estimate theta.hat of theta = the integral from 0 to 0.5 of e^-x dx by sampling from Uniform(0, 0.5), and estimate the variance of theta.hat. Find another Monte Carlo estimator theta.star by sampling from the exponential distribution. Which of the variances is smaller, and why?

    n <- 10000
    upper <- 0.5
    lower <- 0
    x <- runif(n, lower, upper)
    gx <- (exp(-x))
    theta.hat <- (upper-lower)*(mean(gx))
    var.theta.hat <- mean((gx - mean(gx))^2)/(4*n)

    y <- rexp(n)
    gy <- (y < 0.5)
    theta.star <- mean(gy)
    var.theta.star <- mean((gy - mean(gy))^2)/n

    m <- rbind(c((round(theta.hat,4)),(round(theta.star,4)),round(1 - exp(-0.5),4),var.theta.hat,var.theta.star))
    m <- rbind(c("theta.hat","theta.star","Theoretical Value","theta hat variance","theta star variance"),m)
    kbl(m, caption = "Comparison of estimate variances") %>%
    kable_styling(bootstrap_options = "striped",
                full_width = F)
Comparison of estimate variances
theta.hat theta.star Theoretical Value theta hat variance theta star variance
0.3943 0.3911 0.3935 3.21710856322738e-07 2.3814079e-05

Both distribution samplings provide good estimates for the theoretical value (1-e^(0.5)) but the variance of theta hat is significantly smaller than the variance of theta star. This aligns with one's intuition, as the variance of a uniform distribution on the range (0,0.5) is 1/48, smaller than the variance of an exponential distribution with mean 1. It's then intuitive that the transformed distributions might maintain this relationship