HW 2.Rmd

---
title: "HW 2"
author: "Emma Johnson"
date: "2/16/26"
output: pdf_document
---

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

# 1) Normal covered by Cauchy

To determine the optimal value for $c$ satisfying $\frac{f(y)}{g(y)} \leq c$, I will use the First Derivative Test on $\frac{f(y)}{g(y)}$ and evaluate the maximal $y$ in the ratio.

\begin{align}
\frac{f(y)}{c\hspace{2pt}g(y)} &\leq 1 \\
\implies \frac{f(y)}{g(y)} &\leq c
\end{align}

\begin{align}
\frac{f(y)}{g(y)} &= \frac{\frac{1}{\sqrt{2\pi}}e^{-y^2 / 2}}{\frac{1}{\pi (1+y^2)}} \\
&= \sqrt{\frac{\pi}{2}}\hspace{2pt}(1+y^2)\hspace{2pt}e^{-y^2 / 2}
\end{align}

Now, to find $\underset y {max} \{\hspace{2pt}(1+y^2)\hspace{2pt}e^{-y^2 / 2} \hspace{2pt}\}$,

\begin{align}
\frac{d}{dy} [\hspace{2pt}(1+y^2)\hspace{2pt}e^{-y^2 / 2}\hspace{2pt} ] &= 
(1+y^2)(-ye^{-y^2 / 2}) + (2y)\hspace{2pt}e^{-y^2 / 2} \\
&= -ye^{-y^2 / 2}(y^2-1) \\
&= 0 \\
\implies y &= 0, \pm 1
\end{align}

Plugging in these critical values for y,

\begin{align}
\frac{f(-1)}{g(-1)} &= \sqrt{\frac{\pi}{2}}\hspace{2pt}(1+(-1)^2)\hspace{2pt}e^{-(-1)^2 / 2} \\
&= \sqrt{2\pi}\hspace{2pt}e^{-1 / 2} \\
&= \frac{f(1)}{g(1)} \\
\\
\frac{f(0)}{g(0)} &= \sqrt{\frac{\pi}{2}}\hspace{2pt}(1+(0)^2)\hspace{2pt}e^{-(0)^2 / 2} \\
&= \sqrt{\frac{\pi}{2}} \\
&< {\sqrt{2\pi}\hspace{2pt}e^{-1 / 2}} \\
\\
\end{align}

The values $y=\pm 1$ yield the absolute maxima of the function, which can be reaffirmed by taking a look at the plot of $\frac{f(y)}{g(y)}$:

```{r plot 1}
curve(sqrt(pi/2)*(1+x^2)*exp(-x^2 /2), from = -5, to = 5, 
      lwd = 3, xlab = "y", ylab = "f(y) / g(y)")
```

So, the optimal $c$ is $\sqrt{2\pi}\hspace{2pt}e^{-1 / 2} \approx 1.520347$.

# 2) Normal covered by Double-Exponential

## a) Finding the optimal c

First, we will find the optimal value of $c$ that satisfies $\frac{f(y)}{g(y)} \leq c$ for $f(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2 /2}$, $-\infty < x < \infty$, and $g(y) = \frac{1}{2} e^{-|y|}$, $-\infty < y < \infty$.


\begin{align}
\frac{f(y)}{g(y)} &=  \frac{\frac{1}{\sqrt{2\pi}}e^{-y^2 /2}}{\frac{1}{2} e^{-|y|}} \\
&=  \sqrt{\frac{2}{\pi}} \hspace{2pt} e^{|y|-y^2 /2} \\
&\leq c \\
\end{align}


So we have $\frac{f(y)}{g(y)} \propto e^{|y|-y^2 /2}$. Note that it suffices to consider $y \geq 0$ here due to symmetry, so we will consider $e^{y-y^2 /2}$ in the following step. To optimize this expression to find the value $y$ that maximizes $\frac{f(y)}{g(y)}$, we will use the First Derivative Test on $e^{y-y^2 /2}$:


\begin{align}
\frac{d}{dy}(e^{y-y^2 /2}) &= (1-y)\hspace{2pt}e^{y-y^2 /2} \\
&= 0 \\
\\
\implies y &= 1
\end{align}


