Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 2

Table of contents
- a 2
  - Solution a
- b 2
  - i b
    - Solution i
  - ii b
    - Solution ii
- c 2
  - Solution c
- d 2
  - Solution d
- e 2
  - Solution e
- f 2
  - Solution f

f (x_{i} ∣ μ) = \frac{1}{2 π} exp [- \frac{1}{2} (x_{i} - μ)^{2}] .

(1)

The observations are independent and identically distributed, so their joint probability density function, viewed as a function of the unknown mean, is the product:

L (μ) = i = 1 \prod n f (x_{i} ∣ μ) = (2 π)^{- n /2} exp [- \frac{1}{2} i = 1 \sum n (x_{i} - μ)^{2}] .

(2)

Independence justifies multiplication; the common mean and known unit variance justify using this same density factor for every observation.

Solution

 0  0

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

The estimate requested is the sample mean, $\overset{x}{ˉ} = n^{- 1} i \sum x_{i}$ . The named electronic file Iris_data.csv is not present in the supplied source archive or this repository, and its measurements are not printed in the paper. Consequently its numerical sample mean cannot be recovered from the paper alone. The familiar built-in iris data must not be silently substituted: the question does not identify those observations as its sample.

The following R code reads the actual file and prints its sample mean. The later code blocks in Question 2 continue this same script and use the resulting x, n and xbar.

dat <- read.csv("Iris_data.csv", check.names = FALSE)
# Recognize common labels, or a single measurement column.
keys <- tolower(gsub("[^[:alnum:]]", "", names(dat)))
j <- which(keys %in% c("sepallength", "sepallengthcm", "x"))
if (length(j) == 0L && ncol(dat) == 1L) j <- 1L
if (length(j) != 1L) {
  stop("Select the sepal-length column explicitly after inspecting names(dat).")
}
x <- dat[[j]]
stopifnot(is.numeric(x), length(x) > 0L, all(is.finite(x)))
n <- length(x)
xbar <- mean(x)
cat("n =", n, " sample mean =", format(xbar, digits = 12), "\n")

This reader assumes a header row and one observation per row. An alternative layout should be identified explicitly rather than dropping measurements or treating an identifier column as a length.

ii

 0  0

Solution

 0  0

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

The identity

\sum_{i} (x_{i} - μ)^{2} = \sum_{i} (x_{i} - \overset{x}{ˉ})^{2} + n (μ - \overset{x}{ˉ})^{2}

(1)

shows that the likelihood function is a bell-shaped curve in

μ

, maximized at the sample mean, with width

1/ n

. Its actual vertical scale also depends on the within-sample sum of squares.

Continue the R script with:

sse <- sum((x - xbar)^2)
loglik <- function(mu) {
  -n/2 * log(2*pi) - 0.5 * (sse + n * (mu - xbar)^2)
}
grid <- seq(3, 6, length.out = 2001L)
ell <- loglik(grid)
png("q2-likelihood.png", width = 900, height = 600)
par(mar = c(4.5, 4.5, 4, 1))
if (max(ell) > log(.Machine$double.xmin)) {
  plot(grid, exp(ell), type = "l", xlab = expression(mu),
       ylab = expression(L(mu)), main = "Likelihood on [3,6]")
} else {
  # Explicitly label a proportional likelihood if the true scale underflows.
  plot(grid, exp(ell - max(ell)), type = "l", xlab = expression(mu),
       ylab = "Likelihood relative to maximum on plotted grid",
       main = "Likelihood on [3,6] (rescaled)")
}
abline(v = xbar, col = "steelblue", lty = 2)
dev.off()

Working first with the log-likelihood avoids premature numerical underflow. The fallback preserves the shape on the requested interval and labels the altered scale. A data-specific plot requires the missing CSV file; no fabricated observation set is used here.

c

 0  0

Solution

 0  0

The exponential distribution prior has rate

λ = 1/3

and support

μ \geq 0

, so its prior density is

(1/3) e^{- μ /3} 1_{{μ \geq 0}}

. By Bayes' theorem,

p (μ ∣ x) = C exp [- \frac{1}{2} i \sum (x_{i} - μ)^{2} - \frac{μ}{3}] 1_{{μ \geq 0}},

(1)

where

C

is a normalizing constant independent of

μ

Completing the square gives a useful exact reference for the numerical sampling. Put

m = \overset{x}{ˉ} - \frac{1}{3 n}, s = \frac{1}{n} .

(2)

Then this normal mean posterior with an exponential prior is a normal distribution with location

m

and variance

s^{2}

, truncated below at zero. With

ϕ

the standard normal distribution density and

Φ

its cumulative distribution function,

p (μ ∣ x) = \frac{ϕ (( μ - m ) / s )}{s Φ ( m / s )} 1_{{μ \geq 0}} .

(3)

The exponential rate is

1/3

, not three; the prior mean is three. The support restriction must be retained in both the formula and the sampler.

d

 0  0

Solution

 0  0

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Use the random-walk Metropolis algorithm with proposal

μ^{'} = μ + ε

, where

ε \sim N (0, 1)

. The proposal is symmetric, so for a nonnegative proposal the acceptance probability is

a (μ, μ^{'}) = min {1, exp [ℓ_{p} (μ^{'}) - ℓ_{p} (μ)]}, ℓ_{p} (u) = - \frac{1}{2} \sum_{i} (x_{i} - u)^{2} - \frac{u}{3} .

(1)

Negative proposals are rejected. They are not redrawn until positive, since that would change the proposal kernel and require a different acceptance ratio. On rejection, store the current state again; keeping accepted moves alone would bias the empirical posterior distribution.

