Problem 1

We are to evaluate the following integral:

\[ \int_{0}^{1}\frac{1}{1+x^2}e^{-x}dx \]

Noting that \(f_X(x) = e^{-x}\) is the density of an exponential random variable with rate \(\lambda=1\) implies that the above integral can be written as the \(E[\frac{1}{(1+X^2)}]\) where \(X \sim exponential(\lambda=1)\).

Define \(h(x) = \frac{1}{(1+x^2)}\). Then we are trying to estimate \(E[h(X)]\) with \(X \sim exponential(\lambda=1)\).

We will use the method of importance sampling in parts (a) and (b) to estimate this quantity using different important functions \(g\).

Problem 1 Part A

We use \(g(x) = \frac{1}{\pi(1+x^2)}\) with \(-\infty < x < \infty\). In other words \(g\) is the density of a standard Cauchy.

Hence, with \(Y_i \sim Cauchy\) we have our estimator.

\[ \hat{\theta}_{IS} = \frac{1}{n}\sum_{i=1}^nh(Y_i)\frac{f(Y_i)}{g(Y_i)} \]

n = 20000
y = rcauchy(n)
gy = dcauchy(y)
fy = dexp(y, rate=1)
hy = 1/(1+y^2)
theta.ICA = (1/n)*sum(hy*(fy/gy))
SE.theta.ICA = sd(hy*(fy/gy))/(sqrt(n))
list.result = list(theta.ISA = theta.ICA, SE.theta.ISA = SE.theta.ICA)
list.result
## $theta.ISA
## [1] 0.6114511
## 
## $SE.theta.ISA
## [1] 0.006558434

Problem 1 Part B

We use \(g(x) = \frac{4}{\pi(1+x^2)}\) with \(0 \leq x \leq 1\). In other words \(g\) is the density of a truncated Cauchy.

Hence, with \(Y_i \sim truncated Cauchy\) we have our estimator.

\[ \hat{\theta}_{IS} = \frac{1}{n}\sum_{i=1}^nh(Y_i)\frac{f(Y_i)}{g(Y_i)} \]

G.inv = function(x) {
  return(tan(pi*x/4))
}

g.trunCauchy = function(x) {
  denom = pi*(1+x^2)
  return(4/denom)
}

u = runif(n)
y = G.inv(u)
gy = g.trunCauchy(y)
fy = dexp(y, rate=1)
hy = 1/(1+y^2)
theta.ICB = (1/n)*sum(hy*(fy/gy))
SE.theta.ICB = sd(hy*(fy/gy))/(sqrt(n))
list.result = list(theta.ISB = theta.ICB, SE.theta.ISB = SE.theta.ICB)
list.result
## $theta.ISB
## [1] 0.5253039
## 
## $SE.theta.ISB
## [1] 0.0009963577

Our estimate using the truncated Cauchy is much more accurate and has lower variability then our estimate using the standard Cauchy.

Problem 2

We are to estimate that probability that of our random point in 3 dimensions falls in some set \(S\).

We define some notation.

\[ \theta = P(X \in S) = E[1_{X\in S}] \]

Hence we are trying to find \(E[h(X)]\) where \(h(X) = 1_{X \in S}\) and $X $ trivariate \(t\).

Using \(n=20,000\) random variates from an appropriate tri-variate normal and using the importantce sampling method (using standardized sampling weights) we have the following estimator:

\[ \hat{\theta}_{ISS} = \frac{\sum_{i=1}^nh(Y_i)w^{*}(Y_i)}{\sum_{i=1}^nw^{*}(Y_i)} \]

where \(w^{*}(Y_i) = \frac{f(Y_i)}{g(Y_i)}\)

MU = c(0,0,0)
SIGMA = matrix(c(1, 3/5, 1/3, 3/5, 1, 11/15, 1/3, 11/15, 1), nrow=3)
SIGMA.INV = solve(SIGMA)


f.tdist <- function(x) {
  x = as.matrix(x, ncol=1)
  val = 5 + t(x) %*% SIGMA.INV %*% x
  return(val^-4)
}

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

g.norm <- function(x) {
  return(dmvnorm(x, MU,SIGMA))
}

