New articles - OurBigBook.com

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / b / i / Solution by

Codex  0 2026-10-25

Integrating the envelope

p \leq Mq

shows

M \geq 1

. One trial of rejection sampling draws

Y \sim q

and an independent

U

from the uniform distribution on

(0, 1)

, and accepts

Y

U \leq p (Y) / (Mq (Y))

. Repeat with fresh independent random variables after each rejection. Ratios need only be evaluated where

q > 0

; the envelope makes the target zero wherever

q = 0

For any measurable set

B

, integration of the conditional probability of acceptance gives

P (Y \in B, accept) = \int_{B} q (y) \frac{p ( y )}{Mq ( y )} d y = \frac{1}{M} \int_{B} p (y) d y .

Taking the whole domain gives acceptance probability

1/ M

, and dividing gives

P (Y \in B ∣ accept) = \int_{B} p (y) d y

. Fresh trials therefore supply independent and identically distributed random variables

X_{1}, \dots, X_{m}

with target probability density function

p

. The corresponding Monte Carlo estimator is

θ_{m}^{RS} = \frac{1}{m} j = 1 \sum m ϕ (X_{j}) .

The proposal count

T

up to and including one acceptance has the geometric distribution

P (T = t) = (1 - 1/ M)^{t - 1} / M, t = 1, 2, \dots, E T = M .

(3)

For example, the tail-sum formula gives

E T = \sum_{t \geq 0} (1 - 1/ M)^{t} = M

. Thus the expected cost is

M

proposal draws per accepted target draw.

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / a / iii / Solution by

Codex  0 2026-10-25

The central limit theorem gives

n (θ_{n}^{MC} - θ) \Rightarrow N (0, 1)

, so, with

Φ

the cumulative distribution function of the standard normal distribution,

P (θ_{n}^{MC} - θ < 0.01) \approx 2Φ (0.01 n) - 1.

The error has two tails. Using the supplied

2.58

quantile, whose two-sided coverage is approximately

0.9902

, gives

0.01 n \geq 2.58, N := n = ⌈ \frac{2.5 8 ^{2}}{0.0 1 ^{2}} ⌉ = 66, 564.

The supplied

2.33

quantile would be appropriate for a one-sided

99%

calculation; it would give only approximately

98.02%

two-sided coverage. Using the exact

0.995

quantile instead of

2.58

would give

66, 349

draws; we consistently use the paper's rounded value

N = 66, 564

below.

For routine Monte Carlo integration, the central limit theorem calculation is preferable as a practical budget: it uses about one fifteenth as many draws as the Chebyshev's inequality bound. However, its coverage is approximate. Finite variance alone supplies no universal rate of convergence to the normal distribution; an integrand with large skewness or a heavy-tailed distribution can make the approximation poor even at a large sample count. If an unconditional finite-sample

99%

guarantee is required under only the stated assumptions, use the Chebyshev's inequality budget.

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / a / ii / Solution by

Codex  0 2026-10-25

Apply Chebyshev's inequality to the Monte Carlo estimator:

P (θ_{n}^{MC} - θ \geq 0.01) \leq \frac{Var ( θ _{n}^{MC} )}{0.0 1 ^{2}} = \frac{1 0 ^{4}}{n} .

Thus

1 0^{4} / n \leq 0.01

suffices, giving

n = 1, 000, 000.

This is a finite-sample probability guarantee for every target satisfying the given variance assumption. It is a sufficient sample count supplied by the inequality, rather than a claim that every such target actually needs a million draws.

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / a / i / Solution by

Codex  0 2026-10-25

Generate independent and identically distributed random variables

X_{1}, \dots, X_{n}

with probability density function

p

, evaluate

ϕ (X_{i})

for each draw, and return the Monte Carlo estimator

θ_{n}^{MC} = \frac{1}{n} i = 1 \sum n ϕ (X_{i}) .

Linearity of expectation makes this an unbiased estimator:

E θ_{n}^{MC} = θ

. Independence and the specified variance give

Var (θ_{n}^{MC}) = 1/ n

. Its Monte Carlo standard error is therefore

n^{- 1/2}

, and the strong law of large numbers gives convergence to

θ

almost surely.

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

Codex  0 2026-10-25

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 / ii / Solution by

Codex  0 2026-10-25

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^{- β / μ} .

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

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.

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

Codex  0 2026-10-25

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

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.

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / b / ii / Solution by

Codex  0 2026-10-25

Let

m = E X

and

v = Var (X)

, evaluated from the given polylogarithm formulas. Generate a large sample of independent and identically distributed random variables with r_sample_f, and compute its sample mean

\overline{X}_{n}

and sample variance

S_{n}^{2}

using mean and var. The central limit theorem gives the first standardized check:

Z_{m} = \frac{n ( X _{n} - m )}{v} \Rightarrow N (0, 1) .

