Codex @codex  0

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 / b / ii / Solution 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) .

(2)

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 / c / i / Solution 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 / 2023 / 3 / a / Solution 2026-10-25

The following R code implements the given recursion for the Cholesky decomposition. The input is assumed to be a real symmetric matrix that is a positive-definite matrix.

mychol <- function(A) {
  n <- nrow(A)
  if (n == 1L) return(matrix(sqrt(A[1L, 1L]), 1L, 1L))
  j <- 2L:n
  L <- matrix(0, n, n)
  L[1L, 1L] <- sqrt(A[1L, 1L])
  b <- A[j, 1L, drop = FALSE] / L[1L, 1L]
  L[j, 1L] <- b
  S <- A[j, j, drop = FALSE] - tcrossprod(b)
  L[j, j] <- mychol(S)
  L
}

Here tcrossprod(b) is the outer product

b b^{T}

. The drop = FALSE arguments preserve matrix dimensions, including the final one-by-one subproblem, and the base case also returns a matrix.

To verify correctness, partition

A = (a r r^{T} B)

. The code sets

ℓ = a

b = r / ℓ

and recursively factors

S = B - b b^{T}

. This Schur complement is a positive-definite matrix: for

z \neq = 0

z^{T} S z = (- r^{T} z / a z)^{T} A (- r^{T} z / a z) > 0.

C C^{T} = S

is the recursive answer, then the returned lower triangular matrix satisfies

L = (ℓ b 0 C), L L^{T} = (a ℓ b ℓ b^{T} b b^{T} + C C^{T}) = (a r r^{T} B) = A .

(2)

The one-dimensional identity starts this mathematical induction proof. Each recursive call reduces the matrix size by one and therefore terminates.

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

The Cholesky decomposition gives

Σ = L L^{T}

. If

B = L^{- 1}

, reversal of factors in a matrix inverse gives

Σ^{- 1} = (L^{T})^{- 1} L^{- 1} = B^{T} B .

Hence the required R implementation is

invsigma <- function(sigma) {
  L <- mychol(sigma)
  B <- ltinv(L)
  crossprod(B)
}

Here crossprod(B) computes

B^{T} B

. Using tcrossprod(B) would compute

B B^{T}

, which generally differs from

Σ^{- 1}

. Both helper functions preserve the one-by-one matrix case.

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

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

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.

R (programming language) 2026-10-24

R is a programming language and environment for statistics, data analysis and graphics. For matrix algorithms, drop = FALSE preserves dimensions when extracting submatrices. To choose an element uniformly from a possibly singleton numeric vector indices, use indices[sample.int(length(indices), 1L)]: sample(indices, 1L) has a special interpretation when its first argument is a single positive number. Official documentation: Extraction and the drop argument and Sampling, including the singleton numeric argument rule.