Past exam of the simulation and statistical programming course of the University of Oxford / 2023 (A12, A13766W1)

a

i

Solution

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.

ii

Solution

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.

iii

Solution

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.

b

i

Solution

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 .

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.

ii

Solution

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

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.

iii

Solution

With a fixed budget of

n

independent proposal trials, the accepted count

K_{n}

has binomial distribution with success probability

1/ M

. Hence

E K_{n} = \frac{n}{M}, Var (K_{n}) = \frac{n ( M - 1 )}{M ^{2}} .

For

M > 1

, the central limit theorem gives the one-sided approximation

P (K_{n} \geq N) \approx Φ (\frac{n - MN}{n ( M - 1 )}) .

Use

z = 2.33

, since now only the lower tail of the accepted count matters. Requiring

n - MN \geq z n (M - 1)

and solving the resulting quadratic equation in

t = n

gives the one-sided normal approximation for a rejection proposal budget

n_{0} = (\frac{z M - 1 + z ^{2} ( M - 1 ) + 4 MN}{2})^{2} .

Equivalently, before rounding, this is

MN + z^{2} (M - 1) /2 + (z /2) (M - 1) [z^{2} (M - 1) + 4 MN]

. For the value in (a)(iii),

n_{0} = (\frac{2.33 M - 1 + 5.4289 ( M - 1 ) + 266256 M}{2})^{2} .

When

M = 1

, every proposal is accepted and

n_{0} = N = 66, 564

exactly. Otherwise the positive safety margin above

MN

allows for random fluctuations in the number accepted. The average cost in (b)(ii) describes stopping at the

N

th acceptance; this answer chooses a deterministic budget that reaches that acceptance with approximately

99%

probability.

The approximation requires a sufficiently accurate normal distribution approximation to the binomial distribution. An exact

99%

budget can instead be found from its tail probabilities. Also, a

99%

chance of enough accepted draws and approximate

99%

precision conditional on obtaining them are separate claims: their intersection is not automatically a

99%

event.

c

i

Solution

Draw independent and identically distributed random variables

Y_{1}, \dots, Y_{n} \sim q

and use the importance weight

w = p / q

. The importance sampling estimator is

θ_{n}^{IS} = \frac{1}{n} i = 1 \sum n w (Y_{i}) ϕ (Y_{i}) .

The envelope supplies the support condition for importance sampling. Since

E_{q} [w (Y) ϕ (Y)] = E_{p} ϕ (X) = θ

, this is an unbiased estimator. Moreover,

0 \leq w \leq M

, so

E_{q} [w (Y)^{2} ϕ (Y)^{2}] = \int w (x) ϕ (x)^{2} p (x) d x = E_{p} [w (X) ϕ (X)^{2}] \leq M E_{p} [ϕ (X)^{2}] = M (1 + θ^{2}) .

Here

E_{p} ϕ^{2} = Var_{p} (ϕ) + θ^{2} = 1 + θ^{2}

. Subtracting the squared expected value and dividing by

n

proves the bounded-weight importance-sampling moment bound in the requested form:

Var (θ_{n}^{IS}) = \frac{E _{p} [ w ϕ ^{2} ] - θ ^{2}}{n} \leq \frac{M + ( M - 1 ) θ ^{2}}{n} .

ii

Solution

The definition of

a

and (c)(i) give

a \leq M

. Nonnegativity of the variance strengthens the lower bound to

\frac{θ ^{2}}{1 + θ ^{2}} \leq a \leq M .

In fact the variance of one weighted observation is exactly

v_{IS} := Var_{q} (w (Y) ϕ (Y)) = a (1 + θ^{2}) - θ^{2} .

It is finite by the envelope. For

v_{IS} > 0

, the central limit theorem gives approximate two-sided coverage

2Φ (0.01 n / v_{IS}) - 1

, so the paper's rounded quantile yields

n_{IS} = ⌈ 66, 564 [a (1 + θ^{2}) - θ^{2}] ⌉ .

v_{IS} = 0

, the weighted observation equals

θ

almost surely, and one draw is sufficient. These are theoretical budgets:

a

and

θ

may themselves be unknown and need preliminary estimation in an implementation.

Compare with the average rejection sampling cost

66, 564 M

. Ignoring integer rounding, importance sampling uses fewer proposal draws at the same approximate precision precisely when

a (1 + θ^{2}) - θ^{2} < M, a < \frac{M + θ ^{2}}{1 + θ ^{2}} .

For literal integer budgets, the comparison is

⌈ 66, 564 v_{IS} ⌉ < 66, 564 M

. The simpler strict inequality is the intended leading-order condition. The bound

