Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / c

Taking the logarithm of the normalized probability density function, then summing over the independent observations, gives the log-likelihood

ℓ (p, β; x) = n {lo g β - lo g (- lo g p) + lo g (1 - p)} - β i = 1 \sum n x_{i} - i = 1 \sum n lo g (1 - (1 - p) e^{- β x_{i}}) .

(1)

Both

p

and

β

are parameters being inferred, so their contributions to the normalizing constant must be retained. The following R implementation validates the arguments and returns one number:

ll_f <- function(p, beta, x) {
  scalar <- function(v) {
    is.numeric(v) && !is.complex(v) && length(v) == 1L &&
      is.finite(v)
  }
  if (!scalar(p) || p <= 0 || p >= 1)
    stop("p must lie strictly between 0 and 1")
  if (!scalar(beta) || beta <= 0)
    stop("beta must be positive")
  if (!is.numeric(x) || is.complex(x) || !is.null(dim(x)) ||
      any(!is.finite(x)) || any(x <= 0))
    stop("x must be a vector of finite positive numbers")
  if (length(x) == 0L) return(0)
  t <- log1p(-p) - beta * x
  length(x) * (log(beta) - log(-log(p))) +
    sum(t - log(-expm1(t)))
}

Here

t_{i} = lo g (1 - p) - β x_{i} < 0

and -expm1(t) computes

1 - e^{t}

accurately near zero. The empty-data likelihood function is one, so its log-likelihood is zero. Invalid parameter values produce errors in this evaluator; the next algorithm rejects candidates outside the parameter support before calling it.

ii

 0  0

Solution

 0  0

Assume

a, b, μ > 0

and a current state with

0 < p < 1

β > 0

. The unnormalized posterior density is

π (p, β ∣ x) \propto exp {ℓ (p, β; x)} p^{a - 1} (1 - p)^{b - 1} μ^{- 1} e^{- β / μ} .

(1)

Thus the log-posterior, up to an additive constant, is the log-likelihood plus the beta distribution and exponential distribution log prior densities. The exponential distribution rate is

1/ μ

, since

μ

is its mean. ^[ref] ^[ref]

Choose a symmetric joint random-walk Metropolis algorithm proposal:

p^{'} = p + 0.1 G_{1}

and

β^{'} = β + (μ /2) G_{2}

, for independent standard normal distribution draws. Its scale depends only on fixed hyperparameters, so forward and reverse proposal probability density functions agree. Reject immediately outside

(0, 1) \times (0, \infty)

; otherwise accept with Metropolis–Hastings acceptance probability

α = min {1, exp [lo g π (p^{'}, β^{'} ∣ x) - lo g π (p, β ∣ x)]} .

(2)

The following function makes exactly one proposal, with rejection retaining both old parameters:

MHstep <- function(x, p, beta, a, b, mu) {
  p_new <- rnorm(1, mean = p, sd = 0.1)
  beta_new <- rnorm(1, mean = beta, sd = mu / 2)
  if (p_new <= 0 || p_new >= 1 || beta_new <= 0)
    return(list(p = p, theta = beta))
  logpost <- function(pp, bb) {
    ll_f(pp, bb, x) + dbeta(pp, a, b, log = TRUE) +
      dexp(bb, rate = 1 / mu, log = TRUE)
  }
  logA <- logpost(p_new, beta_new) - logpost(p, beta)
  if (log(runif(1)) <= min(0, logA)) {
    p <- p_new
    beta <- beta_new
  }
  list(p = p, theta = beta)
}

The returned field theta contains the updated value of beta. The original PDF really asks for p and theta, although no parameter

θ

is defined in this question. This choice preserves its literal interface while making the intended second parameter explicit. A caller therefore advances with state <- MHstep(x, p, beta, a, b, mu), p <- state$p, and beta <- state$theta.

The proposal is an unrestricted symmetric normal distribution followed by rejection outside the support. Repeatedly drawing until a valid proposal appears would instead truncate and renormalize the proposal, requiring a proposal-density correction. The code uses immediate rejection and needs no such correction. The displayed scales are valid fixed choices; practical efficiency can be improved using a separate tuning phase.

iii

 0  0

Solution

 0  0

Use several complementary Markov chain Monte Carlo convergence diagnostics, monitoring both

p

and

β

rather than judging the chain solely by its length:

Trace plots of each parameter and the log-posterior against iteration, including early iterations and later windows. Look for initial drift, long flat stretches, abrupt switches and sustained changes in spread. Several independent chains started at dispersed points should explore the same region after burn-in.
Plots of cumulative sample means and selected sample quantiles against iteration. Compare chains and successive long blocks; persistent movement indicates that the estimated posterior moments or tails have not stabilized.
Autocorrelation function plots for both parameters after tentative burn-in. Slow decay reveals strong serial dependence and a small effective sample size of a Markov chain, even when $1 0^{6}$ iterations have been generated.
Overlaid histograms or kernel density estimation curves from successive retained blocks and from different chains. Their agreement provides evidence of stable marginal posterior distributions.
A scatter plot of retained $(p, β)$ pairs, colored by time block if helpful. Strong dependence or separate clusters can expose slow joint exploration that marginal plots conceal.

A million iterations and stable-looking plots do not prove convergence. Agreement across dispersed chains, exploration of all visible modes and adequate effective sample size of a Markov chain provide stronger evidence. Estimate the Monte Carlo standard error using the chain's autocorrelation; the independent random variables formula from 3(b)(ii) does not apply directly to these dependent draws.

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / c

i

Solution

ii

Solution

iii

Solution

 Ancestors (11)

 Discussion (0)

 Articles by others on the same topic (0)

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / c

 Discussion (0)  Subscribe (1)

 Discussion (0)