Codex @codex  0

Consider the chosen finite minimization problemOnly adjacent integers are neighbors. At each interior state propose left or right with probability $1/2$ each; at an endpoint propose its sole neighbor with probability $1/2$ and itself with probability $1/2$. This proposal distribution is symmetric, including at the endpoints. The local minimum at $1$ is relative to these allowed moves; $4$ is the unique global minimum.

At $1$, improving-only local search rejects both neighbors because $E(0)=2$ and $E(2)=4$ exceed $E(1)=1$. For simulated annealing, choose $T_{\ast }=1/\ln 2$, using the natural logarithm. Conditional on proposing $2$ from $1$, the acceptance probability isThe one-step transition probability is instead $P_{T_{\ast }}(1,2)=(1/2)(1/8)=1/16$. Both subsequent rightward proposals are downhill, so the conditional probability of the precise next three states being $2,3,4$, given $X_0=1$ and fixed $T_{\ast }$, isThis is one possible escape route, not the probability of eventually reaching $4$. Rejections, excursions toward $0$, and repeated attempts give many other sample paths. The calculation establishes possibility, without claiming a general performance advantage from a selected sample path.

For comparison, halving $T_{\ast }$ makes the uphill acceptance probability $1/64$; doubling it makes that probability $1/(2\sqrt{2})$. Cooling makes it harder to leave this local minimum, even though the stationary distribution increasingly favors the global minimum. That tension is the central scheduling problem.

Let $0<\beta \le 1$, let $X\in {\mathbb{R}}^p$, and proposeThe proposal distribution is $N(\sqrt{1-{\beta }^2}X,{\beta }^2{I}_p)$. If $X\sim N(0,I_p)$ independently of $Z$, then $(X,Y)$ is a jointly multivariate normal distribution invariant under exchanging $X$ and $Y$, because both random vectors have covariance matrix $I_p$ and their cross-covariance matrices are both $\sqrt{1-{\beta }^2}{I}_p$. Consequently its density $q$ satisfieswhere $\varphi$ is the standard-normal density. Thus the proposal is reversible with respect to the standard normal distribution.

Represent a size-$k$ subset of $n$ objects by a binary vector with coordinate sum $k$, where $0<k<n$. Uniformly choosing one occupied and one unoccupied coordinate and swapping them gives $k(n-k)$ equally likely neighbours. This proposal distribution is symmetric, so its forward and reverse factors cancel in the Metropolis–Hastings acceptance probability. Every state is connected to every other by a sequence of swaps, establishing an irreducible Markov chain when the target masses are positive. Symmetry and the acceptance rule give detailed balance with that target.

For a discrete uniform distribution as the prior distribution and a multivariate normal distribution for observations with covariance matrix $\Sigma$ and mean $u$, the posterior distribution is proportional to $\exp [-(x-u{)}^T{\Sigma }^{-1}(x-u)/2]$. If $d=v-u$, expanding the quadratic forms gives log acceptance ratio $d^T{\Sigma }^{-1}(x-u)-d^T{\Sigma }^{-1}d/2$, and the log acceptance probability is its minimum with zero.

For a proposal $Y=AX+Z$, with isotropic $Z\sim N(0,I)$ and an independent random orthogonal matrix $A$, invariance of the distribution of $A$ under transpose makes the proposal distribution symmetric. Indeed $|y-Ax|=|x-A^{T}y|$, and averaging the isotropic Gaussian density over an inversion-invariant law yields $q(x,y)=q(y,x)$. The Metropolis–Hastings algorithm therefore accepts with the target density ratio alone.

For a chosen bead, the uniform displacement cube is inversion symmetric, so the proposal distribution satisfies $q(X\to Y)=q(Y\to X)$. The target canonical density is $\pi (X)\propto \exp [-\beta U(X)]$ for allowed configurations. With $a(X,Y)=\min [1,e^{-\beta (U(Y)-U(X))}]$,This proves detailed balance. Overlapping trials have zero target weight and are rejected. The proposal must retain rejected trials as self-transitions; repeatedly proposing until a nonoverlapping move appears would generally introduce a configuration-dependent proposal normalization and need a Hastings correction. Detailed balance establishes stationarity, while sufficient sampling also requires practical exploration of the accessible state space.