a \leq M

alone does not ensure an improvement when

θ \neq = 0

; this is a proposal-budget variance comparison of importance and rejection sampling.

2

a

i

Solution

Write

w = p / q

where

q > 0

, assigning zero elsewhere. The envelope ensures

0 \leq w / M \leq 1

, and

π

is nonnegative. Integrating its joint probability density over the auxiliary coordinate gives

\int_{0}^{1} π (x, u) d u = Mq (x) \frac{w ( x )}{M} = p (x) .

Consequently

\int_{Ω} \int_{0}^{1} π (x, u) d u d x = \int_{Ω} p (x) d x = 1

. This proves normalization of a probability distribution; the change of integration order is justified by Tonelli's theorem.

The coordinates are generally dependent. For

p (x) > 0

, their conditional probability density is

π (u ∣ x) = \frac{M}{w ( x )} 1_{[0, w (x) / M]} (u) .

Thus

U ∣ X = x

has a uniform distribution on an interval whose endpoint depends on

x

. They are independent exactly when

w (x)

is constant for

p

-almost every

x

, since only then do these conditional laws agree. In particular,

p = q

is an exception:

X \sim p

and

U \sim Unif (0, 1/ M)

are independent. An unconditional answer that they can never be independent would therefore be too strong.

ii

Solution

Integrating the joint probability density over the auxiliary coordinate, as in (i), gives the marginal distribution

π_{X} (x) = \int_{0}^{1} π (x, u) d u = p (x) .

In particular, augmenting with

U

preserves the target probability distribution of

X

. The other marginal distribution, though not needed to identify

X

, is

π_{U} (u) = M \int_{Ω} q (x) 1_{{w (x) \geq M u}} d x

for

0 < u < 1

. It need not be a uniform distribution on

(0, 1)

iii

Solution

For the proposed joint probability density

ρ

, define the acceptance indicator function

A_{i} := 1_{{V_{i} \leq w (Y_{i}) / M}} .

On the proposal's support the importance weight is

π (Y_{i}, V_{i}) / ρ (Y_{i}, V_{i}) = M A_{i}

. Therefore the requested augmented-space importance sampling estimator is

θ_{n}^{IS} = \frac{M}{n} i = 1 \sum n A_{i} ϕ (Y_{i}) .

Its expected value is

θ

, provided

E_{p} ∣ ϕ (X) ∣ < \infty

. This is not the ordinary estimator

n^{- 1} \sum_{i} w (Y_{i}) ϕ (Y_{i})

from question 1. Taking the conditional expectation of this augmented estimator given all

Y_{i}

produces that ordinary estimator, because

E [A_{i} ∣ Y_{i}] = w (Y_{i}) / M

. This is the Rao-Blackwell identity for a fixed-denominator rejection estimator.

iv

Solution

Use the same acceptance indicator functions and set

K_{n} = \sum_{i = 1}^{n} A_{i}

S_{n} = \sum_{i = 1}^{n} A_{i} ϕ (Y_{i})

. The rejection sampling estimator based on these

n

proposals is the mean over the accepted observations:

θ_{n}^{RS} = \frac{S _{n}}{K _{n}} (K_{n} > 0) .

For

K_{n} = 0

, any fixed convention, such as zero, is harmless asymptotically. The

A_{i}

are independent and identically distributed random variables with Bernoulli distribution and success probability

1/ M

, so the strong law of large numbers gives

K_{n} / n \to 1/ M

almost surely. In particular

K_{n} > 0

eventually almost surely. The exact identity

θ_{n}^{IS} = \frac{M K _{n}}{n} θ_{n}^{RS}

then holds whenever

K_{n} > 0

. Wherever

S_{n} \neq = 0

, it follows that

\frac{θ _{n}^{IS}}{θ _{n}^{RS}} = \frac{M K _{n}}{n} ⟶ 1 almost surely .

A qualification is missing from the printed claim: the quotient must be defined eventually. A convenient sufficient assumption is

E_{p} ∣ ϕ ∣ < \infty

and

θ \neq = 0

. Indeed, another application of the strong law of large numbers gives

S_{n} / n \to θ / M

, ensuring that

S_{n}

is eventually nonzero. This proves the stated quotient limit under that assumption.

As printed, the assertion is false as a universally defined quotient. For example,

ϕ \equiv 0

makes both estimates zero for every

n

. Even excluding that example is insufficient: take

p = q

to be the uniform distribution on

(0, 1)

M = 1

, and

ϕ (y) = 1

for

y < 1/2

and

- 1

otherwise. Then

S_{n}