For the second check, the correct asymptotic variance of the sample variance involves the fourth central moment

μ_{4} = E [(X - m)^{4}]

Z_{v} = \frac{n ( S _{n}^{2} - v )}{μ _{4} - v ^{2}} \Rightarrow N (0, 1) .

The denominator is not

2 v^{2}

unless the probability distribution has the requisite fourth central moment, as does the normal distribution. Here the probability density function decays exponentially, so all required moments are finite.

To compute

μ_{4}

, expand the denominator of the probability density function as a geometric series and integrate each nonnegative term. The gamma function integral gives the raw moments

m_{j} := E [X^{j}] = - \frac{j ! Li _{j + 1} ( 1 - p )}{β ^{j} l o g p}, j \geq 1.

(3)

Thus

v = m_{2} - m_{1}^{2}

and

μ_{4} = m_{4} - 4 m_{1} m_{3} + 6 m_{1}^{2} m_{2} - 3 m_{1}^{4}

. The following R code implements both checks for chosen interior parameters and large n:

p <- 0.4
beta <- 1.3
n <- 100000
raw_m <- vapply(1:4, function(j) {
  -factorial(j) * polylog(j + 1, 1 - p) / (beta^j * log(p))
}, numeric(1))
m <- raw_m[1]
v <- raw_m[2] - m^2
mu4 <- raw_m[4] - 4 * m * raw_m[3] +
  6 * m^2 * raw_m[2] - 3 * m^4
x <- r_sample_f(n, p, beta)
c(mean_check = sqrt(n) * (mean(x) - m) / sqrt(v),
  variance_check = sqrt(n) * (var(x) - v) / sqrt(mu4 - v^2))

For each check separately, an absolute value below

1.96

is consistent with its approximate

95%

confidence interval. Repeat with many independent samples: the two standardized statistics should each resemble a standard normal distribution, with about

5%

outside that range. Also vary

p

β

and n; a systematic discrepancy that grows with

n

suggests a bias. The two checks from the same sample are generally dependent, so passing both is not itself a joint

95%

statement, nor can a finite test prove that a generator is correct.

The second central limit theorem follows from the central limit theorem for the sample variance:

n (S_{n}^{2} - v)

differs by a quantity tending to zero in probability from

n^{- 1/2} \sum_{i} [(X_{i} - m)^{2} - v]

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / b / i / Solution by

Codex  0 2026-10-25

A numerical power series evaluation needs a stopping rule. For

r = ∣ z ∣ < 1

and integer

a \geq 0

, the absolute remainder after

N

terms satisfies the polylogarithm series remainder bound

Li_{a} (z) - k = 1 \sum N \frac{z ^{k}}{k ^{a}} \leq \frac{r ^{N + 1}}{( N + 1 ) ^{a} ( 1 - r )} .

Indeed

k^{- a} \leq (N + 1)^{- a}

for

k > N

, and the remaining powers sum as a geometric series. This also covers negative

z

. The code sums terms in vectors of at most 4096 entries, checking the bound after each block; there is no loop over individual terms.

polylog <- function(a, z, tol = 1e-10, max_terms = 1000000) {
  scalar <- function(v) {
    is.numeric(v) && !is.complex(v) && length(v) == 1L &&
      is.finite(v)
  }
  if (!scalar(a) || a < 0 || a != floor(a))
    stop("a must be a nonnegative integer")
  if (!scalar(z) || abs(z) >= 1)
    stop("z must be real with abs(z) < 1")
  if (!scalar(tol) || tol <= 0 || tol >= 1)
    stop("tol must lie strictly between 0 and 1")
  if (!scalar(max_terms) || max_terms < 1 ||
      max_terms != floor(max_terms))
    stop("max_terms must be a positive integer")
  if (z == 0) return(0)
  r <- abs(z)
  N <- 0
  ans <- 0
  repeat {
    if (N >= max_terms)
      stop("requested truncation tolerance not reached")
    B <- min(4096, max_terms - N)
    k <- seq.int(N + 1, N + B)
    ans <- ans + sum(z^k * exp(-a * log(k)))
    N <- N + B
    log_bound <- (N + 1) * log(r) - a * log(N + 1) - log1p(-r)
    if (log_bound <= log(tol)) return(ans)
  }
}

This accepts every positive integer order in the question, and also order zero if

0 \in N

is the convention adopted. The bound controls absolute truncation error; floating-point rounding is a separate source of error. Close to

∣ z ∣ = 1

, convergence can be slow. If the resource limit is reached, the function reports failure instead of silently claiming the requested accuracy. Useful exact checks are

Li_{1} (z) = - lo g (1 - z)

and

Li_{0} (z) = z / (1 - z)

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / a / iii / Solution by

Codex  0 2026-10-25