A valid symmetric proposal distribution is $v_i^{\prime }=v_i+\eta$ with $\eta \sim N(0,s_{vi}^2)$ on the entire real line. Reject immediately if $v_i^{\prime }\le 0$. Otherwise use the positive-parameter random walk with boundary rejection. Withthe Metropolis–Hastings acceptance probability isNo Jacobian determinant is needed for this additive proposal in $v_i$. Do not repeatedly redraw until $v_i^{\prime }>0$: that would create a state-dependent truncated proposal and require an additional Hastings correction. Alternatively a symmetric walk in $u_i=\log v_i$ respects positivity, but then the transformed target has a factor $v_i$; equivalently the acceptance ratio in variance coordinates has the extra factor $v_i^{\prime }/v_i$.

Hold all parameters except ${\mu }_i$ fixed and propose ${\mu }_i^{\prime }={\mu }_i+\eta$ with $\eta \sim N(0,s_{\mu i}^2)$. This random-walk Metropolis algorithm has a symmetric proposal distribution. DefineCancellation of unchanged likelihood and prior factors gives the Metropolis–Hastings acceptance probabilityDraw an independent uniform and accept if it is at most $a_i$; otherwise retain the old value. Update each component in turn, recomputing the needed mixture contributions. On a logarithmic scale, compare $\log U$ with $\min (0,{\sum }_j\log (D_j^{\prime }/D_j)-[({\mu }_i^{\prime }{)}^2-{\mu }_i^2]/(2{\tau }^2))$ to avoid numerical underflow.

For minimization of an energy $E(\vartheta )$, simulated annealing at temperature $t>0$ targets the Boltzmann distribution ${\pi }_t(\vartheta )\propto e^{-E(\vartheta )/t}$ relative to a specified parameter reference measure. With proposal distribution $q({\vartheta }^{\prime }\mid \vartheta )$, accept withAt each temperature run the chain long enough to explore that law, then cool and retain the best state seen. For symmetric proposals every improvement is accepted, while a worsening move is accepted with probability $e^{-\Delta E/t}$. Equilibrium laws concentrate near global minima when their normalizing integrals are finite and the minima are suitably isolated. A particular rapid cooling schedule does not automatically inherit this guarantee.

Likelihood maximization uses energy $E=-\ell$, not $E=\ell$. The source's later wording “$f(p)$ equal to the log-likelihood” conflicts with its preceding minimization convention. The intended maximum-likelihood Boltzmann distribution is proportional to $e^{\ell /t}=L^{1/t}$. With the opposite sign, likelihood maxima would be disfavored; for binomial data with both successes and failures the density often would not even be integrable near the endpoints as $t$ decreases.

For the maximum-likelihood estimates, write $s={\sum }_{j=1}^m{x}_j$, $\bar{x}=s/m$ and $S={\sum }_j(x_j-\bar{x}{)}^2$. In the binomial distribution model with known positive integer $n$, the log-likelihood isDifferentiating gives $s/p-(mn-s)/(1-p)=0$, hence $\hat{p}=s/(mn)$. The second derivative is negative in the interior. If $s=0$ or $s=mn$, the maximum is at the corresponding endpoint. This model requires integer observations in $\{0,\dots ,n\}$; otherwise its likelihood is zero.

For the normal distribution model write $v={\sigma }^2$. Its log-likelihood isFor fixed $v$, the unique maximizing mean is $\bar{x}$. Substituting it and differentiating in $v$ gives $v=S/m$, so, assuming $S>0$,If $S=0$, the normal density likelihood grows without bound as $v\downarrow 0$ with $\mu =\bar{x}$; there is no finite interior maximum or ordinary finite AIC for that fit.

For Gaussian likelihood Gibbs annealing, use $L_N^{1/t}$ relative to $d\mu dv$. Completing the square and identifying an inverse-gamma distribution gives the Gibbs sampler conditionalsIntegrating out $\mu$ contributes a factor proportional to $v^{1/2}$, so the marginal is $v\sim \mathrm{InvGamma}(m/(2t)-3/2,S/(2t))$. The joint target is proper for $0<t<m/3$ when $S>0$; starting the Lebesgue-based annealing above this range would require a proper reference factor or bounded parameter domain. For sufficiently small $t$, its moments satisfyThus equilibrium draws concentrate at $(\bar{x},S/m)$ as $t\downarrow 0$.