is a simple symmetric random walk, which returns to zero infinitely often almost surely, leaving the quotient undefined infinitely often.

With integrability but no assumption

θ \neq = 0

, the exact multiplicative identity remains valid, and both estimates converge to

θ

almost surely. In particular

θ_{n}^{IS} - θ_{n}^{RS} \to 0

almost surely. One may alternatively define the quotient's removable

0/0

cases as

M K_{n} / n

; that explicitly extended quantity always converges to one. These are valid repairs of the source's unrestricted statement.

b

i

Solution

Fix a proposal

y

with

q (y) > 0

, and abbreviate

B = M q (y)

H = h (y)

and

P = p (y)

. The envelope says

0 \leq H \leq P \leq B

. Assuming fresh independent uniform distribution draws at both stages, the law of total probability gives

P (accept ∣ Y = y) = \frac{H}{B} + (1 - \frac{H}{B}) \frac{P - H}{B - H} = \frac{H}{B} + \frac{P - H}{B} = \frac{p ( y )}{M q ( y )} .

The second term includes both reaching the residual test and passing it. If

H = B

, the envelope forces

P = B

and the first test accepts with probability one; the second-stage

0/0

ratio is never evaluated. This is squeeze rejection sampling.

ii

Solution

Average the conditional probability in (i) over

Y \sim q = q / Z_{q}

P (accept) = \int_{Ω} \frac{q ( y )}{Z _{q}} \frac{p ( y )}{M q ( y )} d y = \frac{1}{M Z _{q}} \int_{Ω} p (y) d y = \frac{Z _{p}}{M Z _{q}} .

This is exactly the acceptance mass for an unnormalized rejection envelope. In normalized form the envelope is

p \leq (M Z_{q} / Z_{p}) q

, giving the same acceptance probability. Here

M

bounds the ratio of the unnormalized functions, so its value need not equal the envelope constant for the normalized probability density functions.

The squeeze changes the number of expensive target evaluations, while preserving the per-proposal acceptance probability. The expected number of proposals per accepted draw is

M Z_{q} / Z_{p}

iii

Solution

For a measurable set

C

, (i) gives the joint acceptance probability

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

Divide by

Z_{p} / (M Z_{q})

from (ii). The resulting conditional probability is

P (Y \in C ∣ accept) = \frac{1}{Z _{p}} \int_{C} p (y) d y = \int_{C} p (y) d y .

Thus every accepted value has target probability density function

p

. Repeating rejected trials with fresh independent random variables does not change this law: the contribution from acceptance on trial

j

is multiplied by

(1 - Z_{p} / (M Z_{q}))^{j - 1}

, and summing that geometric series gives the same conditional law. Independent repetitions produce independent target draws.

iv

Solution

The proposal's normalization of a probability distribution is

Z_{q} = \int_{R} e^{- ∣ x ∣} d x = 2, q (x) = \frac{1}{2} e^{- ∣ x ∣} .

This is the Laplace distribution with location zero and scale one. For example, generate its magnitude from the exponential distribution with rate one and give it an independent fair sign. The target is the standard normal distribution, with

Z_{p} = 2 π

To obtain the tight envelope, let

t = ∣ x ∣

. By completing the square,

\frac{p ( x )}{q ( x )} = e^{- t^{2} /2 + t} = e^{1/2} e^{- (t - 1)^{2} /2} \leq e^{1/2} .

Equality occurs at

x = \pm 1

, so

M_{*} = e

is the smallest valid envelope constant. It maximizes the acceptance probability for this fixed proposal. The supplied lower bound must be clipped at zero to satisfy the algorithm's nonnegativity requirement:

h (x) = max {0, 1 - x^{2} /2} .

This is valid because

e^{- s} \geq 1 - s

for

s \geq 0

, while

p

is everywhere positive. Its support is

[- 2, 2]

Here is the fully specified squeeze rejection sampling procedure. Draw

Y

with probability density function

q (y) = e^{- ∣ y ∣} /2

, independently draw

U \sim Unif (0, 1)

, and compute

H = max (0, 1 - Y^{2} /2)

and

B = e e^{- ∣ Y ∣}

. If

U \leq H / B

, return

Y

immediately. Otherwise evaluate

P = e^{- Y^{2} /2}

, independently draw

V \sim Unif (0, 1)

, and return

Y

V \leq (P - H) / (B - H)

. On rejection restart with fresh draws. In this example

B > H

everywhere, so the residual denominator is positive.

Avoiding a target evaluation means accepting at the first test. Its unconditional probability for one proposal is

