Problem 1 Part A

#Problem 1 Part A
n = 20000 #20,000
g.function <- function(x) {
    val = 1/(1+x^2)
    return(val*exp(-x))
}

unifs = runif(n, min=0, max=1)

theta.hat = mean(g.function(unifs))


#Approximate the standard error
SE.theta.hat = sd(g.function(unifs))/sqrt(n)

results.list = list(theta.hat = theta.hat, SE = SE.theta.hat)
results.list
## $theta.hat
## [1] 0.5274981
## 
## $SE
## [1] 0.001735493

Problem 1 Part B and C

#Problem 1 Part b
half.n = n/2
u = runif(half.n)
theta.anti = (1/n)*(sum(g.function(u) + g.function(1-u)))
SE.theta.anti = sd(g.function(u) + g.function(1-u))*(1/sqrt(n))
list.result = list(theta.anti = theta.anti, SE.theta.anti = SE.theta.anti)
list.result
## $theta.anti
## [1] 0.5250573
## 
## $SE.theta.anti
## [1] 0.0004689239
#Problem 1 Part C
cor(g.function(u), g.function(1-u))
## [1] -0.9636159

Problem 2 Part a and Part b

n = 20000
F.inv <- function(x) {
    val = 1 - (1-exp(-1))*x
    return(-1*log(val))
}

unifs = runif(n, min=0, max=1)
x.vals = F.inv(unifs)

f.function <- function(x) {
    val = (exp(-x))/(1-exp(-1))
    return(val)
}

theta.estimate = mean(g.function(x.vals)/f.function(x.vals))
SE.estimate = sd(g.function(x.vals)/f.function(x.vals))/sqrt(n)

results.list = list(theta.hat = theta.estimate, SE = SE.estimate)
results.list
## $theta.hat
## [1] 0.525092
## 
## $SE
## [1] 0.0006821589

Problem 2 Part C

half.n = n/2
u.half = runif(half.n, min=0, max=1)
x.vals1 = F.inv(u.half)
x.vals2 = F.inv(1-u.half)
theta.antiX = (1/n)*(sum(g.function(x.vals1)/f.function(x.vals1) + g.function(x.vals2)/f.function(x.vals2)))

SE.theta.antiX = (1/sqrt(n))*sd((g.function(x.vals1)/f.function(x.vals1) + g.function(x.vals2)/f.function(x.vals2)))
list.result.anti = list(theta.antiX = theta.antiX, SE.theta.antiX = SE.theta.antiX)
list.result.anti
## $theta.antiX
## [1] 0.5244441
## 
## $SE.theta.antiX
## [1] 0.0003378534

Problem 3 Part A

library(MASS)
n = 20000 
mu = c(0,0,0)
sigma = matrix(c(1, 3/5, 1/3, 3/5, 1, 11/15, 1/3, 11/15, 1), 3, 3)
trivar = mvrnorm(n, mu, sigma)

indicator <- function(x) {
    x1 = x[1]
    x2 = x[2]
    x3 = x[3]
    if(x[1] <= 1 && x[2] <= 4 && x[3] <= 2) {
        return(1)
    } else {
        return(0)
    }
}

hit.miss = apply(trivar, 1, indicator)
p.region = mean(hit.miss) #approximates the probability of being in that
p.region
## [1] 0.8346

Problem 3 Part B

#standard error Part B
variance = (1/n)*p.region*(1-p.region)
se.estimate = sqrt(variance)
se.estimate
## [1] 0.002627193

Problem 3 Part C

z.star = qnorm(.975, 0, 1)
lower = p.region - z.star*(se.estimate)
upper = p.region + z.star*(se.estimate)
print(paste("(", lower, ", ", upper, ")"))
## [1] "( 0.82945079673468 ,  0.83974920326532 )"
true.val = 0.8279849897
if(true.val > lower && true.val < upper) {
  print("Yes the true value is captured by our 95% CI")
} else {
  print("No, our CI did not capture the true value")
}
## [1] "No, our CI did not capture the true value"

Problem 4 Part A

# Problem 4 PartA
n = 20000
u = runif(n, min=0, max=1)
h.function <- function(x) {
  return(exp(-0.5)/(1+x^2))
}

MU = pi/(4*exp(0.5))

gx = g.function(u)
hx = h.function(u)
c.star = -cov(gx, hx)/var(hx)

theta.CV = (1/n)*sum(g.function(u) + c.star*(hx-MU))
theta.CV
## [1] 0.5250983

