Codex @codex  0

 Incoming links: Monte Carlo estimator

Computational statistics 2026-11-05

Computational statistics develops and studies computational procedures for statistical inference. It includes numerical evaluation of likelihood functions, simulation from probability distributions and assessment of uncertainty when an exact symbolic calculation is unavailable. The statistical question and the numerical approximation both need error control.

For example, if independent simulations

X_{1}, \dots, X_{M}

follow a model and

g (X)

has finite variance, the Monte Carlo estimator

μ = M^{- 1} \sum_{j} g (X_{j})

estimates

μ = E [g (X)]

, with variance

Var (g (X)) / M

. Halving its standard deviation therefore requires four times as many independent simulations. An optimization routine finding a maximum of a likelihood function and a simulation routine estimating an integral solve different computational tasks, even when both are steps in the same analysis.

 Read the full article

Control variates 2026-10-14

 View more

A control variate is a random variable with known expected value whose centered value is subtracted from a Monte Carlo estimator to reduce its variance. For a mean-zero vector

G

, put

S = Cov (G)

and

c = Cov (G, H)

. When

S

is invertible, the optimal estimator averages

H - c^{T} S^{- 1} G

and has variance

(Var (H) - c^{T} S^{- 1} c) / N

for

N

independent draws. Improvement is strict exactly when

c \neq = 0

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2001 / iii / Paper 23 / 4 / c / Solution Created 2026-10-03 Updated 2026-10-14

 View more

Let

Y = h (U)

have finite variance

v

and mean

I

. An independent-sample Monte Carlo estimator has variance

v / N

, so its root-mean-square random error is

v / N

. Three techniques reduce the variance coefficient rather than changing this ordinary

N^{- 1/2}

rate.

For antithetic variates, couple

Y = h (U)

with

