Chapter 4

Week 6 lecture 2

Challenges I

  1. Approximate the derivative of \(f(x) = \exp(-x^2)\) by finite-differencing, and compare your approximation to the true derivative over \(x \in [-1, 1]\).
Solution 
x <- seq(-1, 1, by = .1)
deriv0 <- - 2 * x * exp(-x^2)
plot(x, deriv0, type = 'l')
delta1 <- sqrt(.Machine$double.eps)
f <- function(x) exp(-x^2)
deriv1 <- (f(x + delta1) - f(x)) / delta1
lines(x, deriv1, col = 2, lwd = 3, lty = 2)

Challenges II

  1. Consider the exponential distribution with pdf \[ f(y \mid \lambda) = \lambda \exp(-\lambda y) \] for \(y, \lambda > 0\). Its log-likelihood, for an independent sample \(\mathbf{y} = (y_1, \ldots, y_n)\) is given by \[ \log f(\mathbf{y} \mid \lambda) = n \log \lambda -\lambda \sum_{i = 1}^n y_i. \]

Write a function, loglik(lambda, y), to evaluate the log-likelihood above, where \(\lambda =\) lambda and \(\mathbf{y} =\) y.

Solution 
loglik <- function(lambda, y) {
  # function to evaluate exponential log-likelihood
  # lambda is a scalar
  # y is a vector
  # returns a scalar
  n <- length(y)
  n * log(lambda) - lambda * sum(y)
}
  1. Use rexp(10) to generate a sample of ten values from the Exp(1) distribution, and then use loglik(lambda, y) to evaluate \(\log f(\mathbf{y} \mid 2)\) for the sample.
Solution 
y0 <- rexp(10, 1)
loglik(2, y0)
  1. Then write functions to evaluate the gradient and Hessian matrix of the log-likelihood.
Solution 
# gradient, i.e. first derivative
loglik_d1 <- function(lambda, y) {
  # function to evaluate gradient of exponential log-likelihood w.r.t. lambda
  # lambda is a scalar
  # y is a vector
  # returns a scalar
  n <- length(y)
  n / lambda - sum(y)
}

# Hessian, i.e. second derivative
loglik_d2 <- function(lambda, y) {
  # function to evaluate Hessian of exponential log-likelihood w.r.t. lambda
  # lambda is a scalar
  # y is a vector
  # returns a 1 x 1 matrix
  n <- length(y)
  matrix(- n / lambda^2, 1, 1)
}
  1. Evaluate the gradient and Hessian for the sample generated above with \(\lambda = 2\).
Solution 
loglik_d1(2, y0)
loglik_d2(2, y0)
  1. Check your gradient calculation against an approximation based on finite-differencing.
Solution 
(loglik(2 + 1e-6, y0) - loglik(2, y0)) / 1e-6
loglik_d1(2, y0)

Week 6 lecture 3

Challenges I

  1. Use the midpoint rule to approximate \[ I = \int_0^2 \text{exp}\left(-\frac{x^2}{2}\right) \text{d}x \] with \(N = 10\) integration nodes.
Solution 
N <- 10
a <- 0
b <- 2
bins <- seq(a, b, length = N + 1)
h <- bins[2] - bins[1]
mids <- bins[-1] - .5 * h
f <- function(x) exp(-x^2/2)
f_mids <- f(mids)
I_hat <- h * sum(f_mids)
  1. Find the relative absolute error of your estimate. [Hint: notice the resemblance of the integrand to the Normal(0, 1) distribution’s pdf.]
Solution 
I0 <- sqrt(2 * pi) * (pnorm(2) - .5)
rel_err <- abs((I0 - I_hat) / I0)
100 * rel_err

Week 8 lecture 1

Challenges I

  1. Repeat Challenge I using Simpson’s rule; i.e. approximate \[ I = \int_0^2 \text{exp}\left(-\frac{x^2}{2}\right) \text{d}x \] with \(N = 10\) and find the relative absolute error of your estimate.
Solution 
f <- function(x) exp(-x^2/2)
N <- 10
a <- 0
b <- 2
h <- (b - a) / N
simpson <- f(a) + f(b)
x1i <- a + h * (2 * (1:N) - 1) / 2
simpson <- simpson + 4 * sum(f(x1i))
x2i <- a + h * (1:(N - 1))
simpson <- simpson + 2 * sum(f(x2i))
simpson <- h * simpson / 6

true <- sqrt(2 * pi) * (pnorm(b) - .5)
rel_err <- abs((true - simpson) / true)
rel_err

Challenges II

  1. Use Gaussian quadrature to approximate \[ I = \int_0^2 \text{exp}(-x^2 / 2) \text{d} x \] with \(N = 7\) integration nodes.
Solution 
f <- function(x) exp(-x^2/2)
N <- 7
gq <- pracma::gaussLegendre(N, 0, 2)
I_hat <- sum(gq$w * f(gq$x))
  1. Calculate the relative absolute error of your approximation.
Solution 
true <- sqrt(2 * pi) * (pnorm(2) - .5)
rel_err <- abs((true - I_hat) / true)

Week 8 lecture 2

Challenges I

  1. Consider independent Exponential random variables \(X_1\) and \(X_2\), each with pdf \[ f_{X_i}(x \mid \lambda) = \lambda \text{exp}(-\lambda x)~~~\text{for}~x > 0, \] and where \(\lambda > 0\). Estimate \(\text{Pr}(X_1 \leq 1, X_2 \leq 1)\) for \(\lambda = 2\) by numerical integration using the midpoint rule and \(N = 10\) for \(X_1\) and \(X_2\); i.e. approximate \[ \text{Pr}(X_1 \leq 1, X_2 \leq 1) = \int_0^1 \int_0^1 \lambda^2 \text{exp}\left\{-\lambda(x_1 + x_2)\right\} \text{d}x_1\text{d}x_2. \]
