1 (1)

1.1 (a) Prove dist of param given data

1.2 (b) Express data given param in terms of data avr.

1.3 (c)

x <- seq(0, 8, length = 200)

prior <- dgamma(x, shape = 6, rate = 2)
post1 <- dgamma(x, shape = 26, rate = 7)

prior_mean <- 6 / 2
sample_mean <- 20 / 5
posterior_mean <- 26 / 7

ymax <- max(prior, post1)

plot(x, prior, type = "l", lwd = 2, col = "blue",
     ylim = c(0, ymax * 1.1),
     ylab = "Density", xlab = expression(lambda),
     main = "Prior and Posterior: Case 1")

lines(x, post1, col = "red", lwd = 2)

# Vertical lines for the means
abline(v = prior_mean, col = "blue", lwd = 2, lty = 2)
abline(v = sample_mean, col = "darkgreen", lwd = 2, lty = 2)
abline(v = posterior_mean, col = "purple", lwd = 2, lty = 3)

legend("topright",
       legend = c("Prior density",
                  "Posterior density",
                  "Prior mean = 3",
                  "Sample mean = 4",
                  "Posterior mean = 3.71"),
       col = c("blue", "red", "blue", "darkgreen", "purple"),
       lwd = 2,
       lty = c(1, 1, 2, 2, 3),
       cex = 0.5)

x <- seq(0, 8, length = 200)

prior <- dgamma(x, shape = 6, rate = 2)
post2 <- dgamma(x, shape = 86, rate = 22)

prior_mean <- 6 / 2
sample_mean <- 80 / 20
posterior_mean <- 86 / 22

ymax <- max(prior, post2)

plot(x, prior, type = "l", lwd = 2, col = "blue",
     ylim = c(0, ymax * 1.1),
     ylab = "Density", xlab = expression(lambda),
     main = "Prior and Posterior: Case 2")

lines(x, post2, col = "red", lwd = 2)

# Vertical lines for the means
abline(v = prior_mean, col = "blue", lwd = 2, lty = 2)
abline(v = sample_mean, col = "darkgreen", lwd = 2, lty = 2)
abline(v = posterior_mean, col = "purple", lwd = 2, lty = 3)

legend("topright",
       legend = c("Prior density",
                  "Posterior density",
                  "Prior mean = 3",
                  "Sample mean = 4",
                  "Posterior mean = 3.91"),
       col = c("blue", "red", "blue", "darkgreen", "purple"),
       lwd = 2,
       lty = c(1, 1, 2, 2, 3),
       cex = .5)

notice with more samples we become more certain on out belief of the parameter

2 (2)

Intuition Reference : Importance Sampling - Why Where You Sample Beats How Often

2.1 (a)

2.1.1 1. Shifted Exponential Distribution (Rate = 1)

Choose the importance density

\[ g_1(x)=e^{-(x-1)}, \qquad x>1. \]

To generate samples,

\[ Y = 1 + E, \qquad E \sim \text{Exp}(1). \]

The corresponding importance weight is

\[ w_1(x) = \frac{\phi(x)}{g_1(x)} = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}+x-1\right). \]

2.1.2 2. Shifted exp Distribution (Rate = 2)

Choose the importance density

\[ g_2(x)=2e^{-2(x-1)}, \qquad x>1. \]

To generate samples,

\[ Y = 1 + E, \qquad E \sim \text{Exp}(2). \]

The corresponding importance weight is

\[ w_2(x) = \frac{\phi(x)}{g_2(x)} = \frac{1}{2\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}+2x-2\right). \]

2.1.3 3. Truncated Normal Distribution

Suppose

\[ Z \sim N(1,1). \]

Since

\[ P(Z>1)=\frac12, \]

the conditional density given \(Z>1\) is

\[ g_3(x) = 2\phi(x-1), \qquad x>1. \]

An equivalent sampling method is

\[ Y = 1 + |Z|, \qquad Z \sim N(0,1). \]

The corresponding importance weight is

\[ w_3(x) = \frac{\phi(x)}{2\phi(x-1)} = \frac12 \exp\left(-x+\frac12\right). \]

2.1.4 Variance comparison: intuition and theory

For an importance density \(g\), the estimator of

\[ \theta=P(X>1)=\int_1^\infty \phi(x)\,dx \]

is

\[ \widehat{\theta}_g = \frac{1}{n}\sum_{i=1}^n \frac{\phi(Y_i)}{g(Y_i)}, \qquad Y_i\sim g. \]

Define the importance weight

\[ W_g(Y)=\frac{\phi(Y)}{g(Y)}. \]

Because the \(Y_i\)’s are independent,

\[ \operatorname{Var}(\widehat{\theta}_g) = \frac{1}{n}\operatorname{Var}_g(W_g(Y)). \]