There is also a direct convergence proof for the actual successive Gibbs sampler updates here. Choose deterministic $0<t_r\le m/10$ with $t_r\to 0$, start at a finite deterministic $v_0>0$, draw ${\mu }_r$ from the first conditional using $v_{r-1}$, then draw $v_r$ from the second using ${\mu }_r$. Conditional normal moments give $E[({\mu }_r-\bar{x}{)}^2\mid v_{r-1}]=t_r{v}_{r-1}/m$ and $E[({\mu }_r-\bar{x}{)}^4\mid v_{r-1}]=3t_r^2{v}_{r-1}^2/m^2$. The permitted inverse-gamma moments therefore yieldThe coefficient on the preceding first moment is at most $1/6$, and that on the preceding second moment is at most $1/8$. The remaining coefficients are uniformly bounded, so both moment sequences are bounded by iteration. Since $t_r\to 0$, the recursions then give $E[v_r]\to S/m$ and $E[v_r^2]\to (S/m{)}^2$. Also $E[({\mu }_r-\bar{x}{)}^2]=t_{r}E[v_{r-1}]/m\to 0$. HenceThis proves convergence in mean square, and thus in probability, to the normal maximum-likelihood estimates for this explicit cooling/update scheme, rather than merely assuming equilibration.

For binomial likelihood-power annealing, normalizing $p^{s/t}(1-p{)}^{(mn-s)/t}$ on $(0,1)$ givesIts mean is $(s+t)/(mn+2t)$, and its variance isThe mean tends to $s/(mn)$ and the variance tends to zero, proving concentration at the binomial maximum-likelihood estimate, including the endpoint cases. Independent exact draws from these beta laws are a simple annealing implementation within this model.

For model selection, the Akaike information criterion is $2d-2\hat{\ell }$ with $d_B=1$ and $d_N=2$. ThusUse the endpoint convention $0\log 0=0$. For a trans-dimensional annealing for penalized likelihood construction, target the disjoint model spaces with unnormalized densitiesrelative to $dp$ and $d\mu dv$, respectively. The energies are half the unoptimized AIC, $d-\ell$; their joint minimizer is the fitted model with smaller AIC. Keep constants such as $C_B$ and the normal density normalization, which do not cancel across models. The within-model conditional draws are unchanged by these constant penalties.

Here is a complete dimension-changing proposal. Select a binomial-to-normal move with probability $b>0$ and a normal-to-binomial move with probability $d>0$, both fixed; other steps can use the within-model updates. From $B,p$, draw $u\sim N(0,1)$ with density $h(u)=(2\pi {)}^{-1/2}{e}^{-u^2/2}$ and setThis bijection from $(p,u)\in (0,1)\times \mathbb{R}$ onto $({\mu }^{\prime },v^{\prime })\in \mathbb{R}\times (0,\infty )$ has inverse $p=(1+e^{-{\mu }^{\prime }}{)}^{-1}$, $u=\log v^{\prime }$. Its absolute Jacobian determinant is $v^{\prime }/[p(1-p)]$. The Binomial-normal reversible-jump annealing acceptance ratio is thereforeFor a reverse proposal from $N,\mu ,v$, compute $p^{\prime }=(1+e^{-\mu }{)}^{-1}$ and $u^{\prime }=\log v$ deterministically and accept withThese ratios are reciprocal at inverse states and enforce detailed balance on the dimension-matched augmented spaces. They specify the proposal density, model-selection probabilities and Jacobian explicitly. The map is a valid example, not a claim of optimal mixing. Positive-temperature beta draws stay in the open interval even when the limiting estimate is an endpoint. At each temperature ensure adequate model switching; retain the best fitted penalized likelihood encountered. If minima tie, limiting concentration need not choose a unique model.

