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 . 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 . The code sets , and recursively factors . This Schur complement is a positive-definite matrix: for ,If is the recursive answer, then the returned lower triangular matrix satisfiesThe one-dimensional identity starts this mathematical induction proof. Each recursive call reduces the matrix size by one and therefore terminates.
An inverse matrix exists only if the lower triangular matrix has no zero diagonal entries; this necessary assumption is implicit in the request for . PartitionWriting its inverse in the same block form and multiplying gives the recursive inverse of a lower triangular matrix:Indeed the lower-left block of the product is , 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 . 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.
Each state has exactly neighbours: choose one occupied coordinate and one unoccupied coordinate to swap. The normalized proposal distribution is thereforeThe swap relation is symmetric and the number of neighbours is constant, so . The discrete uniform distribution used as the prior distribution is also constant on , with mass . Both factors cancel in the Metropolis–Hastings acceptance probability.
Let and . The displayed multivariate normal distribution presupposes that the covariance matrix is a positive-definite matrix, so is a symmetric matrix. Expanding the difference of quadratic forms givesConsequentlyThe 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 be the coordinate changing from zero to one and the coordinate changing from one to zero. Then , and the untruncated scalar simplifies further toOnly two rows and four entries of are needed for this evaluation. This is Metropolis sampling on fixed-cardinality subsets.
The source's function name suggests a log return value, and (c)(i) asks for , 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 , and (c)(iv) compares it directly with a uniform distribution draw. Taking the minimum before exponentiating prevents overflow. The arguments 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.The Cholesky decomposition gives . If , reversal of factors in a matrix inverse givesHence the required R implementation is
invsigma <- function(sigma) {
L <- mychol(sigma)
B <- ltinv(L)
crossprod(B)
}Here
crossprod(B) computes . Using tcrossprod(B) would compute , which generally differs from . Both helper functions preserve the one-by-one matrix case.Choose distinct coordinates uniformly for the starting ones; each element of then has starting probability . 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 , 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 or . 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 transitions, including repetitions after rejection. Its answer is the vector Monte Carlo estimatorThis 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 , so the estimated means lie in and sum to .
For neighbouring states, the symmetry of the proposal and the Metropolis–Hastings acceptance probability give detailed balance:Thus the posterior distribution is a stationary distribution. Any two states in 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 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 gives a period-two alternating chain. This does not invalidate either the code or the state-average estimator.
Articles by others on the same topic
There are currently no matching articles.