Adaptive MCMC resamples the proposal until it lands in the box, then accepts with a plain Metropolis ratio — the truncation's Hastings factor is missing, so the chain systematically under-samples near the bounds

[wshlavacek]

What happens

pick_new_pset proposes with while num != 10000*len_params: around a fresh delta_vector = self.chain_rngs[idx].multivariate_normal(...) draw, applies it with oldpset.get_param(k.name).add(delta_vector_add[k.name], False) (reflect=False), and on except OutOfBoundsException: num += 1; pass simply draws a brand-new delta and tries again. FreeParameter.set_value (pybnf/pset.py:2409-2411) raises OutOfBoundsException exactly when new_value < self.lower_bound or new_value > self.upper_bound and reflect is False. So the effective proposal is the Gaussian renormalized over the box: q(x->y) = phi(y-x; C) / Z(x), with Z(x) = the proposal mass inside the box from x. Z(x) depends on x, so q is NOT symmetric. But got_result accepts with a bare Metropolis ratio and no Hastings term: self.alpha[index] = np.exp((lnposterior-self.ln_current_P[index])) / if self.chain_rngs[index].random() < self.alpha[index] (lines 184-185) — no log(Z(y)/Z(x)). The stationary density of that kernel is pi(x)Z(x), not pi(x) (substitute pi_s(x)=cZ(x) into the balance equation: int c*Z(x)phi(y-x)/Z(x) dx over the box = cZ(y)). Z is smallest at the walls, so am — the sampler the docs call PyBNF's recommended Bayesian method — is biased away from the box boundary. The same pattern is in the pre-adaptive branch (lines 568-586) and the continue_run branch (lines 547-563). For contrast, DREAM gets this right: _calculate_de_pset uses add(..., False) with the comment "Do not reflect the parameter (need to reject if outside bounds)" and return None, cr_idx so the chain stays put (a valid rejection), and basic_mcmc uses reflect=True (a symmetric fold). Only am resamples. (After 10000 failed draws it switches to reflect=True, so the bias is confined to the normal path, which is the path every proposal takes.)

Reproduction

Direct simulation of the two kernels on a flat target (uniform_var p on [0,1], constant likelihood, step 0.2 — i.e. the pre-adaptive else branch of pick_new_pset with step_size=0.2, population_size=1).

Script at /private/tmp/claude-503/-Users-l119605-Code-PyBNF/eb55c5a3-1006-4c5d-9d1e-9d14fe8045f0/scratchpad/k.py, run with uv run python .../k.py from /Users/l119605/Code/PyBNF:

rng=np.random.default_rng(0); s=0.2; x=0.5; N=400000
for i in range(N):
    while True:
        y=x+rng.normal(0,s)
        if 0.0<=y<=1.0: break   # am: resample until in box (reflect=False)
    x=y                         # flat posterior -> always accept
    out[i]=x

Observed, density-normalized 10-bin histogram:
resample-until-in-box (am): [0.708 0.910 1.055 1.140 1.171 1.171 1.145 1.072 0.913 0.715]
reflect (basic_mcmc, mh/pt): [0.994 0.991 0.996 0.993 1.008 1.007 1.008 1.000 1.000 1.003]

Expected for both: flat at 1.000 (the target is exactly Uniform(0,1)).
The am kernel is bowed — 0.708 in the outermost bin vs 1.171 in the center, a ~1.65x deficit at the walls, matching the predicted stationary density proportional to Z(x)=Phi((1-x)/0.2)-Phi(-x/0.2). Far outside Monte Carlo error at N=4e5. The reflecting kernel that mh/pt use is flat to within 1%, confirming the bias is specific to the reject-and-redraw proposal and not to the acceptance test.

To see it inside PyBNF rather than in the abstract: fit_type=am, one uniform_var on [0,1], an objfunc returning a constant score, population_size=1, step_size=0.2, small burn_in, max_iterations 1e6, sample_every=1 — the recorded histogram and the credible68/95 intervals from combine_chains_params come out correspondingly bowed and too narrow.