Also,

\[ E_g[W_g(Y)] = \int_1^\infty \frac{\phi(x)}{g(x)}g(x)\,dx = \int_1^\infty \phi(x)\,dx = \theta. \]

Therefore,

\[ \boxed{ \operatorname{Var}(\widehat{\theta}_g) = \frac{1}{n} \left[ \int_1^\infty \frac{\phi(x)^2}{g(x)}\,dx -\theta^2 \right] } \]

where

\[ \theta=1-\Phi(1). \]

The variance depends on how much the importance weights

\[ \frac{\phi(x)}{g(x)} \]

change across sampled values. An importance density that produces nearly constant weights will have a smaller variance.

2.1.5 Variance using \(g_1\)

For

\[ g_1(x)=e^{-(x-1)},\qquad x>1, \]

the second moment of the importance weight is

\[ E_{g_1}[W_1^2] = \int_1^\infty \frac{\phi(x)^2}{g_1(x)}\,dx. \]

Since

\[ \phi(x)^2 = \frac{1}{2\pi}e^{-x^2}, \]

we obtain

\[ E_{g_1}[W_1^2] = \frac{1}{2\pi} \int_1^\infty e^{-x^2+x-1}\,dx. \]

Thus,

\[ \operatorname{Var}(\widehat{\theta}_{g_1}) = \frac{1}{n} \left[ \int_1^\infty \frac{\phi(x)^2}{g_1(x)}\,dx -\theta^2 \right]. \]

2.1.6 Variance using \(g_2\)

For

\[ g_2(x)=2e^{-2(x-1)},\qquad x>1, \]

we have

\[ E_{g_2}[W_2^2] = \int_1^\infty \frac{\phi(x)^2}{g_2(x)}\,dx. \]

Therefore,

\[ E_{g_2}[W_2^2] = \frac{1}{4\pi} \int_1^\infty e^{-x^2+2x-2}\,dx, \]

and

\[ \operatorname{Var}(\widehat{\theta}_{g_2}) = \frac{1}{n} \left[ \int_1^\infty \frac{\phi(x)^2}{g_2(x)}\,dx -\theta^2 \right]. \]

2.1.7 Variance using \(g_3\)

For

\[ g_3(x)=2\phi(x-1),\qquad x>1, \]

the second moment is

\[ E_{g_3}[W_3^2] = \int_1^\infty \frac{\phi(x)^2}{2\phi(x-1)}\,dx. \]

The corresponding variance is

\[ \operatorname{Var}(\widehat{\theta}_{g_3}) = \frac{1}{n} \left[ \int_1^\infty \frac{\phi(x)^2}{g_3(x)}\,dx -\theta^2 \right]. \]

2.2 R implementation

phi <- function(x) dnorm(x)

g1 <- function(x) {
  ifelse(x > 1, exp(-(x - 1)), 0)
}

g2 <- function(x) {
  ifelse(x > 1, 2 * exp(-2 * (x - 1)), 0)
}

g3 <- function(x) {
  ifelse(x > 1, 2 * dnorm(x - 1), 0)
}

x <- seq(1.001, 5, length.out = 1000)

par(mfrow = c(1, 2))

# Compare the target integrand and importance densities
plot(x, phi(x), type = "l", lwd = 3, col = "black",
     xlab = "x", ylab = "Density",
     main = "Target and importance densities")
lines(x, g1(x), col = "red", lwd = 2)
lines(x, g2(x), col = "blue", lwd = 2)
lines(x, g3(x), col = "darkgreen", lwd = 2)

legend("topright",
       legend = c("phi(x)", "g1: shifted Exp(1)",
                  "g2: shifted Exp(2)", "g3: truncated N(1,1)"),
       col = c("black", "red", "blue", "darkgreen"),
       lty = 1, lwd = c(3, 2, 2, 2), cex = 0.3)

# Compare importance weights
plot(x, phi(x) / g1(x), type = "l", lwd = 2, col = "red",
     xlab = "x", ylab = expression(phi(x) / g(x)),
     main = "Importance weights")
lines(x, phi(x) / g2(x), col = "blue", lwd = 2)
lines(x, phi(x) / g3(x), col = "darkgreen", lwd = 2)

abline(h = 1 - pnorm(1), lty = 2, col = "black")

legend("topright",
       legend = c("g1", "g2", "g3", "Ideal constant weight"),
       col = c("red", "blue", "darkgreen", "black"),
       lty = c(1, 1, 1, 2), lwd = 2, cex = 0.3)

2.3 (b)