This R continuation makes exactly 10,000 proposals, starting at

μ_{0} = 4

set.seed(20240528)
B <- 10000L
mu <- numeric(B + 1L)
mu[1L] <- 4
logpost <- function(u) {
  if (!is.finite(u) || u < 0) return(-Inf)
  loglik(u) - u/3
}
lp <- logpost(mu[1L])
accepted <- 0L
for (i in seq_len(B)) {
  proposal <- rnorm(1L, mean = mu[i], sd = 1)
  lp_new <- logpost(proposal)
  if (log(runif(1L)) < min(0, lp_new - lp)) {
    mu[i + 1L] <- proposal
    lp <- lp_new
    accepted <- accepted + 1L
  } else {
    mu[i + 1L] <- mu[i]
  }
}
cat("Acceptance fraction =", accepted/B, "\n")
# A provisional discard; assess it with the trace in part (e).
burn <- 1000L
post <- mu[(burn + 2L):(B + 1L)]
m <- xbar - 1/(3*n)
s <- 1/sqrt(n)
posterior_density <- function(u) {
  ans <- numeric(length(u))
  ok <- u >= 0
  ans[ok] <- exp(dnorm(u[ok], m, s, log = TRUE) -
                   pnorm(m/s, log.p = TRUE))
  ans
}
png("q2-posterior.png", width = 900, height = 600)
hist(post, breaks = 40, probability = TRUE, xlab = expression(mu),
     main = "Retained Metropolis draws and exact posterior")
u <- seq(0, max(post) + 3*s, length.out = 4001L)
lines(u, posterior_density(u), col = "firebrick", lwd = 2)
dev.off()

The exact truncated normal distribution density overlays the histogram as a check; it does not replace the requested Markov chain Monte Carlo procedure. The 9,000 retained states exclude the initial state and the first 1,000 transitions. The particular histogram and acceptance fraction remain dependent on the missing measurements.

e

 0  0

Solution

 0  0

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Plot all 10,001 stored states, including the initial value, and indicate the provisional burn-in boundary. Also display the retained chain's autocorrelation function:

png("q2-trace.png", width = 1000, height = 700)
par(mfrow = c(2, 1), mar = c(4, 4.5, 2.5, 1))
plot(0:B, mu, type = "l", xlab = "Iteration", ylab = expression(mu),
     main = "Complete Metropolis trace")
abline(v = burn, col = "firebrick", lty = 2)
abline(h = m, col = "steelblue", lty = 3)
acf(post, main = "Autocorrelation after provisional burn-in")
dev.off()
# Check sensitivity to the chosen discard and to the retained time window.
for (discard in c(500L, 1000L, 2000L)) {
  kept <- mu[(discard + 2L):(B + 1L)]
  cat("Discard", discard, ": posterior mean", mean(kept), "\n")
}
cat("Means of retained halves:", mean(post[1:4500]),
    mean(post[4501:9000]), "\n")

A satisfactory trace plot should show an initial transient, if any, followed by repeated exploration of one stable range with no continuing drift. Repeated values are expected because rejections are stored. If

n

is large, the prescribed proposal standard deviation one can be much wider than the posterior scale

1/ n

, giving low acceptance and conspicuous plateaus. A flat-looking trace is therefore not by itself evidence of good exploration.

Compare the retained histogram with the exact posterior density, inspect autocorrelation, compare separate retained windows, and assess sensitivity to the discarded count. Agreement supports the adequacy of this run; neither 10,000 iterations nor a chosen burn-in of 1,000 proves convergence. The effective sample size of a Markov chain can be much smaller than the number of stored states.

The actual trace cannot be diagnosed without Iris_data.csv. The code produces the required trace when the measurements are available; the comments here specify what to examine rather than asserting an unseen chain has converged.

f

 0  0

Solution

 0  0

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Estimate the posterior mean by averaging the retained states:

E [μ ∣ x] = \frac{1}{9000} t = 1001 \sum 10000 μ_{t} .

(1)

The exact truncated normal distribution reference is

E [μ ∣ x] = m + s \frac{ϕ ( m / s )}{Φ ( m / s )}, m = \overset{x}{ˉ} - \frac{1}{3 n}, s = n^{- 1/2} .

(2)

For a posterior well separated from zero, the truncation correction is tiny and the posterior mean is approximately

\overset{x}{ˉ} - 1/ (3 n)

, slightly below the sample mean.

The R continuation prints the Monte Carlo estimator, its approximate batch-based Monte Carlo standard error, and the exact value for comparison:

posterior_mean_mc <- mean(post)
posterior_mean_exact <- m + s * exp(dnorm(m/s, log = TRUE) -
                                    pnorm(m/s, log.p = TRUE))
# Nonoverlapping batches account approximately for serial dependence.
batch_length <- floor(sqrt(length(post)))
number_batches <- floor(length(post)/batch_length)
used <- post[seq_len(batch_length * number_batches)]
batch_averages <- colMeans(matrix(used, nrow = batch_length))
mcse <- sd(batch_averages)/sqrt(number_batches)
cat("Posterior mean, MCMC =", posterior_mean_mc,
    " approximate MCSE =", mcse, "\n")
cat("Posterior mean, exact =", posterior_mean_exact, "\n")

The batch length must be long enough relative to the chain's dependence for this uncertainty estimate to be reliable. No numerical data-specific posterior mean is determined by the supplied paper alone.

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 2

a

Solution

b

i

Solution

ii

Solution

c

Solution

d

Solution

e

Solution

f

Solution

 Ancestors (11)

 Discussion (0)

 Articles by others on the same topic (0)

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 2

 Discussion (0)  Subscribe (1)

 Discussion (0)