From 47f23ea58a1a6dabeacf85f1e32bc7ebc96eb303 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Wed, 29 Jul 2026 08:55:29 -0700 Subject: [PATCH] mcsamplerGPU: fix adaptive-proposal support truncation biasing lnZ low On a mild 2D Gaussian test (T4b, floor_level unset), the adapted sampler was biased by -0.32 nats while every within-run error estimate read ~0.02: the 1D marginals of the drawn samples matched the claimed p_s exactly, but E[prior/p_s] = 0.37, i.e. ~63% of the prior volume sat in histogram bins with exactly zero proposal probability. A zero bin is an absorbing state (it can never be re-drawn), so the sampled support shrinks irreversibly and the integral silently loses the mass outside it. Four fixes: - compute_hist: clamp the uniform-mixture floor to HIST_FLOOR_LEVEL_MIN=1e-2 so no bin can reach zero probability (production ILE already passes 0.1; the clamp only binds for smaller/unset floors). - integrate_log: build adaptation weights from the stored tempered importance weights exp(tempering_exp*lnL + ln p - ln p_s) so the weighted histogram estimates the fixed target L^beta * prior -- the documented contract (see integrate() and the ILE driver comment). The old lnL+max(maxlnL,200) weights ignored tempering_exp and 1/p_s; being near-flat, each histogram replayed the previous proposal's sampling noise, a multiplicative random walk that collapsed the proposal onto a comb of surviving bins (61/100 per dim in the test) and drove the truncation above. - integrate_log: n_adapt freeze test double-multiplied by n (n_adapt was already scaled at parse time), so adaptation never froze regardless of the requested chunk count; also scale the no-kwarg default consistently. - pdf_from_hist: clamp bin index to n_bins-1 (right-edge sample previously indexed out of range). Validation (2D Gaussian, n=2000, neff=1000, 16 runs/config): bias -0.321 +/- 0.011 -> -0.003 +/- 0.004 for all n_adapt in {5..100} and tempering_exp in {0, 0.1, 1.0}; final proposal support 56% -> 100%; runs now reach neff~1000 instead of exhausting nmax at neff~350. E[prior/p_s] over fresh draws from the adapted proposal: 0.373 -> 1.0002. Co-Authored-By: Claude Fable 5 --- .../Code/RIFT/integrators/mcsamplerGPU.py | 30 ++++++++++++++----- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/integrators/mcsamplerGPU.py b/MonteCarloMarginalizeCode/Code/RIFT/integrators/mcsamplerGPU.py index b03ef004f..bff6f14fd 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/integrators/mcsamplerGPU.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/integrators/mcsamplerGPU.py @@ -93,6 +93,12 @@ def profile(fn): rosDebugMessages = True +# Minimum uniform-mixture fraction applied to every adapted histogram (see +# compute_hist): guarantees no bin has exactly zero sampling probability, since a +# zero bin can never be re-drawn (absorbing state) and silently truncates the +# integration domain. Override at module level for controlled experiments. +HIST_FLOOR_LEVEL_MIN = 1e-2 + class NanOrInf(Exception): def __init__(self, value): self.value = value @@ -297,7 +303,12 @@ def compute_hist(self, x_samples, param,weights=None,floor_level=0): # Smooth the histogram # kernel_size =3 # histogram_values = self.xpy.convolve( histogram_values, self.xpy.ones(kernel_size)/kernel_size,mode='same') - # Mix with a uniform sampling + # Mix with a uniform sampling. A bin with exactly zero probability is an + # absorbing state: it can never be drawn again, so the sampled support is + # permanently truncated and the integral is systematically biased LOW by the + # mass outside the support -- a bias no within-run error estimate can see. + # Enforce a minimal floor so every bin stays reachable. + floor_level = max(floor_level, HIST_FLOOR_LEVEL_MIN) histogram_values = histogram_values*(1-floor_level)+floor_level*self.xpy.ones(len(histogram_values))/len(histogram_values) # Evaluate the CDF by taking a cumulative sum of the histogram. @@ -344,7 +355,7 @@ def pdf_from_hist(self, x, param): y = (x - self.x_min[param]) / self.x_max_minus_min[param] # Compute the indices of the histogram bins that `x` falls into. indices = self.xpy.trunc(y / self.dx[param], out=y).astype(np.int32) - indices = self.xpy.minimum(indices,self.n_bins[param]) # prevent being out of range due to rounding ! + indices = self.xpy.minimum(indices,self.n_bins[param]-1) # prevent being out of range due to rounding (x == right edge maps to last bin) # Return the value of the histogram. return self.histogram_values[param][indices] @@ -645,7 +656,7 @@ def integrate_log(self, lnF, *args, xpy=xpy_default,**kwargs): print(" Note: cannot adapt, no history ") tempering_exp = kwargs["tempering_exp"] if "tempering_exp" in kwargs else 0.0 - n_adapt = int(kwargs["n_adapt"]*n) if "n_adapt" in kwargs else 1000 # default to adapt to 1000 chunks, then freeze + n_adapt = int(kwargs["n_adapt"]*n) if "n_adapt" in kwargs else 1000*n # default to adapt to 1000 chunks, then freeze. NOTE: scaled by n, matching integrate() floor_integrated_probability = kwargs["floor_level"] if "floor_level" in kwargs else 0 temper_log = kwargs["tempering_log"] if "tempering_log" in kwargs else False tempering_adapt = kwargs["tempering_adapt"] if "tempering_adapt" in kwargs else False @@ -815,8 +826,7 @@ def integrate_log(self, lnF, *args, xpy=xpy_default,**kwargs): # The total number of adaptive steps is reached # # FIXME: We need a better stopping condition here - if self.ntotal > n_adapt*n: - print(n_adapt,self.ntotal) + if self.ntotal > n_adapt: # n_adapt already scaled by n above; the old test (n_adapt*n) double-counted n and never froze continue # @@ -828,8 +838,14 @@ def inner(arg): return f(arg, p) return inner - weights_alt = self._rvs["log_integrand"][-n_history:]+np.max([maxlnL, 200]) # try to make sure we have some dynamic range here - weights_alt = self.xpy.maximum(weights_alt, 1e-5) # prevent negative weights. NOTE THIS IS IMPORTANT: if you are integrating a function with lnL<0, use an offset! + # Tempered importance weights exp(tempering_exp*lnL + ln p - ln p_s), so the + # weighted histogram of draws estimates the FIXED target L^tempering_exp * prior. + # (The old lnL + max(maxlnL,200) weights ignored tempering_exp and the 1/p_s + # correction: near-flat weights made each histogram replay the previous + # proposal's sampling noise, a multiplicative random walk that collapses the + # proposal onto a comb of surviving bins.) + weights_alt = self._rvs["log_weights"][-n_history:] + weights_alt = self.xpy.exp(weights_alt - self.xpy.max(weights_alt)) weights_alt = weights_alt/(weights_alt.sum()) if weights_alt.dtype == RiftFloat: weights_alt = weights_alt.astype(numpy.float64,copy=False)