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)
}

2. 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)

3. Then write functions to evaluate the gradient and Hessian matrix of the log-likelihood for derivatives w.r.t. \(\lambda\).

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)
}

4. Evaluate the gradient and Hessian for the sample generated above with \(\lambda = 2\).

Solution

loglik_d1(2, y0)
loglik_d2(2, y0)

5. 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)

2. 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