IS_estimate <- function(n, proposal = c("g1", "g2", "g3")) {

  proposal <- match.arg(proposal)

  if (proposal == "g1") {
    # Y = 1 + Exp(1)
    y <- 1 + rexp(n, rate = 1)
    w <- phi(y) / g1(y)

  } else if (proposal == "g2") {
    # Y = 1 + Exp(2)
    y <- 1 + rexp(n, rate = 2)
    w <- phi(y) / g2(y)

  } else if (proposal == "g3") {
    # Y = 1 + |Z|, Z ~ N(0,1)
    y <- 1 + abs(rnorm(n))
    w <- phi(y) / g3(y)
  }

  mean(w)
}
set.seed(123)

IS_estimate(10000, "g1")
## [1] 0.1583584
IS_estimate(10000, "g2")
## [1] 0.1587922
IS_estimate(10000, "g3")
## [1] 0.1578342
print("actual good est : ")
## [1] "actual good est : "
1 - pnorm(1)
## [1] 0.1586553

2.4 (c)

bootstrap_IS <- function(n = 10000, B = 1000, seed = 123) {

  set.seed(seed)

  # Generate importance samples
  y1 <- 1 + rexp(n, rate = 1)
  y2 <- 1 + rexp(n, rate = 2)
  y3 <- 1 + abs(rnorm(n))

  # Compute importance weights using existing objects
  w1 <- phi(y1) / g1(y1)
  w2 <- phi(y2) / g2(y2)
  w3 <- phi(y3) / g3(y3)

  # Bootstrap the importance sampling estimates
  boot1 <- replicate(B, mean(sample(w1, size = n, replace = TRUE)))
  boot2 <- replicate(B, mean(sample(w2, size = n, replace = TRUE)))
  boot3 <- replicate(B, mean(sample(w3, size = n, replace = TRUE)))

  # Store bootstrap estimates
  boot_estimates <- data.frame(
    g1 = boot1,
    g2 = boot2,
    g3 = boot3
  )

  # Summarize estimates and variances
  summary <- data.frame(
    Importance_Function = c("g1", "g2", "g3"),
    Estimate = c(mean(w1), mean(w2), mean(w3)),
    Bootstrap_Mean = colMeans(boot_estimates),
    Bootstrap_Variance = apply(boot_estimates, 2, var),
    Bootstrap_SE = apply(boot_estimates, 2, sd)
  )

  list(
    summary = summary,
    bootstrap_estimates = boot_estimates
  )
}
results <- bootstrap_IS(n = 10000, B = 1000, seed = 123)

results$summary
true_theta <- 1 - pnorm(1)
paste0("true_theta : ", true_theta)
## [1] "true_theta : 0.158655253931457"
boxplot(
  results$bootstrap_estimates,
  col = c("red", "blue", "darkgreen"),
  names = c("g1", "g2", "g3"),
  ylab = expression(hat(theta)),
  main = "Bootstrap Importance Sampling Estimates"
)

abline(h = true_theta, lty = 2, lwd = 2)

  • notice G2 is the best performing

3 (3)

3.1 (a)

The target density is proportional to

\[ q(x)=e^{-x}\sqrt{x}\sin^2(x), \qquad x>0, \]

so that

\[ f(x)=\frac{q(x)}{C}, \]

where

\[ C=\int_0^\infty q(x)\,dx \]

is an unknown normalizing constant.

We wish to estimate

\[ E_f(X^2) = \frac{\int_0^\infty x^2q(x)\,dx} {\int_0^\infty q(x)\,dx}. \]

Since the normalizing constant is unknown, we use the self-normalized importance sampling estimator

\[ \widehat{E_f(X^2)} = \frac{\sum_{i=1}^{n}Y_i^2w(Y_i)} {\sum_{i=1}^{n}w(Y_i)}, \]

where

\[ w(x)=\frac{q(x)}{r(x)}, \]

and \(Y_i\) are sampled from the normalized version of the proposal density.

3.1.1 Define the target density

q <- function(x) {
  exp(-x) * sqrt(x) * sin(x)^2
}

3.1.2 Proposal 1

The first proposal is

\[ r_1(x)=e^{-2x}. \]

Since

\[ \int_0^\infty e^{-2x}\,dx=\frac12, \]

the normalized proposal is

\[ h_1(x)=2e^{-2x}, \]

which is an exponential distribution with rate \(2\).

The corresponding importance weight is

\[ w_1(x) = \frac{q(x)}{r_1(x)} = e^x\sqrt{x}\sin^2(x). \]

E1 <- function(n) {

  # Sample from Exp(rate = 2)
  y <- rexp(n, rate = 2)

  # Importance weights
  w <- exp(y) * sqrt(y) * sin(y)^2

  # Self-normalized importance sampling estimate
  sum(y^2 * w) / sum(w)
}

3.1.3 Proposal 2

The second proposal is