Values $y < 1$ have a derivative $>0$ and values $y > 1$ have a derivative $<0$, implying that $y=1$ is a local max. We can confirm this with a graph to see if it is the absolute max:

```{r plot 2}
curve(sqrt(2/pi)*exp(abs(x) - x^2 / 2), from = -5, to = 5, 
      lwd = 3, xlab = "y", ylab = "f(y) / g(y)")
```

It is confirmed that $y=1$ (and $y=-1$; symmetry) is an absolute max. Therefore, evaluating $y=1$ in $\frac{f(y)}{g(y)}$ will yield the optimal value of $c$.


\begin{align}
\frac{f(1)}{g(1)} &= \sqrt{\frac{2}{\pi}} \hspace{2pt} e^{|1|-1^2 /2} \\
&= \sqrt{\frac{2}{\pi}} \hspace{2pt} e^{1 /2} \\
&\approx 1.315489
\end{align}


Later on, in part c), I will use $c=\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{1 /2}$ to ensure that $\frac{f(y)}{c \hspace{2pt}g(y)} \leq 1$.

## b) General double-exponential

Now, we will consider a more general case of double-exponential with pdf $g(y) = \frac{\alpha}{2} e^{-\alpha|y|}$, $-\infty < y < \infty$. To find the optimal values of $\alpha$ and $c$, we will begin by maximizing $\frac{f(y)}{g(y)}$ with respect to $y$:


\begin{align}
\frac{f(y)}{g(y)} &=  \frac{\frac{1}{\sqrt{2\pi}}e^{-y^2 /2}}{\frac{\alpha}{2} e^{-\alpha|y|}} \\
&=  \frac{1}{\alpha}\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{\alpha|y|-y^2 /2} \\
&\leq c \\
\end{align}


Maximizing this quantity is easiest by considering $\frac{f(y)}{g(y)} \propto e^{\alpha|y|-y^2 /2}$. Once again, it suffices to consider $y \geq 0$ due to symmetry, so we will consider $e^{\alpha y-y^2 /2}$ in the following step. To optimize this expression to find the value $y$ that maximizes $\frac{f(y)}{g(y)}$, we will use the First Derivative Test on $e^{\alpha y-y^2 /2}$:


\begin{align}
\frac{d}{dy}(e^{\alpha y-y^2 /2}) &= (\alpha-y)\hspace{2pt}e^{\alpha y-y^2 /2} \\
&= 0 \\
\\
\implies y &= \alpha
\end{align}


Since the maximum $y$ of $\frac{f(y)}{g(y)}$ is a function of $\alpha$, we can use this to determine the optimal value first by plugging in $y = \alpha$ into the ratio, then by taking the Maximum Likelihood approach to solve for $\alpha$:


\begin{align}
\frac{f(\alpha)}{g(\alpha)} &=  \frac{1}{\alpha}\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{\alpha|\alpha|-\alpha^2 /2} \\
&= \frac{1}{\alpha}\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{\alpha^2-\alpha^2 /2}\\
&= \frac{1}{\alpha}\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{\alpha^2 /2}\\
\\
\log{\frac{f(\alpha)}{g(\alpha)}} &= -\log{\alpha} + \frac{1}{2}\log{\frac{2}{\pi}} + \frac{1}{2}{\alpha ^2}\\
\\
\frac{d}{d\alpha}\log{\frac{f(\alpha)}{g(\alpha)}} &= -\frac{1}{\alpha}+\alpha \\
&= 0\\
\\
\implies& \alpha = \frac{1}{\alpha}\\
\\
\implies& \alpha^2 = 1\\
\\
\implies& \alpha = \pm 1
\end{align}


In practice, $\alpha >0$, so we will let $\alpha = 1$ be the optimal quantity. Note that $\alpha = 1$ makes the pdf $g(y) = \frac{\alpha}{2} e^{-\alpha|y|} = \frac{1}{2} e^{-|y|}$, exactly equivalent to the pdf used in part a). So, we already know the optimal value for $y$, and thus, $c$, from the previous work. Therefore, $\alpha = 1$ and $c = \sqrt{\frac{2}{\pi}} \hspace{2pt} e^{1 /2} \approx 1.315489$ are optimal. 