Verification notes

I tried to refute this three ways and it survived all three.

1. Is the path reachable? Yes. FreeParameter.set_value (/Users/l119605/Code/PyBNF/pybnf/pset.py:2399-2411) raises OutOfBoundsException exactly when new_value < self.lower_bound or new_value > self.upper_bound and reflect is False. lower_bound/upper_bound are finite whenever self.bounded (pset.py:2266-2267, 2295-2296), which is the case for uniform_var/loguniform_var — the standard variable type for a Bayesian fit. All three branches of pick_new_pset in /Users/l119605/Code/PyBNF/pybnf/algorithms/samplers/adaptive_mcmc.py (the adaptive branch ~line 509-524, the continue_run branch ~line 547-563, the pre-adaptive else branch ~line 566-586) start with num = 0, and every one of them selects reflect=False while num < 10000. So the very first proposal — and in practice every proposal — takes the reject-and-redraw path. On except OutOfBoundsException: num += 1; pass the while body re-enters, new_vars is reset to [], and a brand-new multivariate_normal delta is drawn. It is a fresh draw, not a retry of the same one.

2. Is there a correction downstream? No. grep -n "correction\|log_corr\|ln_corr" pybnf/algorithms/samplers/adaptive_mcmc.py returns nothing. got_result (adaptive_mcmc.py:176-185) accepts with a bare ratio: self.alpha[index] = np.exp((lnposterior-self.ln_current_P[index])) then if self.chain_rngs[index].random() < self.alpha[index]. No log(Z(y)/Z(x)) anywhere. For contrast, dream.py:375-383 and :482 explicitly build and apply a log_correction Hastings term, so the codebase knows the concept and applies it where it belongs.

3. Is the math right in context anyway? No, and the sibling sampler's own comment is the clincher. basic_mcmc.py:287-291 says, verbatim, that FreeParameter.add defaults to reflect=True so a step leaving the box "is folded back inside rather than rejected. The fold is symmetric, so the plain Metropolis ratio in got_result still targets the correct bound-restricted posterior." That is the documented justification for mh/pt using a bare ratio — and it is precisely the property am gives up by passing reflect=False. Detailed balance for q(x->y)=phi(y-x)/Z(x) with alpha=min(1,pi(y)/pi(x)) holds for pi(x)Z(x), not pi(x): LHS = phimin(pi(x),pi(y))/Z(x), RHS = phimin(pi(x),pi(y))/Z(y), equal only if Z(x)=Z(y).

4. Does a test pin the current behavior as correct? No. tests/test_adaptive_mcmc.py has 23 tests; the proposal ones (test_random_walk_is_gaussian_step_size, test_proposal_is_current_plus_scaled_draw_from_adapted_covariance, test_seeded_from_chain_sample_covariance) pin the shape and scale of the Gaussian step and the covariance adaptation, and the acceptance ones (test_worse_proposal_alpha_is_metropolis_ratio, test_acceptance_frequency_matches_alpha) pin the bare Metropolis ratio in isolation. Nothing tests the two together against a known target, which is exactly why the mismatch went unnoticed. No test asserts the truncated proposal is symmetric.

I then ran the kernel directly and it reproduces the predicted bow. One correction to the claim: severity should be medium, not high. In the adapted regime the proposal scale is diff*sqrt(diffMatrix) = (2.38^2/d)*posterior_sd, so Z(x)≈1 except within a couple of proposal sds of a wall — the bias is a thin boundary shell, not a global distortion, and it affects only am (mh/pt reflect, dream/pdream reject-and-stay with a Hastings term, hmc reparameterizes). It is still a real defect, because parameters railing against a bound are common in this domain and the shell is exactly where the reported credible interval's edge lives; it also worsens with dimension, since the whole vector is redrawn when any one component leaves the box.

This is committed code, not part of the dirty working tree: git status --porcelain pybnf/algorithms/samplers/adaptive_mcmc.py is empty.

Where

pybnf/algorithms/samplers/adaptive_mcmc.py:518 — severity medium, bug class wrong-math.