The following R function performs inverse transform sampling with vectorized numerical computation. It checks finite real scalar parameters, including a nonnegative integer sample count, and gives beta its requested default.

r_sample_f <- function(n, p, beta = 1) {
  scalar <- function(v) {
    is.numeric(v) && !is.complex(v) && length(v) == 1L &&
      is.finite(v)
  }
  if (!scalar(n) || n < 0 || n != floor(n))
    stop("n must be a nonnegative integer")
  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 (n == 0) return(numeric(0))
  v <- runif(n)
  (log1p(-p) - log(-expm1(v * log(p)))) / beta
}

Here

V = 1 - U

is again a uniform distribution draw, so the expression is the quantile function from (ii). The function expm1(t) evaluates

e^{t} - 1

accurately when

t

is small, avoiding cancellation in

1 - p^{V}

; log1p(-p) evaluates

lo g (1 - p)

. For default bounds, runif excludes both endpoints, so this implementation does not deliberately generate the infinite endpoint of the quantile function. It returns a numeric vector of length n, including an empty vector when n = 0. ^[ref] ^[ref]

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / a / ii / Solution by

Codex  0 2026-10-25

Solving

u = F (x)

gives

1 - (1 - p) e^{- β x} = p^{1 - u}

. Hence the quantile function is

F^{- 1} (u) = - \frac{1}{β} lo g (\frac{1 - p ^{1 - u}}{1 - p}), 0 < u < 1.

For such

u

p < p^{1 - u} < 1

, making the logarithm's argument lie in

(0, 1)

and the output positive. Its endpoint limits are

F^{- 1} (0) = 0

and

F^{- 1} (1) = + \infty

. Substitution into (i)'s cumulative distribution function recovers

u

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / a / i / Solution by

Codex  0 2026-10-25

Let

r = 1 - p

. Since

\frac{d}{d x} lo g (1 - r e^{- β x}) = \frac{β r e ^{- β x}}{1 - r e ^{- β x}},

integration from zero to

x > 0

gives

F (x) = (C / β) [lo g (1 - r e^{- β x}) - lo g p]

. Letting

x \to \infty

and imposing unit total probability gives

C = - \frac{β}{lo g p} > 0, F (x) = ⎩ ⎨ ⎧ 0, 1 - \frac{lo g ( 1 - ( 1 - p ) e ^{- β x} )}{lo g p}, x \leq 0, x > 0.

The probability density function is positive for

x > 0

because

0 < (1 - p) e^{- β x} < 1

. Consequently the cumulative distribution function is increasing, starts at zero and tends to one. This verifies that the derived normalizing constant defines a probability distribution. The family is the exponential-logarithmic distribution with the complementary mixing parameter

r = 1 - p

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 2 / b / ii / Solution by

Codex  0 2026-10-25

Define

a_{y} = min {1, λ^{y} /7}

. Given

Y = y

, the direct-return branch has probability

a_{y}

, while the fallback branch distributes its remaining mass uniformly over

0, \dots, y

. The law of total probability therefore gives, for integer

x \geq 0

P (X = x) = e^{- λ} \frac{λ ^{x}}{x !} min {1, \frac{λ ^{x}}{7}} + e^{- λ} y = x \sum \infty \frac{λ ^{y}}{y ! ( y + 1 )} (1 - min {1, \frac{λ ^{y}}{7}}) .

The summation includes

y = x

: a failed direct-return test can still return

x

through the discrete uniform distribution draw. There is no rejection loop and no additional normalization. Summing the displayed probability mass function over

x

gives

\sum_{y} p_{y} [a_{y} + (1 - a_{y}) (y + 1) / (y + 1)] = 1

. At

λ = 1

, for example,

a_{y} = 1/7

at every

y

; it does not become one. This is an instance of the accept-or-fallback mixture identity.

Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 2 / b / i / Solution by

Codex  0 2026-10-25

Because

U \in [0, 1]

, the conditional acceptance probability is

a_{k} = min (1, λ^{k})

, even when the printed threshold exceeds one. Repeating rejected trials gives the clipped acceptance reweighting identity:

P (X = k) = \frac{e ^{- λ} λ ^{k} a _{k} / k !}{\sum _{j \geq 0} e ^{- λ} λ ^{j} a _{j} / j !} .

For

0 < λ < 1

a_{k} = λ^{k}

, including

a_{0} = 1

. The denominator is

e^{- λ} \sum_{k} λ^{2 k} / k! = e^{λ^{2} - λ}

, so

X \sim {Poisson (λ^{2}), Poisson (λ), 0 < λ < 1, λ \geq 1.

For

λ \geq 1

every proposal is accepted, so reweighting by

λ^{k}

would be incorrect. At

λ = 1

both formulas agree. Termination has probability one since the acceptance probability is

e^{λ^{2} - λ} > 0

in the first case and one in the second.