HW 4.Rmd

---
title: "HW 4"
author: "Emma Johnson"
date: "4/1/26"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# 1) 

Suppose we would like to estimate

\begin{align*}
\theta = \int_0^1 \frac{e^{-x}}{1 + x^2} \, dx
\end{align*}


First, consider the graph of the integrand $g(x) = \frac{e^{-x}}{1 + x^2}$ accompanied by the envelope density $f_3 (x) = \frac{e^{-x}}{1 - e^{-1}}$ for $0<x<1$:

```{r graph1}
x <- seq(0, 1, .01)
f3 <- function(x) exp(-x) / (1 - exp(-1))
g <- function(x) exp(-x) / (1 + x^2)

par(mfrow = c(1, 2))

#figure (a)
plot(x, g(x), type = "l", main = "(a)", ylab = "", ylim = c(0,2), lwd = 3)
lines(x, f3(x), lwd = 2, col = "red")
legend("topright", legend = c("g(x)", "f3(x)"), lwd = c(3,2), col = c("black", "red"),
       inset = 0.02)

#figure (b)
plot(x, g(x)/f3(x), type = "l", main = "(b)", ylab = "", ylim = c(0,2), lty = 5, lwd = 3)
legend("topright", legend = c("g(x)/f3(x)"), lty = 5, lwd = 2, inset = 0.02)
```

The importance function $f_3 (x)$ has a similar shape to $g(x)$, as seen in **(a)**. Further, **(b)** shows how the ratio $\frac{g(x)}{f_3 (x)}$ is near constant over $0 < x < 1$. Therefore, $f_3 (x)$ appears to be a decent choice for importance sampling here.


To find an importance sampling estimator $\hat \theta_1$, we can manipulate how we write $\theta$:


\begin{align*}
\theta &= \int_0^1 \frac{e^{-x}}{1 + x^2} \, dx \\
&= \int_0^1 \frac{\frac{e^{-x}}{1 + x^2}}{\frac{e^{-x}}{1 - e^{-1}}} \, (\frac{e^{-x}}{1 - e^{-1}})\, dx \\
&= \int_0^1 \frac{g(x)}{f_3 (x)} \, (f_3 (x))\, dx \\
&= E_{f_3}[\frac{g(x)}{f_3 (x)}] \\
\\
\implies \hat \theta_1 &= \frac{1}{m} \sum_{i = 1}^m \frac{g(x_i)}{f_3 (x_i)} \\
&= \frac{1}{m} \sum_{i = 1}^m \frac{\frac{e^{-x_i}}{1 + x^2}}{\frac{e^{-x_i}}{1 - e^{-1}}} \\
&= \frac{1}{m} \sum_{i = 1}^m \frac{1 - e^{-1}}{1 + x_i^2}
\end{align*}


where $X_i \sim f_3 (x)$.

## Employing Importance Sampling in R

As shown in class, we can generate $m = 1000$ random variables for $X \sim f_3 (x)$ using the inverse transform method where $X = F_X^{-1}(U)$ for $U \sim \text{UNIF}(0,1)$. First, we must find the cdf $F_X(x)$ based on $f_3 (x)$:


\begin{align*}
F_X(x) &= \int_0^x f_3 (t) \, dt \\
&= \int_0^x \frac{e^{-t}}{1 - e^{-1}} \, dt \\
&= -\frac{e^{-t}}{1 - e^{-1}} /_0^{\,x} \\
&= \frac{1 - e^{-x}}{1 - e^{-1}} ,\,\,\,\,\,0<x<1
\end{align*}

For our sample, let $u = F_X(x)$. Then,

\begin{align*}
u &= \frac{1 - e^{-x}}{1 - e^{-1}} \\
\\
\implies x &= -\log (1 - u(1 - e^{-x}))
\end{align*}



Now, we can calculate $\hat \theta_1$ in R as follows.

```{r f3}
set.seed(5)
m <- 1000
u <- runif(m)
x <- -log(1-u*(1-exp(-1)))
theta.hat1 <- mean(g(x) / f3(x))
theta.hat1

# exact value
integrate(g, 0, 1)
```

The estimator $\hat \theta_1$ comes pretty close to the actual value of the integral, only about 0.001 off. Let's see how this compares with an antithetic variate implementation.

## Antithetic Approach

Using the antithetic variate approach, we will generate $m/2 = 500$ random variables $U \sim \text{UNIF}(0,1)$ where