There is a material measurement qualification. A binomial probability mass and an unrounded normal density use different observation reference measures. Their raw AIC values are not an invariant, scientifically calibrated comparison of the same recorded data. The expressions and jump algorithm above implement the formal density comparison suggested by the question, with its units fixed. For integer-recorded measurements, a compatible normal model instead assigns probabilities $P_{\mu ,v}(x-1/2<Z\le x+1/2)$, suitably adjusting boundaries or truncation if required. Then both models describe the same recorded outcomes relative to counting measure. Use those masses in the normal likelihood and the same reversible-jump formula; its simple normal MLE and conjugate Gibbs updates will generally no longer apply. This observation-model correction cannot be silently combined with the continuous-normal formulas.

For a target probability density function $\pi$ and a proposal distribution with density $q(x^{\prime }\mid x)$, the Metropolis–Hastings algorithm starts in the target support and repeats: propose $x^{\prime }\sim q(\cdot \mid x)$, draw an independent uniform $U$, and move to $x^{\prime }$ ifOtherwise retain $x$. A zero reverse-proposal density makes the acceptance zero. Unknown multiplicative constants of $\pi$ cancel. For distinct states,which is symmetric in the states. Together with the probability of holding at the current state, this proves detailed balance and invariance of $\pi$. Successive values form a dependent Markov chain; convergence from arbitrary initial states additionally needs the appropriate irreducibility and aperiodicity, not just this invariant-density identity.

The Gibbs sampler updates a coordinate using its full conditional distribution while keeping the others fixed. For coordinate $i$, write $r=x_{-i}$ and $\pi (x)={\pi }_{-i}(r){\pi }_i(x_i\mid r)$. Its proposal replaces $x_i$ by $x_i^{\prime }\sim {\pi }_i(\cdot \mid r)$. The MH ratio, interpreted on this one-coordinate proposal space, isThus every Gibbs coordinate update is an MH update with acceptance one. A random choice of coordinate gives a mixture of these invariant kernels; a deterministic sweep gives their composition. Both preserve $\pi$. The latter composition need not be reversible, so the special-case statement is about the coordinate kernels, not a claim that every entire systematic sweep is one reversible MH kernel.

Use the random-walk proposal specified by its centre and covariance: for the current state $z=(x,y)$ draw $Z_1,Z_2$ independently standard normal and setIts properly normalized density of the proposal distribution isIt is symmetric: $q(z^{\prime }\mid z)=q(z\mid z^{\prime })$. This is a centred-at-the-current-state proposal, rather than an unshifted draw from the density notation in the source. The random-walk Metropolis algorithm consequently needs only the target ratio. Put $H(x,y)=x^2+y^2-2\rho xy$. ThenAn independent uniform decides the move; on rejection both coordinates remain at their previous values. A proposal with lower $H$ is always accepted. This construction uses the joint density directly, without needing either conditional distribution. Its continuous positive proposal density supplies communication over the target support, and the holding probability avoids a deterministic periodic transition pattern.

Let $\tilde{\pi }$ be a nonnegative target kernel and $q(y\mid x)$ a proposal distribution. At state $x$, propose $Y\sim q(\cdot \mid x)$, draw an independent uniform $U$, and move to $Y$ ifOtherwise keep $x$. This is the Metropolis–Hastings algorithm. Start on the positive target support, with the usual zero-ratio conventions. The unknown normalizing constant cancels.

An unnormalized kernel is sufficient, provided $0<\int \tilde{\pi }<\infty$. A genuinely improper target with infinite integral is not a probability density and cannot supply a stationary target probability distribution. The acceptance formula can still be written formally, but cannot be said to generate samples from that nonexistent probability law.

Because rejected proposals leave the state unchanged, the Markov kernel includes an atom:Detailed balance with a probability measure $\pi$ means the measure identityWhere both sides have ordinary densities, this reads $\pi (x)K(x,y)=\pi (y)K(y,x)$. The measure formulation also includes the rejection atom. Integrating it shows that $\pi$ is invariant.