Problem 4 Part B

We will first estimate the correlation between \(\frac{e^{-U}}{1+U^2}\) and \(\frac{e^{-0.5}}{1+U^2}\). Call these two expressions \(g(x)\) and \(h(x)\) respectively.

We first note that \(Cor(g(x), h(x)) = \frac{Cov(g(x), h(x))}{SD(h(x))SD(g(x))}\).

Then, by combing a derivation seen in Lecture with the formula for Correlation above we see that

\[ \begin{align} Var(\hat{\theta}_{CV}) &= Var(\hat{\theta}_{MC}) - \frac{1}{n}\frac{[Cov(g(x), h(x))]^2}{Var(h(x))}\\ &=Var(\hat{\theta}_{MC}) - \frac{1}{n}[Cor(g(x), h(x))]^2Var(g(x)) \end{align} \]

See the code below for an implementation of the above formula.

#Problem 4 Part B
corr = cor(g.function(u), h.function(u))

theta.MC = (1/n)*sum(g.function(u))
SE.theta.MC = sd(g.function(u))/sqrt(n)
Var.theta.MC = SE.theta.MC^2

reduc.term = (corr^2)*var(g.function(u))

Var.theta.CV = Var.theta.MC - (1/n)*reduc.term

SE.theta.CV = sqrt(Var.theta.CV)
SE.theta.CV
## [1] 0.0003951081

Problem 4 Part C

#Problem 4 Part C
SE.theta.CV.direct = (1/sqrt(n))*sd(g.function(u) + c.star*(hx-MU))
SE.theta.CV.direct
## [1] 0.0003951081

As we can see, using our derived formula that reduced the standard error of \(\hat{\theta}_{MC}\) (see Part B) gives use the same result as just calculating directly the standard error of our estimate.

Problem 4 Part D

Our estimator is the following:

\[ \hat{\theta}_{CV.X} = \frac{1}{n}\sum_{i=1}^n[\frac{g(X_i)}{f(X_i)} + c(\frac{h(X_i)}{f(X_i)} - \mu)] \] If \(X_i \sim f_{X}\), then we have the following:

\[ \begin{align} E[\hat{\theta}_{CV.X}] &= E[\frac{g(X)}{f(X)}] + c(E[\frac{h(X)}{f(X)}] - \mu)\\ &= \int_{0}^{1}\frac{g(x)}{f(x)}f(x)dx + c(\int_{0}^{1}\frac{h(x)}{f(x)}f(x)dx - \mu)\\ &= \int_{0}^{1}g(x)dx + c(\int_{0}^{1}h(x)dx - \mu)\\ &= \int_{0}^{1}g(x)dx + c(0)\\ & = \int_{0}^{1}g(x)dx \end{align} \]

Therefore, our estimator is unbiased.

We then have the following variance of our estimator:

\[ Var(\hat{\theta}_{CV.X}) = \frac{1}{n}[Var(\frac{g(X)}{f(X)}) + c^2Var(\frac{h(X)}{f(X)}) + 2cCov(\frac{g(X)}{f(X)}, \frac{h(X)}{f(X)})] \]

As a function of \(c\), this variance is minimized at \(c^{*} = \frac{Cov(\frac{g(X)}{f(X)}, \frac{h(X)}{f(X)})}{Var(\frac{h(X)}{f(X)})}\)

See code below for implementation. Note we first simulate \(n\) values from \(X\) by using the inverse function theorem, i.e., \(X = F^{-1}(U)\), where \(U \sim Unif(0,1)\).

#Problem 4 Part D

f.dens <- function(x) {
  return(2 - 2*x)
}

F.inv2 <- function(x) {
  val = 1 - sqrt(1-x)
  return(val)
}

n = 20000
u = runif(n, min=0, max=1)
X = F.inv2(u)

c.star2 = -cov(g.function(X)/f.dens(X), h.function(X)/f.dens(X))/var(h.function(X)/f.dens(X))
theta.CVX = (1/n)*sum(g.function(X)/f.dens(X) + c.star2*(h.function(X)/f.dens(X) - MU))
SE.theta.CVX = (1/sqrt(n))*sd(g.function(X)/f.dens(X) + c.star2*(h.function(X)/f.dens(X) - MU))
list.result = list(theta.CVX = theta.CVX, SE.theta.CVX = SE.theta.CVX)
list.result
## $theta.CVX
## [1] 0.5242105
## 
## $SE.theta.CVX
## [1] 0.0003399026