Let this new estimator be $\hat \theta_2$, defined as:


\begin{align*}
\hat \theta_2 &= \frac{1}{m} \sum_{j = 1}^{m/2}  
\frac{g(x_j)}{f_3 (x_j)} + \frac{g(x'_j)}{f_3(x'_j)} \\
&= \frac{1}{m} \sum_{j = 1}^{m/2}  \frac{1 - e^{-1}}{1 + x_j^2} + \frac{1 - e^{-1}}{1 + {x'_j}^2}
\end{align*} 


where $x_j = - \log (1-u_j(1-e^{-1}))$ for $j = 1, 2, \ldots , 500$ and $x'_j = - \log (1-(1-u_j)(1-e^{-1}))$ for $j = 501, 502, \ldots , 1000$. 


Now, implementing the antithetic approach in R:

```{r anti time}
set.seed(5)
m <- 1000 # sample size
f3 <- function(x) exp(-x) / (1 - exp(-1))
g <- function(x) exp(-x) / (1 + x^2)

# point estimate for antithetic estimator:
u <- runif(m/2)
x <- c(-log(1-u*(1-exp(-1))), -log(1-(1-u)*(1-exp(-1))))
theta.hat2 <- mean(g(x)/f3(x))
theta.hat2

# to find the variances:
n <- 1000 # number of replicates
MC1 <- MC2 <- numeric(n)
for (i in 1:n){
  # theta.hat1
  y <- -log(1-runif(m)*(1-exp(-1)))
  MC1[i] <- mean(g(y)/f3(y))
  
  # theta.hat2
  u <- runif(m/2)
  x <- c(-log(1-u*(1-exp(-1))), -log(1-(1-u)*(1-exp(-1))))
  MC2[i] <- mean(g(x)/f3(x))
}

var(MC1) # theta.hat1
var(MC2) # theta.hat2

# Reduction in variance:
(var(MC1) - var(MC2))/var(MC1)
```

Using a antithetic variables, we obtain $\hat \theta_2 = 0.5236352$, which deviates from the actual value of $\theta$ by around 0.001, about the same as $\hat\theta_1$. However, one approach appears to out-favor the other due to a small variance. The antithetic approach yields a much smaller variance than implementing importance sampling without it. That is, the simulated reduction in variance for using $\hat \theta_2$ over $\hat \theta_1$ is about 87.6\%. 

# 2)

Given some function $g(x) = x^2$ the standard normal pdf $f(x) = \frac{1}{\sqrt{2\pi}} \, e^{-x^2 /2}$, consider estimating


\begin{align*}
\theta &= \int_1^\infty \frac{x^2}{\sqrt{2\pi}} \, e^{-x^2 /2} \, dx \\
&= \int_1^\infty g(x) \, f(x) \, dx
\end{align*}


## Envelope Densities

The envelope densities I will try out are a truncated Cauchy and a shifted Exponential distribution. Specifically, consider $\varphi_1 (x) = \frac{4}{\pi (1 + x^2)}$ and $\varphi_2 (x) = e^{-(x - 1)}$ for $x > 1$:


\begin{align*}
\int_1^\infty \varphi_1 (x) \, dx &= \int_1^\infty \frac{4}{\pi (1 + x^2)} \, dx \\
&= \frac{4}{\pi} (\arctan x /_1^\infty) \\
&= \frac{4}{\pi} (\frac{\pi}{2} - \frac{\pi}{4}) \\
&= 1 \\
\\
\int_1^\infty \varphi_2 (x) \, dx &= \int_1^\infty e^{-(x - 1)} \, dx \\
&= - e^{-(x - 1)}/_1^\infty \\
&=  0 + 1\\
&= 1
\end{align*}


Both $\varphi_1 (x)$ and $\varphi_2 (x)$ are proper pdfs over $1 < x < \infty$.

Further, consider graphs of the functions:

```{r graph2}
x <- seq(1, 10, .01)
f <- dnorm(x)
g <- x^2
fg <- f*g
phi1 <- 4/(pi*(1+x^2))
phi2 <- exp(-(x-1))

plot(x, fg, type = "l", main = "", ylab = "", xlim = c(1, 10), ylim = c(0,1), lwd = 3)
lines(x, phi1, lwd = 2, col = "red")
lines(x, phi2, lwd = 2, col = "blue")
legend("topright", legend = c("g(x)f(x)", "phi1(x)", "phi2(x)"), lwd = c(3,2,2), 
       col = c("black", "red", "blue"), inset = 0.02)
```

Notice how both $\varphi_1 (x)$ and $\varphi_2 (x)$ behave similarly overall to $g(x) \, f(x)$, peaking to the left and tailing off to the right. Further, it appears that both envelopes have heavier tails than $g(x) \, f(x)$.

To verify that $\varphi_1 (x)$ and $\varphi_2 (x)$ have heavier tails than $g(x) \, f(x) = \frac{x^2}{\sqrt{2\pi}} \, e^{-x^2 /2}$, consider:


\begin{align*}
\lim_{x\to \infty} \frac{g(x) \, f(x)}{\varphi_1 (x)} &= 
\lim_{x\to \infty} \frac{\frac{x^2}{\sqrt{2\pi}} \, e^{-x^2 /2}}{\frac{4}{\pi (1 + x^2)}} \\
&= \lim_{x\to \infty} \frac{\pi}{4 \sqrt{2\pi}} \ x^2 \, (1 + x^2)\, e^{-x^2 /2} \\
&= 0 \\
\\
\lim_{x\to \infty} \frac{g(x) \, f(x)}{\varphi_2 (x)} &= 
\lim_{x\to \infty} \frac{\frac{x^2}{\sqrt{2\pi}} \, e^{-x^2 /2}}{e^{-(x - 1)}} \\
&= \lim_{x\to \infty} \frac{\pi}{4 \sqrt{2\pi}} \ x^2 \, e^{x - 1}\, e^{-x^2 /2} \\
&= 0 
\end{align*}


where, in both cases, $e^{-x^2 /2}$ is the dominating term. This confirms that both $\varphi_1 (x)$ and $\varphi_2 (x)$ converge slower than $g(x) \, f(x)$, thus, having heavier tails.

As for the actual importance sampling procedure, we can generate from both of these envelope densities using the inverse transform method. It was already shown how to integrate each pdf, so let's take it a step further to solve for our transformations using $U = F^{-1}_X (x)$:


\begin{align*}
\int_1^x \varphi_1 (t) \, dt &= \frac{4}{\pi} \arctan t /_1^x \\
&= \frac{4}{\pi} \arctan x -1  = u \\
\implies x &= \tan(\frac{\pi}{4}(u + 1)),\,\,\,\,\,0<u<1\\
\\
\int_1^x \varphi_2 (t) \, dt &= - e^{-(t - 1)}/_1^x \\
&=  1 - e^{-(x - 1)} = u \\
\implies x &= 1 - \log (1- u),\,\,\,\,\,0<u<1
\end{align*}

Using these transformations, let us find the estimators $\hat\theta_1 = \frac{1}{m} \sum_{i = 1}^m \frac{g(x_i) \, f(x_i)}{\varphi_1(x_i)}$ for the truncated Cauchy and for $\hat\theta_2 = \frac{1}{m} \sum_{i = 1}^m \frac{g(x_i) \, f(x_i)}{\varphi_2(x_i)}$ for the shifted exponential. Note that we will use sample size $m = 1000$.


```{r imp samp}
set.seed(5)
m <- 1000
fg <- function(x) dnorm(x)*x^2
phi1 <- function(x) 4/(pi*(1+x^2))
phi2 <- function(x) exp(-(x-1))

# phi1 - truncated Cauchy
u <- runif(m)
x <- tan((pi/4)*(u+1))
theta.hat1 <- mean(fg(x) / phi1(x))
theta.hat1

# phi2 - shifted exponential
u <- runif(m)
x <- 1-log(1-u)
theta.hat2 <- mean(fg(x) / phi2(x))
theta.hat2

# exact value:
integrate(fg, 1, Inf)

# variances:
n <- 1000
MC1 <- MC2 <- numeric(n)
for (i in 1:n){
  # phi1
  x1 <- tan((pi/4)*(runif(m)+1))
  MC1[i] <- mean(fg(x1) / phi1(x1))
  
  # phi2
  x2 <- 1-log(1-runif(m))
  MC2[i] <- mean(fg(x2) / phi2(x2))
}
var(MC1)
var(MC2)

# reduction in variance:
(var(MC1) - var(MC2))/var(MC1)
```

Based on the simulation results, I'd recommend using $\varphi_2$. While both estimators $\hat \theta_1$ and $\hat \theta_2$ come close to the actual value $0.400626$ produced by `integrate(.)`, I'd recommend using $\varphi_2$ since it yields a smaller simulated variance than choosing $\varphi_1$ as the envelope density. In fact, choosing $\varphi_2$ results in an approximate 76.8\% reduction in variance when compared to $\varphi_1$.

# 3)


\begin{align*}
\theta = \int_1^\infty \frac{x^2}{\sqrt{2\pi}} \, e^{-x^2 /2} \, dx
\end{align*}


Now, we will use $\varphi_2 (x) = e^{-(x - 1)}$ as our envelope density to proceed with the antithetic variable approach with sample size $m = 1000$. Let this new estimator be $\hat \theta_3$, defined as:


\begin{align*}
\hat \theta_3 = \frac{1}{m} \sum_{j = 1}^{m/2} 
\frac{g(x_j) \, f(x_j)}{\varphi_2 (x_j)} + \frac{g(x'_j) \, f(x'_j)}{\varphi_2 (x'_j)}
\end{align*} 


where $x_j = 1 - \log (1- u_j)$ for $j = 1, 2, \ldots , 500$ and $x'_j = 1 - \log (u_j)$ for $j = 501, 502, \ldots , 1000$. 


```{r anti 3}
set.seed(5)
n <- 1000
m <- 1000
MC2 <- MC3 <- numeric(n)
# for a single estimate
u <- runif(m/2)
x <- c(1-log(1-u), 1-log(u))
theta.hat3 <- mean(fg(x) / phi2(x))
theta.hat3

# for multiple replicates
for (i in 1:n){
y <- 1-log(1-runif(m))
MC2[i] <- mean(fg(y) / phi2(y)) # original estimate
u <- runif(m/2)
x <- c(1-log(1-u), 1-log(u))
MC3[i] <- mean(fg(x) / phi2(x)) # antithetic
}
var(MC2)
var(MC3)

# reduction in variance:
(var(MC2) - var(MC3))/var(MC2)
```

The antithetic approach does not lead to a reduction in variance. Instead, more variance is accumulated, suggesting against using the antithetic approach here.

This is likely due to the fact that $g(x) \, f(x)$ is not monotone over $1 < x < \infty$. Specifically, we can check if $Cov(\frac{g(X) \, f(X)}{\varphi_2(X)}, \frac{g(X') \, f(X')}{\varphi_2(X')}) > 0$, given that  $Cov(U, 1-U) < 0$ (recall that $X$ is generated using $U$ and $X'$ is generated using $1-U$ by the inverse transform method).

We can test this hypothesis with one replicate, as follows:

```{r neg corrs}
set.seed(5)
m <- 1000

u <- runif(m/2)
# verifying negative correlation b/w u & 1-u
cor(u, 1-u) < 0 

x <- c(1-log(1-u), 1-log(u))
# verifying negative correlation b/w x & x'
cor(x[1:(m/2)], x[(m/2 + 1):m]) < 0

fxn.vals <- fg(x) / phi2(x)
# checking  for negative correlation b/w 
# g(x)f(x)/phi1(x) and g(x')f(x')/phi1(x')
cor(fxn.vals[1:(m/2)], fxn.vals[(m/2 + 1):m]) < 0
```

As expected, $Cov(\frac{g(X) \, f(X)}{\varphi_2(X)}, \frac{g(X') \, f(X')}{\varphi_2(X')})$ is positive, which could explain the extra variance in $\hat \theta_3$. If $g(x) \, f(x)$ was monotone over $x>1$, then we'd expect the covariance using antithetic variables to be negative. However, this is not the case, so we cannot rely on the antithetic approach here.

```{r neg corrs}
set.seed(5)
m <- 1000


u <- runif(m/2)
# verifying negative correlation b/w u & 1-u
cor(u, 1-u) < 0 

x <- c(1-log(1-u), 1-log(u))
# verifying negative correlation b/w x & x'
cor(x[1:(m/2)], x[(m/2 + 1):m]) < 0

fxn.vals <- fg(x) / phi2(x)
# checking  for negative correlation b/w 
# g(x)f(x)/phi1(x) and g(x')f(x')/phi1(x')
cor(fxn.vals[1:(m/2)], fxn.vals[(m/2 + 1):m]) < 0
```