Solution 
N <- 10
x1 <- x2 <- (1:N - .5) / N
h <- 1/N
w <- h^2
lambda <- 2
midpoint <- 0
for (i in 1:N) for (j in 1:N)
  midpoint <- midpoint + w * lambda^2 * exp(-lambda*(x1[i] + x2[j]))
midpoint
  1. Find the relative absolute error of your approximation.
Solution 
true <- pexp(1, lambda)^2
rel_err <- abs((true - midpoint) / true)

Challenges II

  1. Approximate \(\text{Pr}(X_1 \leq 1, X_2 \leq 2)\), where \(X_1\) and \(X_2\) are independent Exp(\(\lambda\)) random variables with \(\lambda = 2\). This time use Gaussian quadrature and \(N = 7\) for \(X_1\) and \(X_2\).
Solution 
N <- 7
xw1 <- pracma::gaussLegendre(N, 0, 1)
xw2 <- pracma::gaussLegendre(N, 0, 2)
gq <- 0
f <- function(x1, x2, lambda = 2)
  lambda^2 * exp(-lambda * (x1 + x2)) 
for (i in 1:N) for (j in 1:N) 
  gq <- gq + xw1$w[i] * xw2$w[j] * f(xw1$x[i], xw2$x[j])
gq
  1. Find the relative absolute error of your approximation.
Solution 
true <- pexp(1, 2) * pexp(2, 2)
rel_err <- abs((true - gq) / true)

Week 8 lecture 3

Challenge I

  1. Using the midpoint rule with \(N = 5\) nodes per variable, approximate \[ I = \int_1^3 \ldots \int_1^3 \prod_{i = 1}^4 \left\{1 + (x_i - 2)^2\right\} \text{d}x_1 \ldots \text{d}x_4 \] and estimate its relative absolute error.
Solution 
N <- 5
d <- 4
b <- 3
a <- 1
h <- (b - a) / N
x <- a + h * (1:N - .5)
xx <- lapply(1:d, function(i) x)
X <- expand.grid(xx)
w <- h^d
f <- apply(1 + (X - 2)^2, 1, prod)
midpoint <- w * sum(f)
true <- (8/3)^d
rel_err <- abs((true - midpoint) / true)
rel_err

Challenge II

  1. Consider \(Y \sim \text{Exp}(\lambda)\), so that its pdf is \[ f(y \mid \lambda) = \lambda \text{exp}(-\lambda y) \quad \text{for}~y > 0, \] where \(\lambda\) is a random variable with pdf \[ f(\lambda) = 4 \lambda \text{exp}(-2 \lambda) \quad \text{for}~\lambda > 0. \] Use Laplace’s method to find \(f(y)\), where \[ f(y) = \int_0^\infty f(y \mid \lambda) f(\lambda) \text{d} \lambda. \]
Solution

  • We first notice that, given Equation (4.4), we take \(f()\) to be \(-\log f(y, \lambda)\) and \(x\) to be \(\lambda\), as that’s what we’re integrating out.

  • Then we get \(f(y, \lambda) = f(y \mid \lambda) f(\lambda) = \lambda \text{e}^{-\lambda y} \times 4 \lambda \text{e}^{-2 \lambda} = 4 \lambda^2 \text{e}^{-\lambda(2 + y)}\)

  • So \(-\log f(y, \lambda) = -2 \log \lambda + \lambda(2 + y) - \log 4\)

  • Then \(-\dfrac{\partial \log f(y, \lambda)}{\partial \lambda} = -\dfrac{2}{\lambda} + (2 + y)\), which is minimised at \(\tilde \lambda = 2 / (2 + y)\)

  • Also \(-\dfrac{\partial^2 \log f(y, \lambda)}{\partial \lambda^2} = 2 / \lambda^2\) and \(2 / \tilde \lambda^2 = 2 / (4 / (2 + y)^2) = (2 + y)^2 / 2\)

  • So, if \(g''(\lambda)\) denotes \(-\partial^2 \log f(y, \lambda) / \partial \lambda^2\), then \[\begin{align*} I_n &= f(y, \tilde \lambda) \sqrt{\dfrac{2 \pi}{g''(\tilde \lambda)}}\\ &= 4 \hat \lambda^2 \text{e}^{-\hat \lambda(2 + y)} \sqrt{2\pi \times 2(2 + y)^{-2}} = \dfrac{32 \sqrt{\pi}}{\text{e}^2(2 + y)^3} \end{align*}\]

Week 9 lecture 1

Challenges I

  1. Use Monte Carlo integration with \(N = 10^3\) Monte Carlo samples to approximate \[ \int_{-1}^1 x^2 \text{d} x \] and then calculate the relative absolute error of your approximation.
Solution 
N <- 1e3
V <- 2
mc <- V * mean(runif(N, -1, 1)^2)
true <- 2/3

Challenges II

  1. Approximate \[ I = \int_1^3 \ldots \int_1^3 \prod_{i = 1}^4 \left\{1 + (x_i - 2)^2\right\} \text{d}x_1 \ldots \text{d}x_4 \] using Monte Carlo sampling and a Monte Carlo sample size of \(N = 10^4\) and estimate the relative absolute error of your approximation.
Solution 
N <- 10^4
d <- 4
x <- matrix(runif(N * d, 1, 3), N)
V <- 2^d
I_hat <- V * mean(apply(1 + (x - 2)^2, 1, prod))
I_hat2 <- V * mean(exp(rowSums(log1p((x - 2)^2))))
true <- (8/3)^d
rel_err <- abs((true - I_hat) / true)
rel_err