\[ r_2(x)=x^{-1/2}e^{-x/2}. \]

This is proportional to the density of

\[ Y\sim\text{Gamma}\left(\frac12,\frac12\right). \]

The corresponding importance weight is

\[ w_2(x) = \frac{q(x)}{r_2(x)} = xe^{-x/2}\sin^2(x). \]

E2 <- function(n) {

  # Sample from Gamma(shape = 1/2, rate = 1/2)
  y <- rgamma(n, shape = 0.5, rate = 0.5)

  # Importance weights
  w <- y * exp(-y / 2) * sin(y)^2

  # Self-normalized importance sampling estimate
  sum(y^2 * w) / sum(w)
}

3.1.4 Example

set.seed(123)

E1(10000)
## [1] 3.309329
E2(10000)
## [1] 3.872407

3.2 (b) Comparison with the theoretical value

The theoretical value of \(E_f(X^2)\) is

\[ E_f(X^2) = \frac{\displaystyle\int_0^\infty x^2q(x)\,dx} {\displaystyle\int_0^\infty q(x)\,dx}, \]

where

\[ q(x)=e^{-x}\sqrt{x}\sin^2(x). \]

The denominator is the normalizing constant of the target density, while the numerator is

\[ \int_0^\infty x^2q(x)\,dx = \int_0^\infty x^{5/2}e^{-x}\sin^2(x)\,dx. \]

We can evaluate both integrals numerically using integrate().

# Normalizing constant
normalizing_constant <- integrate(
  q,
  lower = 0,
  upper = Inf
)$value

# Numerator of E_f(X^2)
numerator <- integrate(
  function(x) x^2 * q(x),
  lower = 0,
  upper = Inf
)$value

# Theoretical value
theoretical_value <- numerator / normalizing_constant

theoretical_value
## [1] 3.814188

We now compare the importance sampling estimates from E1() and E2() with the theoretical value.

set.seed(123)

n <- 100000

estimate_E1 <- E1(n)
estimate_E2 <- E2(n)

comparison <- data.frame(
  Method = c(
    "Proposal r1",
    "Proposal r2",
    "Theoretical"
  ),
  Estimate = c(
    estimate_E1,
    estimate_E2,
    theoretical_value
  )
)

comparison

The absolute errors of the two importance sampling estimates are

error_comparison <- data.frame(
  Method = c("Proposal r1", "Proposal r2"),
  Estimate = c(estimate_E1, estimate_E2),
  Absolute_Error = c(
    abs(estimate_E1 - theoretical_value),
    abs(estimate_E2 - theoretical_value)
  )
)

error_comparison

3.2.1 Repeated Monte Carlo comparison

Because a Monte Carlo estimate changes from one random sample to another, a more informative comparison is obtained by repeating each estimator several times.

set.seed(123)

n <- 10000
B <- 1000

estimates_E1 <- replicate(B, E1(n))
estimates_E2 <- replicate(B, E2(n))

simulation_comparison <- data.frame(
  Method = c("Proposal r1", "Proposal r2"),
  Mean_Estimate = c(
    mean(estimates_E1),
    mean(estimates_E2)
  ),
  Bias = c(
    mean(estimates_E1) - theoretical_value,
    mean(estimates_E2) - theoretical_value
  ),
  Variance = c(
    var(estimates_E1),
    var(estimates_E2)
  ),
  MSE = c(
    mean((estimates_E1 - theoretical_value)^2),
    mean((estimates_E2 - theoretical_value)^2)
  )
)

simulation_comparison

The distributions of the estimates can also be compared graphically.

boxplot(
  estimates_E1,
  estimates_E2,
  names = c("r1", "r2"),
  col = c("steelblue", "orange"),
  ylab = expression(hat(E)[f](X^2)),
  main = "Importance Sampling Estimates"
)

abline(
  h = theoretical_value,
  col = "red",
  lty = 2,
  lwd = 2
)

legend(
  "topright",
  legend = "Theoretical value",
  col = "red",
  lty = 2,
  lwd = 2,
  cex = 0.8
)

4 Curiosity :

4.1 Antithetic Variables : U , 1-U reasoning

U1 <- runif(1000)
U2 <- 1-U1

plot(U1, U2)

cor(U1, U2) # notice cor()
## [1] -1

requirements :

  • X, Y are must be from the same underlying dist.

  • X,Y should be negatively correlated to reduce variance

resource : Antithetic Variates + R Demo

5 Computing corr() of same shape – importance Sampling

x <- runif(100)
y <- sin(x)
cor(y, y + 5)
## [1] 1
plot(y, y + 5)

plot(x, y)

plot(x, y + 5)


x <- rnorm(100)
plot(x, sin(x) + rnorm(x, 1, .1))