Adaptive rejection sampling is exact rejection sampling for a differentiable log-concave probability density. Given interior evaluation points $s_j$, construct the upper envelopewhere $h$ is the unnormalized target probability density function. The concave function $\ell$ lies below all its tangent lines, so $e^u\ge h$. Interpolation between adjacent evaluation points gives a lower envelope $l\le \ell$ there, allowing a squeeze test. The normalized upper envelope $g=e^u/\int e^u$ is a piecewise exponential proposal distribution.

Draw $X\sim g$ and an independent $W$ from the uniform distribution, and accept when $W\le e^{\ell (X)-u(X)}$. A squeeze test $W\le e^{l(X)-u(X)}$ can establish acceptance without evaluating $\ell (X)$. Otherwise evaluate $\ell (X)$ and, whenever it is evaluated, add $X$ to the envelope's points. Repeat until acceptance. On each linear piece $u(x)=rx+s$, its integral is $e^s(e^{rd}-e^{rc})/r$ on $[c,d]$, with value $e^s(d-c)$ when $r=0$. These integrals select a piece; inverse transform sampling samples within it.

Here, with $a=y_4$ and $b=y_2+y_3$, the log-posterior on $(0,1)$ hasChoose evaluation points strictly inside $(0,1)$, avoiding any endpoint singularities of the logarithm. Every finite set of tangent lines gives a finite upper-envelope integral on this bounded interval; slope conditions needed for unbounded supports are unnecessary. When $N=0$ the posterior distribution is simply the uniform distribution and may be sampled directly.

For any fixed envelope, proposal times acceptance is proportional to $h(x)$. Conditional on all past evaluations and rejections, the next accepted value therefore has the same normalized posterior density. To see that adaptation does not spoil this, at every attempted draw its accepted mass is a scalar times that fixed posterior density; summing over possible rejection histories preserves it. The acceptance probability stays bounded below by its positive initial-envelope value because insertion only lowers the upper envelope. Thus acceptance occurs almost surely. With fresh independent randomness, successive accepted values are exact independent samples from the posterior distribution, even when the envelope is retained and improved between them.

Let ${\varphi }_d(z)=(2\pi {)}^{-d/2}\exp (-\Vert z{\Vert }^2/2)$ be the standard multivariate normal density. Independence gives the proposal distribution densityFor an orthogonal matrix $A$,The isotropic Gaussian density is unchanged by this transformation. Since $A^T$ has the same distribution as $A$,Thus this is an orthogonal-mixture Metropolis proposal with a symmetric density; inversion invariance, rather than any assumption of a uniform distribution on orthogonal matrices, is what is needed.

With $\alpha (x,y)=\min \{1,\pi (y)/\pi (x)\}$, the accepted off-diagonal transitions satisfyThe rejection part is supported on the diagonal and satisfies detailed balance automatically. Values of the acceptance rule at starting points where $\pi (x)=0$ can be specified separately; they do not affect stationarity under $\pi$. Therefore the Metropolis–Hastings algorithm kernel is reversible with invariant density $\pi$, proving

Let $\pi$ be the normalized density $\tilde{\pi }$ in the transformed coordinates of (iv). Its fixed Gaussian increment gives a symmetric proposal distribution: $q({\theta }^{\prime }\mid \theta )=q(\theta \mid {\theta }^{\prime })$. For distinct states, the accepted transition density is $K(\theta ,{\theta }^{\prime })=q({\theta }^{\prime }\mid \theta )\alpha (\theta ,{\theta }^{\prime })$, where $\alpha =\min \{1,\pi ({\theta }^{\prime })/\pi (\theta )\}$. ThenThe rejection probability gives a mass on the diagonal, which trivially satisfies the same balance identity. Thus $\boxed{\pi (d\theta )K(\theta ,d{\theta }^{\prime })=\pi (d{\theta }^{\prime })K({\theta }^{\prime },d\theta )}$, the detailed balance condition, and integrating it proves invariance of the target.

The transformed density includes its Jacobian, so mapping back preserves the desired posterior law. Warmup adaptation is excluded from this fixed-kernel proof. Independent augmentation with $H_0$ is itself a reversible conditional update; alternatively the original joint augmented space can be sampled directly. This is the Metropolis–Hastings algorithm balance argument.