P (no target evaluation) = \int_{R} q (y) \frac{h ( y )}{M q ( y )} d y = \frac{1}{2 e} \int_{- 2}^{2} (1 - y^{2} /2) d y = \frac{2 2}{3 e} \approx 0.571843.

For comparison, the overall acceptance probability is

2 π / (2 e) = π / (2 e) \approx 0.760173

. The first-stage saving

0.571843

is the fraction of proposals requiring no target evaluation, including both eventual accepted and rejected proposals in the denominator; it is not conditional on acceptance.

3

a

Solution

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 .

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

b

Solution

An inverse matrix exists only if the lower triangular matrix has no zero diagonal entries; this necessary assumption is implicit in the request for

L^{- 1}

. Partition

L = (a b 0 C) .

Writing its inverse in the same block form and multiplying gives the recursive inverse of a lower triangular matrix:

L^{- 1} = (a^{- 1} - C^{- 1} b / a 0 C^{- 1}) .

Indeed the lower-left block of the product is

b / a - C (C^{- 1} b / a) = 0

, and the diagonal blocks are identities. This gives an implementation with no call to solve:

ltinv <- function(L) {
  n <- nrow(L)
  if (n == 1L) return(matrix(1 / L[1L, 1L], 1L, 1L))
  j <- 2L:n
  Cinv <- ltinv(L[j, j, drop = FALSE])
  B <- matrix(0, n, n)
  B[1L, 1L] <- 1 / L[1L, 1L]
  B[j, j] <- Cinv
  B[j, 1L] <- -as.vector(Cinv %*% L[j, 1L, drop = FALSE]) / L[1L, 1L]
  B
}

The preserved matrix dimensions also handle

n = 2

. A mathematical induction using the block identity proves that every returned matrix is the desired inverse. Applied to the Cholesky decomposition from (a), all diagonal entries are positive, so all divisions are valid.

c

i

Solution

Each state has exactly

k (n - k)

neighbours: choose one occupied coordinate and one unoccupied coordinate to swap. The normalized proposal distribution is therefore

q (v ∣ u) = \frac{1}{k ( n - k )} 1_{{v \sim u}} .

The swap relation is symmetric and the number of neighbours is constant, so

q (v ∣ u) = q (u ∣ v)

. The discrete uniform distribution used as the prior distribution is also constant on

S_{k}

, with mass

(k n)^{- 1}

. Both factors cancel in the Metropolis–Hastings acceptance probability.

Let

Q = Σ^{- 1}

and

d = v - u

. The displayed multivariate normal distribution presupposes that the covariance matrix

Σ

is a positive-definite matrix, so

Q

is a symmetric matrix. Expanding the difference of quadratic forms gives

lo g \frac{π ( v ∣ x )}{π ( u ∣ x )} = - \frac{1}{2} {(x - v)^{T} Q (x - v) - (x - u)^{T} Q (x - u)} = d^{T} Q (x - u) - \frac{1}{2} d^{T} Q d = d^{T} Q (x - \frac{u + v}{2}) .

Consequently

lo g α (u, v) = min {0, (v - u)^{T} Σ^{- 1} (x - \frac{u + v}{2})} .

The truncation at zero is essential: a log acceptance ratio can be positive, but the logarithm of an acceptance probability cannot.

For the particular proposed swap, let

a = i_{1}

be the coordinate changing from zero to one and

b = i_{2}

the coordinate changing from one to zero. Then

d = e_{a} - e_{b}

, and the untruncated scalar simplifies further to

(Q (x - u))_{a} - (Q (x - u))_{b} - \frac{1}{2} (Q_{aa} + Q_{bb} - 2 Q_{ab}) .

Only two rows and four entries of

Q

are needed for this evaluation. This is Metropolis sampling on fixed-cardinality subsets.

ii

Solution

The source's function name suggests a log return value, and (c)(i) asks for

lo g α

, while the wording here and in (c)(iv) asks for the acceptance probability. We use the literal probability-returning convention and explicitly exponentiate the result of (c)(i):

mylogapp <- function(x, u, v, sigma_inv) {
  d <- v - u
  log_ratio <- as.numeric(crossprod(d, sigma_inv %*% (x - (u + v) / 2)))
  log_alpha <- min(0, log_ratio)
  exp(log_alpha)
}

Thus mylogapp returns

α \in [0, 1]

, and (c)(iv) compares it directly with a uniform distribution draw. Taking the minimum before exponentiating prevents overflow. The arguments

u, v

are assumed to be valid neighbouring states. A log-returning convention would instead return log_alpha and compare log(runif(1)) with it; the two conventions implement the same Metropolis–Hastings algorithm, provided the caller uses the matching comparison.

iii

Solution

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.

iv

Solution

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