#-------
library(MASS)
library(mvtnorm)
n = 20000
Y = mvrnorm(n, MU, SIGMA) #trivariate normal in matrix

f.y = apply(Y, 1, f.tdist) 
g.y = apply(Y, 1, g.norm)
h.y = apply(Y, 1, h.function)
w.star = f.y/g.y
w = sum(w.star)

theta.ISS = sum(h.y * w.star)/w
theta.ISS
## [1] 0.7958699

Our estimate is close to the true value of 0.79145379

Problem 3 Part A

As we can see from the plot below, the double-exponential and standard normal intersect at a value of \(y\) greater than 0 and less than 0.

vals = seq(-3, 3, length.out=1000)
plot(1, type='n', xlab='y', xlim=c(-3,3), ylim=c(0,1),  ylab='Density')
lines(vals, exp(-abs(vals))/2, col='blue')
lines(vals, dnorm(vals), col='red')
legend('topleft', legend=c("Double-Exp", "Stand Norm"), text.col=c('blue', 'red'))

Problem 3 Part B

We construct an envelope \(e(x) = \frac{g(x)}{\alpha}\) where \(g \sim\) Standard Normal. A property of envelopes that we will utilize is the fact that \(e(x) > f(x)\) for all \(x\) values in the support. We first seek the \(\alpha\) between 0 and 1 that has this property.

library(nimble)
## nimble version 0.11.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit http://R-nimble.org.
## 
## Attaching package: 'nimble'
## The following object is masked from 'package:stats':
## 
##     simulate
search.alpha = function(alpha, x) {
  gx = exp(-abs(x))/2
  fx = dnorm(x)
  bool = any(gx/alpha < fx)
  while(bool & alpha > 0) {
    alpha = alpha - 0.01
    bool = any(gx/alpha < fx)
  }
  return(alpha)
}
x = seq(-3, 3, length.out=1000)
alpha = 1
feasible.alpha = search.alpha(alpha, x)
feasible.alpha
## [1] 0.76

Thus, the \(\alpha\) we ultimately choose must be less than or equal to 0.76. If \(\alpha\) is larger than 0.76 than we don’t have \(e(x) > f(x)\) for all \(x\) values in the support.

Now we search for the \(\alpha\) that minimizes the number of rejections.

alpha.possible = seq(0.01,feasible.alpha, by=0.01)
rejects = rep(-1, length(alpha.possible))
n = 10000
U = runif(n)
Y = rdexp(n)
fy = dnorm(Y)
for(i in 1:length(alpha.possible)) {
  ey = ddexp(Y)/alpha.possible[i]
  rejects[i] = length(Y[U > (fy/ey)])
}
plot(alpha.possible, rejects)

best.alpha = alpha.possible[which.min(rejects)]
best.alpha
## [1] 0.76

As we can see, the value of \(\alpha\) that will minimize the number of rejections is \(\alpha=0.76\). This can be intuitively understood by the following plot.

vals = seq(-3, 3, length.out=1000)
plot(1, type='n', xlab='y', xlim=c(-3,3), ylim=c(0,1),  ylab='Density')
lines(vals, ddexp(vals)/best.alpha, col='blue')
lines(vals, ddexp(vals)/0.50, col='green')
lines(vals, dnorm(vals), col='red')
legend('topleft', legend=c("Double-Exp/0.76", "Stand Norm", "Double-Exp/0.50"), text.col=c('blue', 'red', 'green'))

We can see that any \(\alpha\) less than 0.76 will just continue to raise the density higher and higher. This will lead to more of a gap between \(g\) and \(f\) where \(f\) is the red standard normal.

Problem 3 Part C

Our algorithm will be the following:

1.) Generate $Y $ Double-Exponential(\(\lambda=1\))
2.) Generate $U $ Unif(0,1)
3.) Reject \(Y\) if \(U > \frac{f_X(Y)}{e(Y)}\) where $f_X(x) $ Standard Normal
3b.) Otherwise accept \(Y\) and set \(X = Y\)

n = 10000
U = runif(n)
Y = rdexp(n)
fy = dnorm(Y)
ey = ddexp(Y)/0.76
X = Y[U <= (fy/ey)]
par(mfrow=c(1,2))
hist(X)
qqnorm(X)