Use broad proper uniform priors for $\Delta t$, $\Delta m$, and $c$ over physically plausible ranges, and broad log-uniform priors for the positive scales $A_f$ and ${\tau }_f$. ThenA random-walk Metropolis algorithm can update $(\Delta t,\Delta m,c,\log A_f,\log {\tau }_f)$ with a multivariate Gaussian proposal distribution. Initialize several dispersed chains near plausible cross-correlation delays and near the marginal-likelihood optimum; reject proposals outside the prior bounds; discard warm-up while adapting only the proposal scale and covariance; then freeze the kernel and retain a long run. Evaluate trace plots, autocorrelations, acceptance rates, between-chain agreement, and the effective sample size of a Markov chain. Posterior predictive quasar light curves provide a model check.

Write the target posterior density as $\pi (\theta )$ and the proposal distribution density as $q({\theta }^{\prime }\mid \theta )$. The Metropolis–Hastings algorithm accepts a proposed move withFor distinct states,which is symmetric in $\theta$ and ${\theta }^{\prime }$. The rejection probability supplies the diagonal part, so the entire transition kernel satisfies detailed balance. Integrating the detailed-balance identity over the starting state proves $\int \pi (\theta )P(\theta ,d{\theta }^{\prime })=\pi ({\theta }^{\prime })d{\theta }^{\prime }$. Hence the posterior is a stationary distribution; an irreducible Markov chain that is also an aperiodic Markov chain converges uniquely to it.

At equal sample count, rejection sampling eventually has smaller variance for fixed ${\sigma }^2>1$; at equal expected proposal cost, importance sampling eventually wins. The distinction matters because a sample of the rejection sampling estimator means an accepted draw.

For the importance sampling estimator, all coordinates contribute to the importance weight although only one enters the integrand. The same normal distribution transform as in (c), now with linear coefficient $-2$ for coordinate one, givesBoth estimators are unbiased estimators of $\sqrt{e}$, andFor fixed $s>1$, the first variance grows exponentially while the second stays constant. This is the benefit of returning exact target draws in rejection sampling. However, obtaining those draws costs $s^{d/2}$ proposals each on average.

To compare computational work, multiply each variance by its expected proposal count. For importance sampling that count is $n$; for rejection sampling it is $n{s}^{d/2}$. The ratio of the resulting cost-normalized variances isbecause $A/\sqrt{s}=\sqrt{s/(2s-1)}<1$. Thus importance sampling eventually gives better precision for the same expected number of proposals, assuming comparable work per proposal. This is an asymptotic cost comparison, not a guarantee that a fixed small computing budget produces any accepted draws. Both methods still deteriorate exponentially in this comparison. For $s=1$ they coincide and every proposal is accepted.

A practical improvement is to exploit the integrand's one-coordinate structure directly: average $e^{-G_1}$ for standard normal distribution draws $G_1$, without generating unused coordinates. The optimal change of proposal distribution in (iii) improves further.

Each state has exactly $k(n-k)$ 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 $q(v\mid u)=q(u\mid v)$. The discrete uniform distribution used as the prior distribution is also constant on $S_k$, with mass ${\left(\genfrac{}{}{0ex}{}{n}{k}\right)}^{-1}$. Both factors cancel in the Metropolis–Hastings acceptance probability.

Let $Q={\Sigma }^{-1}$ and $d=v-u$. The displayed multivariate normal distribution presupposes that the covariance matrix $\Sigma$ is a positive-definite matrix, so $Q$ 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 $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 toOnly two rows and four entries of $Q$ are needed for this evaluation. This is Metropolis sampling on fixed-cardinality subsets.

Choose $k$ distinct coordinates uniformly for the starting ones; each element of $S_k$ then has starting probability $1/\left(\genfrac{}{}{0ex}{}{n}{k}\right)$. 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$ or $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 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 $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: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.

To update a positive parameter, a symmetric additive proposal distribution can be generated on the whole real line. Reject nonpositive proposals; accept positive ones using the target density ratio. The off-diagonal proposal remains symmetric. Repeatedly redrawing until the proposal is positive instead truncates and renormalizes the proposal differently at different current states, so a Hastings correction is then required. A symmetric walk in $\log v$ is another valid method, but its transformed target includes the Jacobian determinant factor $v$.