## c) Generating a sample

Using the optimal values $\alpha = 1$ and $c=\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{1 /2}$, here is the acceptance-rejection method to generate a sample of size $n= 1000$ from $f(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2 /2}$, $-\infty < x < \infty$ using the covering density $g(y) = \frac{1}{2} e^{-|y|}$. 

Note that we have $\frac{f(y)}{c\hspace{2pt}g(y)} = \frac{1}{c}\sqrt{\frac{2}{\pi}} \hspace{2pt} e^{|y|-y^2 /2}$. Further note that in the code, I wrote `y <- sample(c(-1, 1), 1) * rexp(1)`, which translates to $y \in \{ -\infty, \infty \}$ with the `sample()` part, making use of `rexp()`, which would limit $y>0$, otherwise.

```{r part 2 sample}
set.seed(5)

n <- 1000
k <- 0 # counter for accepted
j <- 0 # iterations
x1 <- numeric(n)
c <- sqrt(2/pi) * exp(1/2) # approx 1.315489

while (k < n) {
  u <- runif(1)
  j <- j + 1
  y <- sample(c(-1, 1), 1) * rexp(1)  # random variate from g 
  if (u < 1/c * sqrt(2/pi) * exp(abs(y) - y^2 / 2)) {
    # accepting the current value
    k <- k + 1
    x1[k] <- y
  }
}
j

hist(x1, prob = TRUE, xlim = c(-4, 4)); curve(1/sqrt(2*pi) * exp(- x^2 / 2), from = -4, 
                                                          to = 4, add = TRUE, 
                                                          col = "red", lwd = 2)
```

## ALTERNATIVELY

Consider defining `y` based on the Inverse Transform method. Then, let $Y \sim g(y) = \frac{1}{2} e^{-|y|}$ for $-\infty < y < \infty$ and $U \sim \text{Unif}(0,1)$. First, I will derive the CDF of $Y$:

Case 1: $y < 0$

\begin{align}
G_Y(y) &= \int_{-\infty}^{y} \frac{1}{2} e^{-|t|} \ dt \\
&= \frac{1}{2} \int_{-\infty}^{y} e^{t} \ dt \\
&= \frac{1}{2} (e^{t}) /_{-\infty}^{\ y} \\
&= \frac{1}{2} (e^{y} - 0) \\
&= \frac{1}{2} e^{y} ; \ \ y<0\\
\end{align}

Case 2: $y \geq 0$

\begin{align}
G_Y(y) &= \int_{-\infty}^{y} \frac{1}{2} e^{-|t|} \ dt \\
&= \int_{-\infty}^{0} \frac{1}{2} e^{t} \ dt \ +  \int_{0}^{y} \frac{1}{2} e^{-t} \ dt\\
&= \frac{1}{2} [ \ (e^{t}) /_{-\infty}^{\ 0} \ + \ (-e^{-t}) /_{0}^{\ y} \ ]\\
&= \frac{1}{2} [ \ (1 - 0) \ + \ (1 -e^{-y}) \ ] \\
&= 1- \frac{1}{2} e^{-y} ; \ \ y \geq 0\\
\end{align}

So we have the CDF 

$$
G_Y(y) =
\begin{cases}
\frac{1}{2} e^{y}, & y<0 \\\\
1 - \frac{1}{2} e^{-y}, & y \geq 0
\end{cases}
$$

Now let $U = G_Y(y)$. Then, $u< \frac{1}{2} \ \ {\small(= \frac{1}{2} e^{0})}$ when $y<0$ and $u \geq \frac{1}{2} \ \ {\small(=1 - \frac{1}{2} e^{0})}$ when $y \geq 0$.

#### FIX Y or y?
The inverse $Y =G_Y^{-1}(u)$ can be found by setting $u=\frac{1}{2} e^{y}$ and $u = 1 - \frac{1}{2} e^{-y}$ for each respective case, summarized below:

$$
y =
\begin{cases}
\log(2u), & 0<u<\frac{1}{2} \\\\
-\log(2(1-u)), & \frac{1}{2} \leq u < 1
\end{cases}
$$

Using this transform method, we can generate random variables $y$ from $g(y) = \frac{1}{2} e^{-|y|}$ by utilizing $u$ from UNIF(0,1).

```{r part 2 pt 2}
set.seed(5)

n <- 1000
k <- 0 # counter for accepted
j <- 0 # iterations
x2 <- numeric(n)
c <- sqrt(2/pi) * exp(1/2) # approx 1.315489

while (k < n) {
  j <- j + 1
  
  u1 <- runif(1)
  if (u1 < 0.5) {
    y <- log(2*u1)
  } else {
    y <- -log(2*(1-u1))
  }
  
  u2 <- runif(1)
  if (u2 < 1/c * sqrt(2/pi) * exp(abs(y) - y^2 / 2)) {
    # accepting the current value
    k <- k + 1
    x2[k] <- y
  }
}
j

hist(x2, prob = TRUE, xlim = c(-4, 4)); curve(1/sqrt(2*pi) * exp(- x^2 / 2), from = -4, 
                                                          to = 4, add = TRUE, 
                                                          col = "red", lwd = 2)
```


By the looks of the histogram, the sample appears to be pretty normally ditributed. We can confirm this with a quick normality test:

```{r normality test}
shapiro.test(x1)
shapiro.test(x2)
```
 
The null hypothesis of the Shapiro-Wilk normality test is $H_0:$ the data come from a normal distribution. Since $p$ is large (for both methods), we can accept the null hypothesis and conclude that our samples likely follow a normal distribution.

# 3) Comparing Performances

Note that the acceptance-rejection function I've written omits the return of `j` or `k` to try and optimize performance. 

```{r acc rej alg}
set.seed(5)

acc_rej_beta <- function(n, a, b) {
  
  x1 <- seq(from = 0, to = 1, by = 1 / n) ## by 0.001 or 1/n?
  
  beta_coef <- gamma(a+b)/gamma(a)/gamma(b)
  
  c <- ceiling(max(beta_coef * x1^(a-1) * (1-x1)^(b-1)))
  
  k <- 0
  x <- numeric(n)

  while (k < n) {
    
    u <- runif(1)
    y <- runif(1)
    
    if (u < beta_coef * y^(a - 1) * (1-y)^(b-1) / c) {
      k <- k+1
      x[k] <- y
    }
  }
  return(x)
}

acc.rej1 <- system.time(acc_rej_beta(n = 1e+6,  a = 2, b = 3))
acc.rej2 <- system.time(acc_rej_beta(n = 1e+6,  a = 5, b = 5))
```

```{r transformation gamma beta}
set.seed(5)

mybeta <- function(n, a, b) {
  u <- rgamma(n, a, 1)
  v <- rgamma(n, b, 1)
  x <- u / (u + v)
  return(x)
}

my.beta1 <- system.time(mybeta(n = 1e+6,  a = 2, b = 3))
my.beta2 <- system.time(mybeta(n = 1e+6,  a = 5, b = 5))
```

```{r rbeta}
set.seed(5)

rbeta1 <- system.time(rbeta(n = 1e+6, shape1 = 2, shape2 = 3))
rbeta2 <- system.time(rbeta(n = 1e+6, shape1 = 5, shape2 = 5))
```

```{r comparisons, echo = FALSE}
as.data.frame(rbind(acc.rej1, acc.rej2, my.beta1, my.beta2, rbeta1, rbeta2))[,-c(4:5)]
```

Note that the indices with a suffix `1` refer to $n = 10^6$, $a = 2$, and $b=3$, whereas the suffix `2` refers to $n = 10^6$ and $a = b = 5$.

It appears that the acceptance-rejection method takes much longer than the transformation and `rbeta` methods combined. The transformation method and `rbeta` methods are comparable, though running the transformation method with $a = b = 5$ is by far the quickest. It is interesting to see this one so much more efficient than the `rbeta` method with the same parameters.