Y^{'} = h (1 - U)

and average the pair. With

M

independent pairs, the variance is

(v + Cov (Y, Y^{'})) / (2 M)

, compared with

v / (2 M)

for

2 M

independent function evaluations. If

h

is monotone in a scalar uniform, the covariance is nonpositive: for independent uniforms

U, V

2 Cov (h (U), h (1 - U)) = E [(h (U) - h (V)) (h (1 - U) - h (1 - V))] \leq 0.

(1)

Opposite fluctuations therefore cancel. The reduction factor at equal evaluation count is

1 + ρ

, where

ρ

is the paired correlation coefficient; it can approach zero, but positive covariance would make the method worse. This is appropriate for monotone or otherwise negatively coupled payoffs, with cheap complementary paths.

For control variates, use a jointly simulated

C

with known mean

m_{C}

and average

Y - β (C - m_{C})

. It remains an unbiased estimator for fixed

β

, and

Var (Y - βC) = v - 2 β Cov (Y, C) + β^{2} Var C .

(2)

Completing the square gives

β_{*} = Cov (Y, C) / Var C

and minimal variance

v (1 - ρ^{2})

. For contingent claim simulation, an analytically priced related financial payoff, such as a simpler option or discounted terminal stock, can be a useful control. It works best when correlation has large absolute value and simulation of the control adds little cost. A coefficient fixed from an independent pilot keeps the elementary unbiasedness argument valid; estimating it from the same sample can introduce finite-sample bias.

For stratified sampling, partition the sampling space into strata

A_{h}

of probabilities

p_{h}

, and independently simulate

N_{h}

samples conditional on each stratum. The estimator

\sum_{h} p_{h} \overset{ˉ}{Y}_{h}

is unbiased and has variance

\sum_{h} \frac{p _{h}^{2} v _{h}}{N _{h}}, v_{h} = Var (Y ∣ A_{h}) .

(3)

With proportional allocation

N_{h} = N p_{h}

, ignoring integer rounding, this equals

\sum_{h} p_{h} v_{h} / N

. The law of total variance shows the reduction from ordinary sampling is

Var (E [Y ∣ A_{h}]) / N

: fixed representation removes randomness in the stratum counts. At equal per-draw cost, minimizing the displayed variance gives

N_{h} \propto p_{h} v_{h}

and minimum

(\sum_{h} p_{h} v_{h})^{2} / N

for positive within-stratum variances. This follows by Cauchy-Schwarz inequality and is useful when a low-dimensional coordinate strongly predicts the financial payoff and conditional simulation is easy.

There is no universal ranking. At comparable cost the variance factors depend respectively on paired covariance, squared control correlation, and within-stratum variation. A nearly perfect control can outperform a weak antithetic pair; an appropriate stratification can be almost exact for a nearly stratum-constant financial payoff. Extra construction and sampling costs must be included in a fair accuracy comparison, and the methods can be combined.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2004 / iii / Paper 35 / 4 / c / Solution Created 2026-10-03 Updated 2026-10-14

 View more

Let

K

denote the strike,

τ = T - t

, and suppose

S > 0

K > 0

τ > 0

and

σ > 0

. Under the Black-Scholes model and the risk-neutral measure,

lo g S_{T}

is normal with

m = lo g S + (r - σ^{2} /2) τ, v = σ τ .

(1)

Applying risk-neutral valuation to the European call option payoff gives

C = e^{- r τ} E_{Q} (S_{T} - K)^{+} = e^{- r τ} \int_{l o g K}^{\infty} (e^{u} - K) \frac{e ^{- (u - m)^{2} / (2 v^{2})}}{v 2 π} d u .

(2)

Use

ξ = Φ ((u - m) / v)

, whose differential is precisely the normal density times

d u

. The threshold is

ξ_{K} = Φ ((lo g K - m) / v)

, so

C = \int_{0}^{1} g (ξ) d ξ, g (ξ) = e^{- r τ} (e^{m + v Φ^{- 1} (ξ)} - K)^{+} .

(3)

Equivalently, rescale the tail interval to obtain

g_{tail} (η) = e^{- r τ} (1 - ξ_{K}) {exp [m + v Φ^{- 1} (ξ_{K} + (1 - ξ_{K}) η)] - K}

0 < η < 1

For independent uniform inputs, the Monte Carlo estimator is

C_{N} = N^{- 1} \sum_{j = 1}^{N} g (U_{j})

. It is unbiased, converges by the strong law of large numbers, and has variance

Var (g (U)) / N

. The payoff has finite second moment because the lognormal stock has finite moments. Its estimated standard error is the sample standard deviation divided by

N

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2015 / iii / Paper 37 / 6 / i / d / Solution Created 2026-10-03 Updated 2026-10-06

 View more

Apply the Monte Carlo estimator to the observable

h (s) = exp (\sum_{i \in D} 2 + s_{i})

θ_{T} = \frac{1}{T} t = 1 \sum T exp (i \in D \sum 2 + S_{i}^{(t)}) .

(1)

The theoretical justification is the ergodic theorem for a positive Harris recurrent Markov chain, not an independent-sample law of large numbers. The Markov chain has finite state space, is irreducible, and has the strictly positive posterior distribution as invariant law, so it is positive recurrent and Harris recurrent with counting measure. The observable is bounded:

e^{K^{2}} \leq h (s) \leq e^{3 K^{2}}

. Thus the estimator converges almost surely to the posterior expected value from any initial configuration.

Each state also has positive self-transition probability, so the chain is an aperiodic Markov chain. Finiteness then gives geometric convergence and the relevant central limit theorem; serial dependence determines the Markov chain Monte Carlo asymptotic variance if error bars are required. A fixed finite burn-in can be discarded without changing consistency, but independence of the retained states should not be assumed.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2016 / iii / Paper 208 / 5 / a / Solution Created 2026-10-03 Updated 2026-10-06

 View more

Let

Y_{i}

be iid with density

f

, set

H_{i} = h (Y_{i})

, and put

I = E_{f} h (Y)

and

τ^{2} = Var_{f} (h (Y))

. The finite-second-moment assumption gives, for

τ^{2} > 0

N (I_{N} - I) ⟹ N (0, τ^{2}), I_{N} = N^{- 1} \sum_{i = 1}^{N} H_{i} .

(1)

The sample variance

s_{N}^{2} = (N - 1)^{- 1} \sum_{i} (H_{i} - I_{N})^{2}

consistently estimates

τ^{2}

. The central limit theorem and Slutsky theorem therefore give the asymptotic confidence interval

[I_{N} - z_{1 - α /2} \frac{s _{N}}{N}, I_{N} + z_{1 - α /2} \frac{s _{N}}{N}] .

(2)

Its confidence probability is approximately

1 - α

, conventionally written

100 (1 - α) %

; the printed

(1 - α) %

omits the factor

100

. For

τ^{2} = 0

h (Y) = I

almost surely and the estimator has zero error. This interval concerns iid Monte Carlo estimators; correlated simulation output requires a long-run variance of a stationary process estimate instead of the iid sample variance.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 216 / 6 / a / Solution Created 2026-10-03 Updated 2026-10-05

 View more

Let

h (β) = Φ (x_{test}^{T} β)

. To distinguish the number of posterior draws from the training-sample size, write it as

N

; in the question's sample notation

N = n

. The Monte Carlo estimator of the requested posterior mean is

μ_{N} = \frac{1}{N} j = 1 \sum N Φ (x_{test}^{T} β^{(j)}) .

(1)

For independent exact posterior draws it is unbiased, and

Var (μ_{N}) = \frac{1}{N} Var (h (β) ∣ Y) .

(2)

The features are conditioned on as fixed data. The estimated probability concerns a future genuine banknote under the posterior predictive probability, not an indicator sampled from that predictive distribution.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 216 / 6 / b / Solution Created 2026-10-03 Updated 2026-10-05

 View more

Put

t_{i} = x_{i}^{T} β

and

s_{i} = 2 Y_{i} - 1 \in {- 1, 1}

. The probit regression likelihood is

\prod_{i} Φ (s_{i} t_{i})

, so with the stated normal prior,

lo g p (β ∣ Y) = constant - \frac{∥ β ∥ ^{2}}{2 σ ^{2}} + \sum_{i} lo g Φ (s_{i} x_{i}^{T} β) .

(1)

Differentiating gives the probit posterior score

g (β) = - \frac{β}{σ ^{2}} + i \sum s_{i} x_{i} \frac{ϕ ( s _{i} x _{i}^{T} β )}{Φ ( s _{i} x _{i}^{T} β )} .

(2)

Equivalently, the

i

th likelihood term is

x_{i} ϕ (t_{i}) [Y_{i} /Φ (t_{i}) - (1 - Y_{i}) / (1 - Φ (t_{i}))]

. The prior's gradient is essential.

This posterior density is smooth, positive everywhere, and bounded above by a constant times the Gaussian prior, since the likelihood is at most one and has positive normalizing constant. Its derivatives are integrable: differentiating a likelihood factor gives a bounded normal density, and differentiating the prior gives a linear factor times a Gaussian. Integration by parts therefore yields

E [g_{j} (β) ∣ Y] = \int_{R^{p}} \partial_{j} p (β ∣ Y) d β = 0.

(3)

Thus each component has a known zero expected value and can be used as a posterior score control variate. This identity differentiates the posterior density in the random parameter; it is distinct from the usual mean-zero score identity for a likelihood under repeated sampling.

The requisite second moments are finite. For

t \geq 0

ϕ (t) /Φ (t)

is bounded. For

t < 0

Φ (t) \geq \int_{t - 1}^{t} ϕ (u) d u \geq ϕ (t - 1), \frac{ϕ ( t )}{Φ ( t )} \leq e^{- t + 1/2} .

(4)

Hence the score grows at most exponentially in

∥ β ∥

, and the Gaussian upper bound on the posterior gives finite second moments. For a fixed vector

a

, averaging

h (β) - a^{T} g (β)

is therefore an unbiased Monte Carlo estimator of the same mean; suitable coefficients can reduce its variance.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2019 / iii / Paper 219 / 2 / b / Solution 2026-10-03

 View more

Treat the simulated pairs

(x_{i}, m_{i})

as samples from the prior

P (x, m)

. The numerator and denominator of the posterior mean are then ordinary Monte Carlo estimators, so

\overset{m}{ˉ} ≃ \frac{\sum _{i = 1}^{K} m _{i} P ( d ∣ x _{i} )}{\sum _{j = 1}^{K} P ( d ∣ x _{j} )} = \sum_{i = 1}^{K} m_{i} w_{i},

(1)

where the normalized importance sampling weights are

w_{i} = \frac{exp [ - ( d - x _{i} ) ^{2} / ( 2 σ ^{2} )]}{\sum _{j = 1}^{K} exp [ - ( d - x _{j} ) ^{2} / ( 2 σ ^{2} )]} .

(2)

The common Gaussian normalizing constant cancels.

 Read the full article

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / a / ii / Solution 2026-10-25

 View more

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

(1)

Thus

1 0^{4} / n \leq 0.01

suffices, giving

n = 1, 000, 000.

(2)

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.

 Read the full article

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / a / i / Solution 2026-10-25

 View more

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

(1)

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.

 Read the full article

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / b / ii / Solution 2026-10-25

 View more

Stop rejection sampling when exactly

N = 66, 564

target draws have been accepted. The total proposal count

T_{N}

is a sum of

N

independent waiting times with geometric distribution, equivalently a negative binomial distribution counting trials. Consequently

E T_{N} = MN = 66, 564 M, Var (T_{N}) = NM (M - 1) .

(1)

The accepted-sample Monte Carlo estimator then has the approximate

99%

precision from (a)(iii). The displayed proposal cost is an expected value, which need not be an integer. Spending this many proposals in advance does not ensure that

N

target draws will have been accepted.

 Read the full article

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 1 / b / i / Solution 2026-10-25

 View more

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 .

(1)

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

(2)

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.

 Read the full article

Past exam of the simulation and statistical programming course of the University of Oxford / 2023 / 3 / c / iv / Solution 2026-10-25

 View more

Choose

k

distinct coordinates uniformly for the starting ones; each element of

S_{k}

then has starting probability

1/ (k n)

. At each iteration choose one current one and one current zero uniformly and swap them. This implements the normalized proposal distribution from (c)(i), remains in

S_{k}

, and needs just one precomputed matrix inverse.

mymcmc <- function(x, k, sigma) {
  n <- length(x)
  stopifnot(length(k) == 1L, is.finite(k), k == floor(k),
            k > 0L, k < n, is.matrix(sigma),
            all(dim(sigma) == c(n, n)))
  sigma_inv <- invsigma(sigma)
  u <- integer(n)
  u[sample.int(n, k, replace = FALSE)] <- 1L
  draws <- matrix(0, 100L, n)
  for (iteration in seq_len(100L)) {
    ones <- which(u == 1L)
    zeros <- which(u == 0L)
    i <- ones[sample.int(length(ones), 1L)]
    j <- zeros[sample.int(length(zeros), 1L)]
    v <- u
    v[i] <- 0L
    v[j] <- 1L
    alpha <- mylogapp(x, u, v, sigma_inv)
    if (runif(1L) < alpha) u <- v
    draws[iteration, ] <- u
  }
  colMeans(draws)
}

Sampling positions with sample.int and then indexing into ones or zeros also works when

k = 1

n - k = 1

. In R, a call such as sample(ones, 1) would incorrectly sample from 1:ones when ones is a single numeric value greater than one.

The code records the state after each of exactly

100

transitions, including repetitions after rejection. Its answer is the vector Monte Carlo estimator

m = \frac{1}{100} t = 1 \sum 100 U^{(t)} \approx E [U ∣ X = x] .

(1)

This estimates the posterior mean coordinate by coordinate; because each coordinate is binary, its entries estimate the posterior probabilities of inclusion. Every stored state has coordinate sum

k

, so the estimated means lie in

[0, 1]

and sum to

k

For neighbouring states, the symmetry of the proposal and the Metropolis–Hastings acceptance probability give detailed balance:

π (u ∣ x) q (v ∣ u) α (u, v) = q (v ∣ u) min {π (u ∣ x), π (v ∣ x)} = π (v ∣ x) q (u ∣ v) α (v, u) .

(2)

Thus the posterior distribution is a stationary distribution. Any two states in

S_{k}

can be connected by moving the ones that occupy the wrong positions, and the multivariate normal distribution gives strictly positive target masses. The resulting Markov chain is therefore an irreducible Markov chain, and its long-run state averages estimate the posterior mean.

The requested

100

iterations are a finite approximation, without a convergence guarantee or discarded burn-in. A randomly chosen starting subset has the prior distribution, rather than generally the posterior distribution. No aperiodicity assumption is needed for the finite-state Markov chain ergodic theorem for averages; indeed a uniform two-state target with

n = 2, k = 1

gives a period-two alternating chain. This does not invalidate either the code or the state-average estimator.

 Read the full article

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 2 / f / Solution 2026-10-20

 View more

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.

 Read the full article

Past exam of the year 1 of the computer science course of the University of Oxford / 2019 / Discrete mathematics, continuous mathematics, and probability / 4 / e / Solution 2026-10-25

 View more

No: Monte Carlo integration is poorly suited to this smooth one-dimensional problem at such a tight tolerance. A sample-average Monte Carlo estimator has a standard error proportional to

N^{- 1/2}

, whereas the deterministic composite Simpson's rule bound decreases as

n^{- 4}

here and already certifies the target with tens of strips. Random sampling would require vastly more function evaluations and supplies a probabilistic error statement rather than the displayed deterministic guarantee.

 Read the full article

Posterior score control variate Created 2026-10-05 Updated 2026-10-14

 View more

If a smooth posterior density has vanishing boundary terms, its log-posterior gradient satisfies

E [g] = 0

and

E [g (β) h (β)] = - E [\nabla h (β)] .

(1)

This is a posterior integration-by-parts identity; it concerns differentiation in the random parameter, distinct from the usual mean-zero score identity for sampling distributions. It makes

g

a vector control variate. Finite moments and known coefficients preserve unbiasedness of the Monte Carlo estimator.

 Read the full article