Paper deep dive
Estimating Great Expectations Under Autoregressive Language Models with Potentials
Francesco I. Re, Shubhangi Ghosh, Tim Vieira, Ryan Cotterell
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 10/12/2026, 5:57:55 AM
Summary
This paper introduces a variance-reduction framework for estimating expectations of test functionals under autoregressive language models. By utilizing next-token conditional probabilities available during sampling, the authors define 'potentials'—real-valued functions on prefixes that decompose test functionals additively. This approach, termed the potential estimator, leverages control variates to significantly reduce variance compared to standard Monte Carlo and importance sampling methods across various applications including loss estimation, divergence calculation, and watermark detection.
Entities (9)
Relation Signals (7)
Potential Estimator → utilizes → Next-token conditional probabilities
confidence 96% · make estimation more efficient by exploiting the next-token conditional probabilities which are available as a by-product of sampling.
Potential Estimator → reducesvarianceof → Monte Carlo Estimation
confidence 95% · Our experiments demonstrate substantial gains in statistical efficiency... yielding up to 17x the effective sample size of Monte Carlo.
Oracle Potential → yields → Zero Variance
confidence 94% · The oracle potential... satisfies... Var_p(hat_theta(phi^*)) = 0.
Potential Estimator → reducesvarianceof → Importance Sampling
confidence 92% · learned potentials selected through nested cross-fitting reduce the variance of the importance-sampling estimator... by a median 12%
Potential Estimator → appliesto → Entropy
confidence 90% · For losses and divergences, simple potentials reduce variance... Section 4.1 illustrates this for entropy.
Potential Estimator → appliesto → KL Divergence
confidence 90% · We estimate seven losses... KL divergence... simple potentials reduce variance in all tested configurations
Potential Estimator → appliesto → Watermark Detection
confidence 88% · our estimator reduces the variance of watermark-detection estimates by a median 36%
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Many applications of language models hinge not on individual samples but on the expectation of a test functional under the model. Estimating such expectations reliably can be computationally expensive. In this paper, we show how to make estimation more efficient by exploiting the next-token conditional probabilities which are available as a by-product of sampling. We do so through potentials: real-valued functions on prefixes that decompose the test functional additively. We construct an estimator whose variance depends on the chosen potential, and derive conditions under which a potential reduces this variance. We then develop practical potentials for several estimands and applications, and demonstrate substantial variance reductions across several estimands at comparable computational cost.
Tags
Links
- Source: https://arxiv.org/abs/2610.11399v2
- Canonical: https://arxiv.org/abs/2610.11399v2
Trouble viewing inline? Open PDF directly →
Full Text
140,516 characters extracted from source content.
Expand or collapse full text
Estimating Great Expectations under Autoregressive Language Models with Potentials Francesco Ignazio Re Affiliation: ETH Zürich Email: mailto:francesco.re@inf.ethz.chfrancesco.re@inf.ethz.ch Shubhangi Ghosh Affiliation: Columbia University Email: mailto:sg4156@columbia.edusg4156@columbia.edu Tim Vieira Affiliation: ETH Zürich Email: mailto:tim.vieira@gmail.comtim.vieira@gmail.com Ryan Cotterell Affiliation: ETH Zürich Email: mailto:ryan.cotterell@inf.ethz.chryan.cotterell@inf.ethz.ch Abstract Many applications of language models hinge not on individual samples but on the expectation of a test functional under the model. Estimating such expectations reliably can be computationally expensive. In this paper, we show how to make estimation more efficient by exploiting the next-token conditional probabilities which are available as a by-product of sampling. We do so through potentials: real-valued functions on prefixes that decompose the test functional additively. We construct an estimator whose variance depends on the chosen potential, and derive conditions under which a potential reduces this variance. We then develop practical potentials for several estimands and applications, and demonstrate substantial variance reductions across several estimands at comparable computational cost. Code Repository 1 Introduction A language model is a probability distribution over strings. The most common application of such a model is to generate free-form text from a prompt (Radford et al., 2019; Brown et al., 2020). However, for some applications, the quantity of interest is not the sample itself but rather the expectation of a test functional under the model. For instance, • In reinforcement learning from human feedback, proximal policy optimization (Schulman et al., 2017, PPO;) regularizes the policy by penalizing its Kullback–Leibler divergence from a reference model (Christiano et al., 2017; Stiennon et al., 2020; Ouyang et al., 2022), a divergence that is itself an expectation under the policy. • In AI safety, one is often interested in the probability that a model emits a harmful or otherwise undesirable output (Jones et al., 2025; Angell et al., 2026). • In distillation, one fits a simpler surrogate, such as an n-gram model, to a large language model by matching expectations such as the cross-entropy under the teacher (Hinton et al., 2015; Agarwal et al., 2024). • In language-model watermarking, one may be interested in the probability that a generated text is detected as watermarked; that can be viewed as an expectation under the model’s distribution (Kirchenbauer et al., 2023). • In interpretability, distributional functionals quantify how much a given context shifts a model’s predictions (Pezeshkpour, 2023; Du et al., 2024; Stoehr et al., 2024). More precisely, given a language model p and a test function f, we wish to compute the expectation of f under p. Unless p and f exhibit exploitable structure11 1 When p is a probabilistic finite-state automaton (Rabin, 1963) or a probabilistic context-free grammar (Booth & Thompson, 1973) and f is additively decomposable, one can compute such expectations with dynamic programming (Eisner, 2002; Li & Eisner, 2009; Zmigrod et al., 2021)., this expectation must be approximated. Monte Carlo methods offer a well known approach: sample strings, evaluate f on them, and average the results. Drawing independent strings directly from p yields the standard Monte Carlo estimator; drawing instead from a proposal distribution and reweighting by the corresponding likelihood ratio yields importance sampling. These estimators are unbiased but can have high variance, so achieving a desired level of precision may require generating many samples at substantial computational cost. In this paper, we introduce a variance-reduction framework that exploits the next-token distributions computed as a by-product of sampling. The idea is to write f as an additive decomposition of token-level increments along the string and to replace each increment with its conditional expectation given the prefix, taken over the next token. Applied to an existing unbiased estimator, this procedure always preserves unbiasedness while making the variance depend on the decomposition. The choice of decomposition is thus the central design choice in the framework. For a given f, we parametrize the family of its additive decompositions using potentials: real-valued functions on prefixes. We borrow this construction from potential-based reward shaping in reinforcement learning, where potentials redistribute rewards along a trajectory to facilitate learning (Ng et al., 1999). Here, we repurpose potentials to redistribute the function value over marginalizable token level increments and refer to the resulting estimator as the potential estimator. We first show that there exists an oracle potential yielding a zero-variance estimator. Since it is generally intractable, we derive an exact criterion for when an arbitrary potential reduces variance. We then use the oracle potential to guide potential design: when f has known additive structure, we use components of the oracle potential that can be computed directly from f; for arbitrary functions over strings, we approximate the oracle potential by learning it from samples. In general, given a set of candidate potentials, we introduce a potential-selection procedure that preserves unbiasedness while using the full sample. Our experiments demonstrate substantial gains in statistical efficiency across four applications and several estimands. For losses and divergences, simple potentials reduce variance in all tested configurations, yielding up to 17×17× the effective sample size of Monte Carlo. Their gradient counterparts also improve on Monte Carlo in every evaluated comparison, with variance reductions of up to 89%89\%. For expected n-gram counts and frequencies, our estimator produces even larger gains: about 99%99\% lower variance than Monte Carlo after accounting for computational cost (median over targets with lower estimation uncertainty). For non-additive event probabilities, approximations to the oracle potential are also effective: our estimator reduces the variance of watermark-detection estimates by a median 36%36\%, while learned potentials selected through nested cross-fitting reduce the variance of the importance-sampling estimator introduced by (Dorman et al., 2026) by a median 12%12\%. In the latter experiment, this corresponds to approximately 900,000900,000 additional baseline generations. Overall, our results establish potential design as a practical strategy for improving expectation estimates after sampling, across both direct Monte Carlo and importance sampling. We discuss related work in Appendix A. 2 Expectations under Autoregressive Language Models 2.1 Language Models Let Σ be a finite, non-empty set of symbols, and let eos∉Σ eos∉ be a distinguished terminal symbol. A string is a finite sequence of symbols from Σ . Let Σ∗ ^* denote the (countably infinite) set of all such strings. For ,′∈Σ∗ y, y ∈ ^*, we write ⪰′ y y when ′ y is a prefix of y, and ′≺ y y when ′ y is a strict prefix of y. A language model p is a probability distribution over complete strings taking values in Σ∗ ^*, i.e., p∈△(Σ∗)p∈ ( ^*), where △() (X) denotes the set of all probability distributions over a countable set X. We write ∼p Y p for a random string drawn from p, so that ⪰ Y y is the event that Y extends y. For any distribution p, we use the following notation for the prefix probability of a string ∈Σ∗ y∈ ^*, which gives the marginal probability that a string sampled from p extends y: p→()=defPr[⪰]=∑∈Σ∗p(). p( y) 0.5pt def=Pr[ Y y]= _ z∈ ^*p( y z). (1) The conditional prefix probability is then p→(∣)=defPr[⪰∣⪰]=p→()/p→() p( z y) 0.5pt def=Pr[ Y y\, z Y y]= p( y\, z)/ p( y), and we extend p→(⋅∣) p(· y) to the terminal symbol eos by p→(eos∣)=defp()/p→() p( eos y) 0.5pt def=p( y)/ p( y). Whenever p→()=0 p( y)=0, we use the convention that p→(∣)=def0 p( z y) 0.5pt def=0 for every ∈Σ z∈ and p→(eos∣)=def1 p( eos y) 0.5pt def=1, so that p→(⋅∣) p(· y) remains a valid conditional distribution over Σ∪eos ∪\ eos\. Using the conditional prefix probabilities, the probability of any string =y1⋯y||∈Σ∗ y=y_1·s y_| y|∈ ^* under a language model admits a left-to-right autoregressive factorization: p()=p→(eos∣)∏n=1||p→(yn∣<n),p( y)= p( eos y)\, _n=1^| y| p(y_n y_<n), (2) where <n=defy1⋯yn−1 y_<n 0.5pt def=y_1·s y_n-1. Each p→(⋅∣<n) p(· y_<n) in Eq. 2 is a probability distribution over Σ∪eos ∪\ eos\. 2.2 Problem Statement Let p∈△(Σ∗)p∈ ( ^*) be a language model and let f:Σ∗→ℝf ^* be a real-valued test function. We study the problem of estimating θ=def∼p[f()]=∑∈Σ∗p()f(),θ 0.5pt def= E _ Y p[f( Y)]= _ y∈ ^*p( y)f( y), (3) assuming throughout that θ is finite. For modern language models, this expectation is rarely available in closed form. The string space Σ∗ ^* is countably infinite, and summing exactly over all strings is generally infeasible. A standard solution is to approximate such expectations using Monte Carlo methods. We review the canonical Monte Carlo estimator and importance sampling in the remainder of this section. In the next section, we introduce potential-based estimation, an additional estimation stage that can reduce the variance of these estimators after generation by reusing the sampled strings and the model’s next-token distributions retained at sampling time. 2.3 Monte Carlo Estimation The most direct estimator of θ is an empirical average over sampled strings. We assume throughout the paper that p[f()2]<∞ E _p\! [f( Y)^2 ]<∞ and let (1),…,(M)∼i.i.d.p Y^(1),…, Y^(M) i.i.d. p . The Monte Carlo estimator and its variance are θ^MC=def1M∑m=1Mf((m)),Varp(θ^MC)=1M(p[f()2]−θ2). θ_MC 0.5pt def= 1M _m=1^Mf( Y^(m)), _p( θ_MC)= 1M ( E _p\! [f( Y)^2 ]-θ^2 ). (4) Because the samples are drawn i.i.d. from p and θ is finite, θ^MC θ_MC is unbiased, [θ^MC]=θ E [ θ_MC]=θ. Moreover, the strong law of large numbers guarantees that θ^MC θ_MC is strongly consistent: it converges almost surely to θ as M→∞M→∞ (Owen, 2013; Robert & Casella, 2004). 2.4 Importance Sampling Importance sampling estimates an expectation under p using samples from a proposal distribution q. Assume that q()>0q( y)>0 whenever p()|f()|>0p( y)|f( y)|>0. For ∈Σ∗ y∈ ^*, define fq()=defp()f()/q()f_q( y) 0.5pt def=p( y)f( y)/q( y) when q()>0q( y)>0, and fq()=def0f_q( y) 0.5pt def=0 otherwise. Then θ=∼p[f()]=∼q[fq()]θ= E _ Y p[f( Y)]= E _ Y q[f_q( Y)]. Thus, for M independent samples (1),…,(M)∼i.i.d.q Y^(1),…, Y^(M) i.i.d. q, the importance-sampling estimator is simply the Monte Carlo estimator of fqf_q under q, such that θ^IS=def1M∑m=1Mp((m))q((m))f((m)),Varq(θ^IS)=1M(q[fq()2]−θ2). θ_IS 0.5pt def= 1M _m=1^M p( Y^(m))q( Y^(m))f( Y^(m)), _q( θ_IS)= 1M ( E _q[f_q( Y)^2]-θ^2 ). (5) This estimator is unbiased for every such proposal distribution q, but its variance may be lower or higher than that of θ^MC θ_MC in Equation 4 (Robert & Casella, 2004; Owen, 2013). Its variance is small when the weighted contribution p()q()f() p( Y)q( Y)f( Y) varies little under q. When p[|f()|]>0 E _p[|f( Y)|]>0, the variance-minimizing proposal is q⋆()∝defp()|f()|q ( y) 0.5pt def p( y)|f( y)|, which assigns probability in proportion to each string’s absolute contribution to the expectation and makes the magnitude of the weighted contribution constant. Proposition 3 makes this precise. 3 Estimating Expectations with Potentials We seek to construct a lower-variance estimator by exploiting the target model’s next-token distributions along sampled strings, which are already computed when sampling from p or evaluating p()p( Y) for importance weighting under a proposal q. For notational simplicity, we assume q=pq=p throughout this section; the general case q≠pq≠ p is given in Section B.5 (B.5 Potentials under a Proposal). 3.1 Additive decompositions over prefixes: potentials Our construction begins with a potential, a prefix-level function ϕ→:Σ∗→ℝ φ ^* assigning a real value to each prefix. We use its increments over successive prefixes to define token-level contributions along each sampled string. Specifically, for the fixed test function f, we define the increment operator Δ , which maps a potential ϕ→ φ to its increment function Δϕ→ φ, given by Δϕ→(a∣)=defϕ→(a)−ϕ→()(a∈Σ),Δϕ→(eos∣)=deff()−ϕ→(). φ(a y) 0.5pt def= φ( ya)- φ( y) (a∈ ), φ( eos y) 0.5pt def=f( y)- φ( y). (6) For notational ease, we write Δϕ→(a∣) φ(a y) for [Δϕ→](a∣)[ φ](a y), and use the same convention whenever an operator acts on a function. Using this notation, given the random string ∼p Y p and its length ||| Y|, we can rewrite f()f( Y) as f()=ϕ→(ε)+∑t=1||Δϕ→(Yt∣<t)+Δϕ→(eos∣).f( Y)= φ( )+ _t=1^| Y| φ(Y_t Y_<t)+ φ( eos Y). (7) This representation expresses f()f( Y) as a sum of random token-level increments Δϕ→(Yt∣<t) φ(Y_t Y_<t) along the sampled string Y. At each realized prefix <t y_<t, we already have the next-token distribution p→(⋅∣<t) p(· y_<t), which provides the probabilities needed to compute the next increment’s conditional expectation by averaging over possible next tokens. In the next subsection, we show how such expected value can be used for variance reduction. 3.2 The Potential Estimator We use the Bellman operator to map a potential to its expected value over the next token: ϕ→()=def∑a∈Σp→(a∣)ϕ→(a)+p→(eos∣)f().B φ( y) 0.5pt def= _a∈ p(a y)\, φ( ya)+ p( eos y)\,f( y). (8) The expected potential increment at prefix y is then Δϕ→()=defA∼p→(⋅∣)[Δϕ→(A∣)]=ϕ→()−ϕ→().B φ( y) 0.5pt def= E _A p(· y) [ φ(A y) ]=B φ( y)- φ( y). (9) Substituting these expectations in Eq. 7 gives the following estimator. Definition 1 (Potential estimator). Let (1),…,(M)∼i.i.d.p Y^(1),…, Y^(M) i.i.d. p. Given a potential ϕ→ φ, define θ^(ϕ→)=defϕ→(ε)+1M∑m=1M∑t=1|(m)|+1Δϕ→(<t(m)),Varp(θ^(ϕ→))=1MVarp(∑t=1||+1Δϕ→(<t)). θ( φ) 0.5pt def= φ( )+ 1M _m=1^M _t=1^| Y^(m)|+1B φ( Y^(m)_<t), _p ( θ( φ) )= 1MVar_p\! ( _t=1^| Y|+1B φ( Y_<t) ). (10) The inner sum includes the final prefix (m) Y^(m), at which eos is drawn. We assume that the sum of potential increments along a sampled string has finite second moment. Under this condition, the estimator is unbiased (Proposition 4), while its variance depends on the choice of potential. This raises the natural question of when a potential reduces variance, and how much reduction is possible. We answer this question in the next proposition, which identifies a zero-variance potential, which we call the oracle potential, and characterizes exactly when an arbitrary potential improves on the unbiased baseline. Proposition 1 (Oracle potential and variance comparison). The oracle potential, defined on prefixes of positive probability, satisfies ϕ→⋆()=p[f()∣⪰]⟹θ^(ϕ→⋆)=θ,Varp(θ^(ϕ→⋆))=0. φ ( y)= E _p\! [f( Y) Y y ] θ( φ )=θ, _p\! ( θ( φ ) )=0. Let θ^∗ θ_* be any unbiased baseline estimator under p. For every fixed potential ϕ→ φ, Varp(θ^(ϕ→))≤Varp(θ^∗)⇔Covp(θ^∗,θ^∗−θ^(ϕ→))≥12Varp(θ^∗−θ^(ϕ→)).Var_p ( θ( φ) ) _p( θ_*) _p ( θ_*, θ_*- θ( φ) )≥ 12Var_p ( θ_*- θ( φ) ). Despite the above condition need not hold for a chosen potential, we can exploit the unbiasdness of θ^∗ θ_* and θ^(ϕ→) θ( φ) to build a control-variate estimator with variance equal or lower than θ^∗ θ_*. Definition 2 (Potential control-variate estimator). Given a potential ϕ→ φ and an unbiased baseline estimator θ^∗ θ_*, define θ^(ϕ→,α)=defθ^∗−α(θ^∗−θ^(ϕ→)),α∈ℝ. θ( φ,α) 0.5pt def= θ_*-α ( θ_*- θ( φ) ), α . (11) The coefficient α controls the contribution of the correction. Setting α=0α=0 recovers the baseline estimator, while α=1α=1 recovers the potential estimator. As is standard for control variates, we choose this coefficient to minimize variance. The next proposition gives the optimal choice. Proposition 2 (Unbiasedness and optimal coefficient). For every fixed potential ϕ→ φ and fixed coefficient α, the estimator is unbiased: p[θ^(ϕ→,α)]=θ E _p[ θ( φ,α)]=θ. If both θ^∗ θ_* and θ^(ϕ→) θ( φ) have finite second moments and Varp(θ^∗−θ^(ϕ→))>0Var_p ( θ_*- θ( φ) )>0, the variance-minimizing coefficient satisfies α⋆=Covp(θ^∗,θ^∗−θ^(ϕ→))Varp(θ^∗−θ^(ϕ→))⟹Varp(θ^(ϕ→,α⋆))≤Varp(θ^∗).α = Cov_p ( θ_*, θ_*- θ( φ) )Var_p ( θ_*- θ( φ) ) _p ( θ( φ,α ) ) _p( θ_*). (12) If the difference has zero variance, the two estimators agree almost surely, and the choice of coefficient has no effect. Having characterized the role of the potential and the optimal coefficient, we next turn to their practical construction and estimation. In the next section, we describe how to design and select potentials and estimate parameters such as α∗α^* from the same samples, without sacrificing unbiasedness. 4 Potential Design in practice: applications and experiments In this section, we show how to construct potentials across several estimands and applications, and ask whether the resulting potential estimators achieve lower variance at comparable computational cost. Depending on the structure of f, we follow one of two complementary approaches: we either exploit known additive structure to derive a potential analytically, or learn a potential from data. We present three applications in the main text and a fourth in Section C.1 (C.1 Analysis of Watermark Detection). 4.1 Estimating proper losses, divergences, and gradients Table 1: Entropy f()=−logp()f( Y)=- p( Y) ϕ→⋆()=−logp→()+p[−logp()p→()|⪰] aligned φ ( [rgb]0,0,0 x)&=- p( [rgb]0,0,0 x)\\ & + E _p\! [- p( Y) p( [rgb]0,0,0 x)\; |\; Y [rgb]0,0,0 x ] aligned ϕ→()=−logp→() φ( [rgb]0,0,0 x)=- p( [rgb]0,0,0 x) Δϕ→(a∣)=−logp→(a∣) φ(a [rgb]0,0,0 x)=- p(a [rgb]0,0,0 x) Δϕ→(eos∣)=−logp→(eos∣) φ( eos [rgb]0,0,0 x)=- p( eos [rgb]0,0,0 x) We estimate seven losses frequently used in language model training and evaluation along with their gradients: entropy (Dai et al., 2018; Braverman et al., 2020), KL divergence (Schulman et al., 2017; Ouyang et al., 2022), Brier loss and power-score regrets (Shao et al., 2024) with exponents β∈1.25,1.5,1.75,2β∈\1.25,1.5,1.75,2\. Defining a potential analytically using the function value at the observed prefix. The above losses and their gradients have a natural additive structure over prefixes; for any such function, the value of the function at each observed prefix yields a potential. Section 4.1 illustrates this for entropy. For a complete string ⪰ y [rgb]0,0,0 x, the chain rule splits −logp()- p( y) into a prefix term −logp→()- p( [rgb]0,0,0 x) and a continuation term −log(p()/p→())- (p( y)/ p( [rgb]0,0,0 x)), which includes termination. The oracle potential adds the expected continuation term to the prefix term; keeping only the prefix term thus approximates the oracle potential by its observed component. Here the increments are the per-token surprisals −logp→(a∣)- p(a [rgb]0,0,0 x) (with a final increment at eos), which are computed during sampling anyway, so the potential is free. We define the potentials for the other quantities and their gradients in the appendix: Section C.2.1 gives those for KL divergence, the Brier loss and the power-score regrets (Sections C.2.1, C.2.1 and C.2.1), and Section C.3.1 those for the gradients (Section C.3.1). For KL, our potential estimator coincides with that of Amini et al. (2025); we reproduce their experiment in Section C.2.4 (C.2.4 Reproducing Amini et al. (2025)). We still include KL in our evaluation because their claimed variance ordering relative to Monte Carlo does not hold in general, as shown by Proposition 1 and the counterexample in Section C.2.4. Across all quantities and gradients, we therefore test how consistently the estimator reduces variance and whether fitting the control-variate coefficient provides additional gains. 4.1.1 Experiments and results (a) variance ratio per quantity (b) DKLD_KL ratio vs. estimated DKLD_KL (c) cross-fitted control variate Figure 1: Variance reduction in the proper-loss experiment (lower is better). (a) Variance ratio of the potential estimator to Monte Carlo across power scores (β→1β→ 1 gives DKLD_KL); entropy and Brier are shown on the right. Each series is one policy–reference pair (see 10): solid lines are the four main policy–reference pairs and dashed lines are consecutive Gemma-2 SPPO pairs. (b) For DKLD_KL, variance ratio against the estimated DKL(p∥q)D_KL(p\|q), with one colored marker per pair and grey markers for individual prompts from all pairs. The shaded regions indicate prompts below and above the Monte Carlo variance. (c) Variance ratio of the cross-fitted control-variate estimator to the potential estimator; hollow markers use the full-sample oracle coefficient. The potential estimator reduces variance in every configuration, by up to 17×17× in effective sample size. We evaluate the estimator on eight configurations from three post-training families (full details in Section C.2.2, C.2.2 Experimental setup). Across all 5656 pair–quantity combinations, the potential reduces variance relative to Monte Carlo (Figure 1a), with ratios between 0.0580.058 and 0.9870.987 and every bootstrap upper bound below one. These correspond to effective sample-size multipliers of 1.01×1.01×–17.3×17.3×. The gains are largest when the policy and reference are close. The variation across policy–reference pairs is larger than across estimands. For KL at horizon 256256, the variance ratio is 0.860.86 on GPT-2 IMDb, 0.770.77 on Tülu-3, and 0.320.32 on Gemma-2, dropping to 0.060.06–0.110.11 between consecutive Gemma SPPO iterations (equivalent to 9.1×9.1×–16.7×16.7× the Monte Carlo sample size). The same pattern appears across prompts: smaller KL is associated with larger variance reductions (Figure 1b). Variance reduction holds for the vast majority of prompts, but not all of them: Monte Carlo has lower variance on 1.5%1.5\% of all 4,0964,096 prompts, including 12.3%12.3\% of Tülu-3 prompts (95% interval 9.6%9.6\%–15.2%15.2\%). This provides empirical evidence that the general KL variance ordering of Amini et al. (2025) does not hold universally; Fitting the control-variate coefficient yields ratios of 0.860.86–1.001.00 relative to the untuned potential estimator (Figure 1c); among KL prompts where the potential estimator increases variance, fitting α lowers the median variance ratio relative to Monte Carlo from 1.1341.134 to 0.9740.974. Section C.3 repeats the comparison for the gradient of each estimand with respect to the policy parameters along a DPO trajectory of the GPT-2 sentiment policy: early in training, the potential removes most of the gradient variance; when it helps little on its own, fitting the control-variate coefficient recovers substantial gains. 4.2 Estimating Expected Counts and Frequencies of n-grams We next study expected counts and frequencies of word n-grams, i.e., contiguous sequences of n words. These statistics arise in applications such as model distillation (Suresh et al., 2021; Krishnan et al., 2023) and computational psycholinguistics (Opedal et al., 2024; Kiegeland et al., 2026). More broadly, n-gram counts are one instance of a general class of functions f()=#occurrences of substring in f( Y)=\#\occurrences of substring w in Y\, encompassing quantities such as the number of forbidden keywords in a string, unsafe API calls, or other keyword-based patterns emitted by a model. Designing potentials that marginalize over multiple future tokens. We start by observing that, for any target substring ∈Σ∗ w∈ , the count function f admits two natural decompositions: f()=∑tZ,tend=∑tZ,tstart,Z,tend= ends at t,Z,tstart= starts at t.f( Y)= _tZ^end_ w,t= _tZ^start_ w,t, Z^end_ w,t= 1\ w ends at t\, Z^start_ w,t= 1\ w starts at t\. Table 2: Substring counts, lookahead occurrence potential f()=#occurrences of in f( Y)=\#\occurrences of w in Y\ f()=∑tZ,tstartf( Y)= _tZ^start_ w,t, Z,tstart=Z,tcan+Z,tnoncanZ^start_ w,t=Z^can_ w,t+Z^noncan_ w,t ϕ→⋆()=p[∑t=1||Z,tstart|⪰] φ ( [rgb]0,0,0 x)= E _p\! [ _t=1^| Y|Z^start_ w,t\; |\; Y [rgb]0,0,0 x ] ϕ→()=∑t≤||p[Z,tcan∣≤t=≤t] φ( [rgb]0,0,0 x)= _t≤| [rgb]0,0,0 x| E _p\! [Z^can_ w,t Y_≤ t= [rgb]0,0,0 x_≤ t ] Δϕ→(a∣)=p[Z,tcan∣≤t=a] φ(a [rgb]0,0,0 x)= E _p\! [Z^can_ w,t Y_≤ t= [rgb]0,0,0 xa ], t=||+1t=| [rgb]0,0,0 x|+1 Δϕ→(eos∣)=∑t≤||(Z,tcan−p[Z,tcan∣≤t])+∑tZ,tnoncan φ( eos [rgb]0,0,0 x)= _t≤| [rgb]0,0,0 x| (Z^can_ w,t- E _p[Z^can_ w,t Y_≤ t] )+ _tZ^noncan_ w,t The ZendZ^end_ w decomposition yields a potential that is completely determined by evaluating f at prefix x. Its expected increment, however, is nonzero only at prefixes that already match the substring up to its final token, so marginalization is useful only at those few prefixes. By contrast, ZstartZ^start_ w is generally still unobserved at the current prefix. We therefore replace it by its full conditional expectation—the probability that the occurrence starting at that position is eventually completed—and take the cumulative expected count over starting positions as the potential. In practice, we marginalize over the canonical tokenization =(w1,…,wL) w=(w_1,…,w_L) with wi∈Σw_i∈ and L=def||L 0.5pt def=| w|. At each sampled prefix, this requires L additional token evaluations: the remaining target tokens after the next token, plus the word boundary. Thus, if C denotes the cost of generating one sample from model p, the additional marginalization cost is O(LC)O(LC). Including generation gives a total per-target token-evaluation cost proxy of (L+1)C(L+1)C. Section 4.2 gives the mathematical definition of the potential using ZstartZ^start_ w and its increments. This example illustrates the flexibility of the telescoping construction in Equation 7: the potential may marginalize only a selected component of f—here, occurrences under the canonical tokenization—while the terminal increment recovers the remainder, so the decomposition still equals f exactly at termination.22 2 In this application, part of the terminal increment has expectation zero under p and can therefore be dropped without changing the estimator’s expectation, yielding Δ~ϕ→(eos∣)=∑tZ,tnoncan φ( eos [rgb]0,0,0 x)= _tZ^noncan_ w,t. 4.2.1 Experiments & Results Setup. We construct 7,6477,647 word n-gram targets (n=1,…,6n=1,…,6) from the English texts of MECO (Siegelman et al., 2022), with canonical tokenization lengths ranging from 11 to 1010 tokens. We ask whether the potential estimator that uses the potential defined in § 4.2 reduces variance at comparable computational cost with MC, and how the gains vary with target rarity; (a) ρ^c(w) ρ_c(w) with log-ratio 95% CI (b) ρ^c(w) ρ_c(w) against target rarity (c) resolved by token length Figure 2: Cost-adjusted variance ratio ρ^c(w) ρ_c(w) for n-gram frequency estimation. (a) Empirical cumulative distribution over 3,7003,700 targets with bootstrap RSE at most 22, with approximate 95%95\% log-ratio confidence bounds. (b) Variance reduction versus estimated frequency for the 2,8312,831 targets with at least 2525 expected occurrences in the independent frequency bank. (c) Fraction of targets whose frequency estimate reaches each RSE threshold, grouped by token length, at the same sample size M=4,000M=4,000 for both estimators. Dashed parity lines in (a–b) mark equal variance at equal cost under the proxy; in (c), line styles distinguish RSE thresholds. For each target, we estimate its expected count and frequency (the ratio estimator of Equation 29, Section C.4.1) from 4,0004,000 GPT-2 samples using the potential estimator. We estimate the Monte Carlo variance from an independent 4×1064× 10^6-sample bank drawn from the same model. We then report the cost-adjusted variance ratio ρ^c(w) ρ_c(w) (Equation 30, Section C.4.2), scaling the potential-estimator variance by the per-target cost proxy Lw+1L_w+1 and the Monte Carlo variance to 4,0004,000 samples. Thus, ρ^c(w)<1 ρ_c(w)<1 indicates lower variance at equal cost under this proxy. To measure target rarity independently, we estimate each target’s frequency from a separate 2×1062× 10^6-sample bank. We report frequency results here because they are easier to interpret, and defer count results to Section C.4; both show the largest gains for rare targets. The appendix also reports results for the simpler observed-prefix potential (Section 4.1), which provides virtually no variance reduction. Half the retained targets reach equal variance on at most 1.3%1.3\% of the Monte Carlo compute. Panel Figure 2a reports the cost-adjusted variance ratio ρ^c ρ_c for the 3,7003,700 of 7,6477,647 targets estimated with relative standard error at most 22, with bootstrap 95%95\% confidence intervals. The median ratio is 1.24%1.24\%, and 74.7%74.7\% of targets require at most 10%10\% of the Monte Carlo compute. Using the upper confidence bound instead gives a median of 2.56%2.56\%, with 95.6%95.6\% of targets still below parity. This shows that the efficiency gains remain substantial even under a conservative treatment of uncertainty. At the same sample size (M=4,000M=4,000), the potential estimator yields substantially more precise frequency estimates, with 79.7%79.7\% of targets reaching RSE ≤20%≤ 20\% and 96.5%96.5\% reaching RSE ≤50%≤ 50\%, compared with 11.4%11.4\% and 19.0%19.0\% for Monte Carlo (see Figure 2c). The gap grows with token length: for targets of length 44 or more, Monte Carlo almost never reaches RSE ≤50%≤ 50\%, whereas the potential estimator does so for at least 90%90\% of targets in every length group. For Lw≥7L_w≥ 7, Monte Carlo does not observe any target even in the 4×1064× 10^6-sample reference bank, while the potential estimator gives a finite estimate for all of them. 4.3 Event analysis with language models Dorman et al. (2026) introduce a framework for estimating event probabilities over any function using an importance-sampling estimator: they sample from distributions tilted toward different function values and reweight the resulting trajectories to estimate probabilities under the original language model. Given a function f on strings, they partition its range into intervals DjD_j and estimate the probability that f()f( Y) falls in each interval. The corresponding event indicators are fj()=f()∈Djf_j( Y)= 1\f( Y)∈ D_j\. We augment their estimator of p[fj()] E _p[f_j( Y)] with potentials for fjf_j. Learning potentials from samples. We note that these indicators need not admit a useful additive decomposition over prefixes. A simple baseline is the terminal potential, which assigns zero to unfinished prefixes and the indicator value fj()f_j( Y) at the string termination. To exploit information from earlier prefixes, we instead learn an approximation to the oracle potential: the conditional probability that fj()=1f_j( Y)=1 given the observed prefix. Each sampled prefix i [rgb]0,0,0 x_i, i=1,…,Bi=1,…,B, is paired with its trajectory’s outcome fj(i)f_j( y_i), i⪰i y_i [rgb]0,0,0 x_i. Using prefix features ()∈ℝd ξ( [rgb]0,0,0 x) ^d, we fit a predictor33 3 When trajectories are obtained through importance sampling as in this application, we weight each prefix’s prediction loss by its trajectory’s importance weight, so that fitting targets conditional event probabilities under the original language model. gj(⋅,):ℝd→[0,1]g_j(·; η) ^d→[0,1] with parameters ∈ℝm η ^m. We can define the potential ϕ→j(,)=gj((),)≈ϕ→j⋆()=Prp[fj()=1∣⪰]. φ_j( [rgb]0,0,0 x; η)=g_j( ξ( [rgb]0,0,0 x); η)≈ φ_j ( [rgb]0,0,0 x)=Pr_p[f_j( Y)=1 Y [rgb]0,0,0 x]. (13) Unbiased estimation with learned potentials. Fitting and evaluating potentials on the same samples can introduce bias. We avoid this through K-fold cross-fitting (Chernozhukov et al., 2018): each potential is trained on the other folds and evaluated on the held-out fold. Averaging the fold-specific estimates preserves unbiasedness under independent sampling and exact importance weights (when needed), while using all samples. Potential selection with nested cross-fitting. Different potentials may achieve different variance reductions, motivating data-driven selection. However, selecting a potential on the same samples used for estimation can introduce bias (Varma & Simon, 2006). We therefore use nested cross-fitting (Cawley & Talbot, 2010): for each held-out fold, we select the potential through inner cross-validation on the remaining data, then apply it to the held-out fold. We detail this procedure in Algorithm 1. 4.3.1 Experiments & Results Setup. We replicate the main experiment of Dorman et al. (2026), estimating the distribution of the Automated Readability Index (ARI; Equation 31), a sequence-level readability score, over texts generated by TinyStories-8M. Here f=ARIf=ARI. During sampling, we additionally retain the next-token distributions for subsequent potential fitting and evaluation. We consider three candidate potentials: the terminal potential, which marginalizes over the final token, and two learned potentials based on a simple logistic regression and a mixture model, respectively. We defer their details to Section C.5.2 (C.5.2 Potential definitions). We also include a control-variate variant of the terminal potential with a learned coefficient α, giving four candidates, among which we select using nested cross-fitting with K=5K=5 outer folds. Throughout, we refer to the original estimator of Dorman et al. (2026) by the acronym MBAR. MBAR uses estimated weights and is asymptotically unbiased under its regularity assumptions (Shirts & Chodera, 2008); we do not claim exact finite-sample unbiasedness for this experiment. Learned potentials reduce variance by ∼12% 12\%, equivalent to ∼0.9×106 0.9× 10^6 additional MBAR samples. Estimator Var. ratio vs. MBAR [95% CI] Fitted selection 0.88 [0.86, 0.91] Terminal 0.92 [0.90, 0.94] Terminal, α 0.92 [0.90, 0.95] Logistic 0.90 [0.86, 0.91] Mixture model 0.87 [0.84, 0.90] Table 3: variance relative to MBAR, median over the 6464 populated ARI bins; 95%95\% intervals from 1,0001,000 paired bootstrap replicas. Full version in 19. Table 3 summarizes the variance ratio of each potential estimator relative to MBAR; lower values indicate a more efficient estimator. The learned potentials give the largest reductions: the logistic and mixture-model potentials attain median variance ratios of 0.900.90 and 0.870.87, respectively, and nested selection attains 0.880.88, closely matching the best individual candidate. The terminal potential gives a smaller but still consistent reduction, with a median ratio of 0.920.92. It exploits cases where membership in a bin DjD_j is determined only at the final token. The experiment uses N=6.68×106N=6.68× 10^6 generations, matching the original experiments. Assuming variance scales inversely with sample size, the median ratio of 0.880.88 corresponds to an MBAR budget requiring 0.9×1060.9× 10^6 more generations than our potential estimator. We used only CPUs for potential estimation, as detailed in Section C.5.4 (C.5.4 Compute costs and additional results). 5 Conclusion We introduced a general framework for reducing the variance of expectation estimates under autoregressive language models using next-token distributions available as a by-product of sampling. By repurposing potentials from reinforcement learning, we construct an unbiased estimator whose variance depends on the chosen potential. Across several estimands and applications, we showed how practical potentials can improve statistical efficiency at comparable computational cost. More broadly, our results establish potential design as a general strategy for constructing lower-variance estimators under language models. AI use statement In this work, we used generative AI tools to provide feedback on research methodology and experimental design, assist with the interpretation of experimental results, and help check and refine mathematical statements and derivations. Additionally, we used generative AI tools to assist with software development, and drafting and editing portions of the manuscript for clarity, structure, and readability. AI-assisted text and code were reviewed and verified by the authors. We take responsibility for the final content of this work, including text, claims, and artifacts produced with the aid of generative AI. References Agarwal et al. (2024) Rishabh Agarwal, Nino Vieillard, Yongchao Zhou, Piotr Stanczyk, Sabela Ramos Garea, Matthieu Geist, and Olivier Bachem. On-policy distillation of language models: Learning from self-generated mistakes. In International Conference on Learning Representations, volume 2024, p. 21246–21263, 2024. Amini et al. (2025) Afra Amini, Tim Vieira, and Ryan Cotterell. Better estimation of the Kullback–Leibler divergence between language models. In Advances in Neural Information Processing Systems (NeurIPS), 2025. URL https://proceedings.neurips.c/paper_files/paper/2025/hash/a2e9be3d64bcc196731d0a4e0e2c33ae-Abstract-Conference.html. Angell et al. (2026) Rico Angell, Raghav Singhal, Zachary Horvitz, Zhou Yu, Rajesh Ranganath, Kathleen McKeown, and He He. Estimating tail risks in language model output distributions. arXiv preprint arXiv:2604.22167, 2026. Booth & Thompson (1973) Taylor L Booth and Richard A Thompson. Applying probability measures to abstract languages. IEEE transactions on Computers, 100(5):442–450, 1973. Boyd et al. (2022) Alex Boyd, Samuel Showalter, Stephan Mandt, and Padhraic Smyth. Predictive querying for autoregressive neural sequence models. Advances in Neural Information Processing Systems, 35:23751–23764, 2022. Braverman et al. (2020) Mark Braverman, Xinyi Chen, Sham Kakade, Karthik Narasimhan, Cyril Zhang, and Yi Zhang. Calibration, entropy rates, and memory in language models. In International Conference on Machine Learning, p. 1089–1099. PMLR, 2020. Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in neural information processing systems, 33:1877–1901, 2020. Cawley & Talbot (2010) Gavin C Cawley and Nicola LC Talbot. On over-fitting in model selection and subsequent selection bias in performance evaluation. The Journal of Machine Learning Research, 11:2079–2107, 2010. Chernozhukov et al. (2018) Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters, 2018. Christiano et al. (2017) Paul F Christiano, Jan Leike, Tom Brown, Miljan Martic, Shane Legg, and Dario Amodei. Deep reinforcement learning from human preferences. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL https://proceedings.neurips.c/paper_files/paper/2017/file/d5e2c0adad503c91f91df240d0cd4e49-Paper.pdf. Cochran (1977) William G. Cochran. Sampling Techniques. John Wiley & Sons, New York, 3 edition, 1977. Dai et al. (2018) Zihang Dai, Qizhe Xie, and Eduard Hovy. From credit assignment to entropy regularization: Two new algorithms for neural sequence prediction. In Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), p. 1672–1682, 2018. Del Moral (2004) Pierre Del Moral. Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications. 2004. URL https://doi.org/10.1007/978-1-4684-9393-1. Dorman et al. (2026) Jake McAllister Dorman, Edward Gillman, Dominic C Rose, Jamie F Mair, and Juan P Garrahan. Rare event analysis of large language models. arXiv preprint arXiv:2602.06791, 2026. Du et al. (2024) Kevin Du, Vésteinn Snæbjarnarson, Niklas Stoehr, Jennifer White, Aaron Schein, and Ryan Cotterell. Context versus prior knowledge in language models. In Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), p. 13211–13235, 2024. Eisner (2002) Jason Eisner. Parameter estimation for probabilistic finite-state transducers. In Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics, p. 1–8, 2002. Hammersley (1956) J. M. Hammersley. Conditional Monte Carlo. Journal of the ACM, 3(2), 1956. URL https://doi.org/10.1145/320825.320827. Hinton et al. (2015) Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. Distilling the knowledge in a neural network. arXiv preprint arXiv:1503.02531, 2015. Jones et al. (2025) Erik Jones, Meg Tong, Jesse Mu, Mohammed Mahfoud, Jan Leike, Roger Grosse, Jared Kaplan, William Fithian, Ethan Perez, and Mrinank Sharma. Forecasting rare language model behaviors. arXiv preprint arXiv:2502.16797, 2025. Kiegeland et al. (2026) Samuel Kiegeland, Vésteinn Snæbjarnarson, Tim Vieira, and Ryan Cotterell. On the proper treatment of units in surprisal theory. In Proceedings of the 64th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), p. 32202–32224, 2026. Kirchenbauer et al. (2023) John Kirchenbauer, Jonas Geiping, Yuxin Wen, Jonathan Katz, Ian Miers, and Tom Goldstein. A watermark for large language models. In International conference on machine learning, p. 17061–17084. PMLR, 2023. Krishnan et al. (2023) Aravind Krishnan, Jesujoba Alabi, and Dietrich Klakow. On the n-gram approximation of pre-trained language models. arXiv preprint arXiv:2306.06892, 2023. Lew et al. (2023) Alexander K. Lew, Tan Zhi-Xuan, Gabriel Grand, and Vikash K. Mansinghka. Sequential Monte Carlo steering of large language models using probabilistic programs. In ICML Workshop on Sampling and Optimization in Discrete Space (SODS), 2023. URL https://openreview.net/forum?id=Ul2K0qXxXy. Li & Eisner (2009) Zhifei Li and Jason Eisner. First-and second-order expectation semirings with applications to minimum-risk training on translation forests. In Proceedings of the 2009 Conference on Empirical Methods in Natural Language Processing, p. 40–51, 2009. Liu et al. (2026) Zixuan Liu, Fangzheng Wu, Brian Summa, and Zizhan Zheng. Adaptive multilevel twisted sequential Monte Carlo for rare events estimation in language models, 2026. URL https://arxiv.org/abs/2608.21736. Loula et al. (2025) João Loula, Benjamin LeBrun, Li Du, Ben Lipkin, Clemente Pasti, Gabriel Grand, Tianyu Liu, Yahya Emara, Marjorie Freedman, Jason Eisner, Ryan Cotterell, Vikash Mansinghka, Alexander K. Lew, Tim Vieira, and Timothy J. O’Donnell. Syntactic and semantic control of large language models via sequential Monte Carlo. In International Conference on Learning Representations (ICLR), 2025. URL https://openreview.net/forum?id=xoXn62FzD0. Ng et al. (1999) Andrew Y Ng, Daishi Harada, and Stuart Russell. Policy invariance under reward transformations: Theory and application to reward shaping. In Icml, volume 99, p. 278–287. Citeseer, 1999. Opedal et al. (2024) Andreas Opedal, Eleanor Chodroff, Ryan Cotterell, and Ethan Wilcox. On the role of context in reading time prediction. In Yaser Al-Onaizan, Mohit Bansal, and Yun-Nung Chen (eds.), Proceedings of the 2024 Conference on Empirical Methods in Natural Language Processing, p. 3042–3058, Miami, Florida, USA, November 2024. Association for Computational Linguistics. doi: 10.18653/v1/2024.emnlp-main.179. URL https://aclanthology.org/2024.emnlp-main.179/. Ouyang et al. (2022) Long Ouyang, Jeffrey Wu, Xu Jiang, Diogo Almeida, Carroll Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, et al. Training language models to follow instructions with human feedback. Advances in neural information processing systems, 35:27730–27744, 2022. Owen (2013) Art B. Owen. Monte Carlo theory, methods and examples. https://artowen.su.domains/mc/, 2013. Pezeshkpour (2023) Pouya Pezeshkpour. Measuring and modifying factual knowledge in large language models. In 2023 international conference on machine learning and applications (ICMLA), p. 831–838. IEEE, 2023. Rabin (1963) Michael O Rabin. Probabilistic automata. Information and control, 6(3):230–245, 1963. Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al. Language models are unsupervised multitask learners. OpenAI blog, 1(8):9, 2019. Rieder (2012) Helmut Rieder. Robust Asymptotic Statistics: Volume I, volume 1. Springer Science & Business Media, 2012. Robert & Casella (2004) Christian P. Robert and George Casella. Monte Carlo Statistical Methods. 2nd edition, 2004. URL https://doi.org/10.1007/978-1-4757-4145-2. Schulman et al. (2015) John Schulman, Philipp Moritz, Sergey Levine, Michael Jordan, and Pieter Abbeel. High-dimensional continuous control using generalized advantage estimation. arXiv preprint arXiv:1506.02438, 2015. Schulman et al. (2017) John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347, 2017. Shao et al. (2024) Chenze Shao, Fandong Meng, Yijin Liu, and Jie Zhou. Language generation with strictly proper scoring rules. arXiv preprint arXiv:2405.18906, 2024. Shirts & Chodera (2008) Michael R. Shirts and John D. Chodera. Statistically optimal analysis of samples from multiple equilibrium states. The Journal of Chemical Physics, 129(12):124105, 2008. doi: 10.1063/1.2978177. Siegelman et al. (2022) Noam Siegelman, Sascha Schroeder, Cengiz Acartürk, Hee-Don Ahn, Svetlana Alexeeva, Simona Amenta, Raymond Bertram, Rolando Bonandrini, Marc Brysbaert, Daria Chernova, et al. Expanding horizons of cross-linguistic research on reading: The multilingual eye-movement corpus (meco). Behavior research methods, 54(6):2843–2863, 2022. Stiennon et al. (2020) Nisan Stiennon, Long Ouyang, Jeffrey Wu, Daniel Ziegler, Ryan Lowe, Chelsea Voss, Alec Radford, Dario Amodei, and Paul F Christiano. Learning to summarize with human feedback. Advances in neural information processing systems, 33:3008–3021, 2020. Stoehr et al. (2024) Niklas Stoehr, Kevin Du, Vésteinn Snæbjarnarson, Robert West, Ryan Cotterell, and Aaron Schein. Activation scaling for steering and interpreting language models. In Findings of the Association for Computational Linguistics: EMNLP 2024, p. 8189–8200, 2024. Suresh et al. (2021) Ananda Theertha Suresh, Brian Roark, Michael Riley, and Vlad Schogol. Approximating probabilistic models as weighted finite automata. Computational Linguistics, 47(2):221–254, 2021. Varma & Simon (2006) Sudhir Varma and Richard Simon. Bias in error estimation when using cross-validation for model selection. BMC bioinformatics, 7(1):91, 2006. Whiteley & Lee (2014) Nick Whiteley and Anthony Lee. Twisted particle filters. The Annals of Statistics, 42(1), 2014. URL https://doi.org/10.1214/13-AOS1167. Wiewiora (2003) Eric Wiewiora. Potential-based shaping and Q-Value initialization are equivalent. Journal of Artificial Intelligence Research, 19, 2003. URL https://doi.org/10.1613/jair.1190. Zhao et al. (2024) Stephen Zhao, Rob Brekelmans, Alireza Makhzani, and Roger Baker Grosse. Probabilistic inference in language models via twisted sequential Monte Carlo. In Proceedings of the International Conference on Machine Learning (ICML), 2024. URL https://proceedings.mlr.press/v235/zhao24c.html. Zmigrod et al. (2021) Ran Zmigrod, Tim Vieira, and Ryan Cotterell. Efficient computation of expectations under spanning tree distributions. Transactions of the Association for Computational Linguistics, 9:675–690, 2021. Appendix Contents Limitations Our framework has several limitations. First, variance reduction depends on the choice of potential. Although the oracle potential yields zero variance under direct sampling, practical potentials only approximate it and are not guaranteed to improve over the baseline for every individual prompt. The control-variate formulation provides a safeguard against poorly chosen potentials, although estimating its coefficient or selecting among learned potentials requires cross-fitting to preserve unbiasedness. Thus, the strong performance of the potentials studied here should not be interpreted as implying that arbitrary potentials will yield similar gains. Second, the computational advantage depends on the cost of evaluating the potential. In several of our applications, the required quantities are already available during sampling or can be computed relatively cheaply, whereas more elaborate potentials may require additional model evaluations or fitting. At the same time, this tradeoff suggests a natural direction for future work: designing potentials that achieve stronger variance reduction while remaining computationally inexpensive. Third, our approach requires access to the model’s next-token probability distribution along sampled prefixes. It therefore applies directly when logits or token probabilities are available, but not to black-box generation APIs that expose only sampled outputs. Finally, our empirical conclusions are necessarily limited to the estimands, models, and experimental configurations considered here; the magnitude of the variance reductions may differ in other settings. Appendix A Related Work Importance sampling and sequential Monte Carlo for language models. Several methods reduce variance by modifying the sampling distribution to favor events of interest and correcting through importance weighting (Angell et al., 2026; Dorman et al., 2026). For example, Boyd et al. (2022) combine query-specific importance sampling with beam search to estimate probabilities of sequence events. Sequential Monte Carlo (Del Moral, 2004) and its twisted variants (Whiteley & Lee, 2014) have also been applied to controlled generation (Lew et al., 2023; Zhao et al., 2024; Loula et al., 2025) and rare-event estimation (Liu et al., 2026). In such methods, a twist assigns a nonnegative weight to each prefix, guiding sampling and resampling toward prefixes judged more likely to produce the desired outcome. For event-probability estimation, the optimal twist weights each prefix according to the probability that completing it under the original language model will satisfy the event. It is therefore proportional to our oracle potential ϕ→⋆ φ , which represents the same conditional probability (Zhao et al., 2024). Our framework uses prefix potentials to construct additive, zero-mean corrections that enable variance reduction after sampling, exploiting the next-token distribution created at generation. We demonstrate its complementarity with importance sampling by augmenting the importance sampling estimator of Dorman et al. (2026), who sample from exponentially tilted sequence distributions using Markov chain Monte Carlo. Conditional Monte Carlo, Rao–Blackwellization, and control variates. Replacing a finite-variance estimator by its conditional expectation cannot increase its variance, the principle underlying conditional Monte Carlo (Hammersley, 1956) and Rao–Blackwellization (Robert & Casella, 2004; Owen, 2013). For KL estimation under language models, Amini et al. (2025) replace each sampled token log-likelihood ratio by the exact KL divergence between the next-token distributions at the sampled prefix, and derive an analogous gradient estimator. Although their KL estimator is unbiased, the variance-reduction guarantee in their Theorem 2 does not hold in general: averaging individual token contributions can alter their covariances and increase the variance of the sum. Our framework recovers their KL estimator as a particular choice of potential and extends the construction to arbitrary test functions, explicitly distinguishing unbiasedness from variance reduction. The difference between a potential estimator and any unbiased baseline defines a zero-mean control variate (Owen, 2013). This correction preserves unbiasedness for any fixed coefficient, while the population-optimal coefficient yields variance no greater than that of the baseline. Potentials and reward shaping in reinforcement learning. Potential-based reward shaping in reinforcement learning uses changes in a function of the state to redistribute rewards while preserving optimal policies (Ng et al., 1999; Wiewiora, 2003). Value functions, which predict future rewards, provide natural potentials: the resulting shaped rewards measure prediction errors and are used in generalized advantage estimation (Schulman et al., 2015). Our construction uses the same telescoping principle, with text prefixes as states and the test function’s value as the terminal reward. We then average the potential increments over the next-token distribution to construct unbiased expectation estimates after sampling. The ideal potential predicts the expected terminal reward given each prefix; its averaged increments vanish, yielding our zero-variance oracle (Proposition 1). Appendix B Additional Propositions and Proofs Throughout, ∼p Y p denotes a single sampled string, ||| Y| its length excluding eos, and <t Y_<t its prefix of length t−1t-1, so that <1=ε Y_<1= and <||+1= Y_<| Y|+1= Y. All expectations are assumed finite. B.1 Importance Sampling and the Optimal Proposal Proposition 3 (Importance-sampling variance and optimal proposal). Let q be a proposal satisfying q()>0q( y)>0 whenever p()|f()|>0p( y)|f( y)|>0, and let θ^IS θ_IS be the importance-sampling estimator of Eq. 5. Then q[θ^IS]=θ E _q[ θ_IS]=θ and Varq(θ^IS)=1M(∑∈Σ∗p()2f()2q()−θ2).Var_q( θ_IS)= 1M ( _ y∈ ^* p( y)^2f( y)^2q( y)-θ^2 ). (14) Terms with q()=0q( y)=0 are read as zero. If p|f()|>0 E _p|f( Y)|>0, the variance is minimized uniquely by q⋆()=p()|f()|p|f()|,Varq⋆(θ^IS)=1M((p|f()|)2−θ2),q ( y)= p( y)\,|f( y)| E _p|f( Y)|, _q ( θ_IS)= 1M ( ( E _p|f( Y)| )^2-θ^2 ), (15) and the minimum is zero if and only if f has constant sign on the support of p. Proof. The summands of θ^IS θ_IS are i.i.d., so it suffices to treat M=1M=1. Using the support assumption, q[p()q()f()] E _q\! [ p( Y)q( Y)f( Y) ] =∑:q()>0q()p()q()f()=∑∈Σ∗p()f()=θ, = _ y:\,q( y)>0q( y) p( y)q( y)f( y)= _ y∈ ^*p( y)f( y)=θ, where the middle equality uses that p()f()=0p( y)f( y)=0 whenever q()=0q( y)=0. The same computation for the second moment gives q[(p()f()/q())2]=∑p()2f()2/q() E _q[(p( Y)f( Y)/q( Y))^2]= _ yp( y)^2f( y)^2/q( y), which yields Eq. 14. For the optimal proposal, write c=defp|f()|=∑p()|f()|c 0.5pt def= E _p|f( Y)|= _ yp( y)|f( y)|. By the Cauchy–Schwarz inequality applied to the vectors (p()|f()|/q()) (p( y)|f( y)|/ q( y) )_ y and (q()) ( q( y) )_ y over the support of q, c2=(∑p()|f()|q()q())2≤∑p()2f()2q()∑q()=∑p()2f()2q(). c^2= ( _ y p( y)|f( y)| q( y) q( y) )^2≤ _ y p( y)^2f( y)^2q( y) _ yq( y)= _ y p( y)^2f( y)^2q( y). Hence Varq(θ^IS)≥(c2−θ2)/MVar_q( θ_IS)≥(c^2-θ^2)/M for every admissible q, with equality if and only if the two vectors are proportional, that is, q()=p()|f()|/c=q⋆()q( y)=p( y)|f( y)|/c=q ( y) on the support of q; since q⋆q sums to one, equality forces q=q⋆q=q . Substituting q⋆q into Eq. 14 confirms the value in Eq. 15. Finally, c2−θ2=(c+|θ|)(c−|θ|)≥0c^2-θ^2=(c+|θ|)(c-|θ|)≥ 0 because c=p|f()|≥|pf()|=|θ|c= E _p|f( Y)|≥| E _pf( Y)|=|θ|, with equality if and only if f()f( Y) has almost surely constant sign under p. ∎ B.2 Unbiasedness of the Potential Estimator Proposition 4 (Unbiasedness and Variance). Let ϕ→ φ be any potential whose increments are absolutely summable in expectation, p[∑t=1||+1|Δϕ→(Yt∣<t)|]<∞ E _p [ _t=1^| Y|+1| φ(Y_t Y_<t)| ]<∞. Then the potential estimator of Eq. 10 is unbiased: p[θ^(ϕ→)]=θ E _p[ θ( φ)]=θ. Proof. It suffices to treat M=1M=1. Write T=def||+1T 0.5pt def=| Y|+1 for the index of the eos step, so that the event t≤T\t≤ T\ is determined by <t Y_<t: it holds exactly when <t Y_<t has not yet been terminated. For each t, conditioning on <t Y_<t and using that Yt∼p→(⋅∣<t)Y_t p(· Y_<t) on this event, p[t≤TΔϕ→(Yt∣<t)]=p[t≤T[Δϕ→(Yt∣<t)∣<t]]=p[t≤TΔϕ→(<t)], E _p [ 1\t≤ T\\, φ(Y_t Y_<t) ]= E _p [ 1\t≤ T\\, E [ φ(Y_t Y_<t) Y_<t] ]= E _p [ 1\t≤ T\\,B φ( Y_<t) ], by Eq. 9. Absolute summability lets us sum over t inside the expectation on both sides, giving p[∑t=1TΔϕ→(Yt∣<t)]=p[∑t=1TΔϕ→(<t)]=p[θ^(ϕ→)]−ϕ→(ε). E _p [ _t=1^T φ(Y_t Y_<t) ]= E _p [ _t=1^TB φ( Y_<t) ]= E _p[ θ( φ)]- φ( ). The left-hand side equals p[f()]−ϕ→(ε)=θ−ϕ→(ε) E _p[f( Y)]- φ( )=θ- φ( ) by the telescoping identity Eq. 7, which proves the claim. ∎ B.3 Proof of Proposition 1 Proof of Proposition 1. Oracle potential. For a prefix y with p→()>0 p( y)>0, the oracle potential is the continuation value ϕ→⋆()=p[f()∣⪰]=1p→()∑∈Σ∗p()f(), φ ( y)= E _p[f( Y) Y y]= 1 p( y) _ z∈ ^*p( y z)\,f( y z), which is finite whenever p|f()|<∞ E _p|f( Y)|<∞. In particular ϕ→⋆(ε)=θ φ ( )=θ. Splitting the continuation z according to its first symbol gives the one-step recursion ϕ→⋆()=∑a∈Σp→(a∣)ϕ→⋆(a)+p→(eos∣)f(), φ ( y)= _a∈ p(a y)\, φ ( ya)+ p( eos y)\,f( y), (16) where terms with p→(a∣)=0 p(a y)=0 are read as zero. The right-hand side is precisely ϕ→⋆()B φ ( y), so the expected increment of Eq. 9 vanishes at every prefix of positive probability: Δϕ→⋆()=ϕ→⋆()−ϕ→⋆()=0.B φ ( y)=B φ ( y)- φ ( y)=0. Every prefix <t(m) Y^(m)_<t visited by a sampled string has positive probability, so all inner sums in Eq. 10 vanish and θ^(ϕ→⋆)=ϕ→⋆(ε)+1M∑m=1M∑t=1|(m)|+10=θ θ( φ )= φ ( )+ 1M _m=1^M _t=1^| Y^(m)|+10=θ pathwise. A constant has zero variance, which proves the first claim. Variance comparison. Writing θ^(ϕ→)=θ^∗−(θ^∗−θ^(ϕ→)) θ( φ)= θ_*- ( θ_*- θ( φ) ) and expanding, Varp(θ^(ϕ→))=Varp(θ^∗)+Varp(θ^∗−θ^(ϕ→))−2Covp(θ^∗,θ^∗−θ^(ϕ→)).Var_p ( θ( φ) )=Var_p( θ_*)+Var_p ( θ_*- θ( φ) )-2Cov_p ( θ_*, θ_*- θ( φ) ). Therefore Varp(θ^(ϕ→))≤Varp(θ^∗)Var_p( θ( φ)) _p( θ_*) holds if and only if Varp(θ^∗−θ^(ϕ→))≤2Covp(θ^∗,θ^∗−θ^(ϕ→))Var_p( θ_*- θ( φ))≤ 2Cov_p( θ_*, θ_*- θ( φ)), which is the stated condition. ∎ B.4 Proof of Proposition 2 Proof of Proposition 2. Unbiasedness. By linearity, p[θ^(ϕ→,α)]=p[θ^∗]−α(p[θ^∗]−p[θ^(ϕ→)])=θ−α(θ−θ)=θ, E _p[ θ( φ,α)]= E _p[ θ_*]-α ( E _p[ θ_*]- E _p[ θ( φ)] )=θ-α(θ-θ)=θ, using that θ^∗ θ_* is unbiased by assumption and θ^(ϕ→) θ( φ) is unbiased by Proposition 4. Optimal coefficient. Write D=defθ^∗−θ^(ϕ→)D 0.5pt def= θ_*- θ( φ). Since both estimators have finite second moments, Varp(θ^(ϕ→,α))=Varp(θ^∗−αD)=Varp(θ^∗)−2αCovp(θ^∗,D)+α2Varp(D).Var_p ( θ( φ,α) )=Var_p( θ_*-α D)=Var_p( θ_*)-2α\,Cov_p( θ_*,D)+α^2\,Var_p(D). (17) When Varp(D)>0Var_p(D)>0, this is a strictly convex quadratic in α, minimized where its derivative vanishes, αVarp(D)=Covp(θ^∗,D)α\,Var_p(D)=Cov_p( θ_*,D), that is, at α⋆=Covp(θ^∗,D)/Varp(D)α =Cov_p( θ_*,D)/Var_p(D) as in Eq. 12. Substituting α⋆α into Eq. 17, Varp(θ^(ϕ→,α⋆))=Varp(θ^∗)−Covp(θ^∗,D)2Varp(D)≤Varp(θ^∗).Var_p ( θ( φ,α ) )=Var_p( θ_*)- Cov_p( θ_*,D)^2Var_p(D) _p( θ_*). When Varp(D)=0Var_p(D)=0, D is almost surely equal to its mean, which is zero by unbiasedness, so θ^(ϕ→,α)=θ^∗ θ( φ,α)= θ_* almost surely for every α. ∎ B.5 Potentials under a Proposal We now drop the assumption q=pq=p of Section 3. Let q be a proposal with respect to which p is absolutely continuous, so that q()=0q( y)=0 implies p()=0p( y)=0, and let (1),…,(M)∼i.i.d.q Y^(1),…, Y^(M) i.i.d. q. The potential construction remains unchanged: given a potential ϕ→ φ, we first form the corresponding single-sample potential estimator under the target distribution p, and then apply the usual importance weight to account for sampling from q. Definition 3 (Importance-weighted potential estimator). Let (1),…,(M)∼i.i.d.q Y^(1),…, Y^(M) i.i.d. q. Given a potential ϕ→ φ, define θ^IS(ϕ→)=def1M∑m=1Mp((m))q((m))[ϕ→(ε)+∑t=1|(m)|+1Δϕ→(<t(m))]. θ_IS( φ) 0.5pt def= 1M _m=1^M p( Y^(m))q( Y^(m)) [ φ( )+ _t=1^| Y^(m)|+1B φ( Y^(m)_<t) ]. (18) The next proposition is the counterpart of Proposition 1 under a proposal. Proposition 5 (Unbiasedness and variance under a proposal). Let ϕ→ φ satisfy the summability condition of Proposition 4, and assume that θ^IS(ϕ→) θ_IS( φ) has a finite second moment. Assume that θ^IS(ϕ→) θ_IS( φ) has a finite second moment. Then q[θ^IS(ϕ→)]=θ E _q[ θ_IS( φ)]=θ. Moreover, for the oracle potential, Varq(θ^IS(ϕ→⋆))=θ2M(∑∈Σ∗p()2q()−1).Var_q ( θ_IS( φ ) )= θ^2M ( _ y∈ ^* p( y)^2q( y)-1 ). In particular, if θ≠0θ≠ 0, the oracle potential yields zero variance if and only if q=pq=p. Proof. Unbiasedness. It suffices to treat M=1M=1. By a change of measure, q[θ^IS(ϕ→)] E _q[ θ_IS( φ)] =q[p()q()(ϕ→(ε)+∑t=1||+1Δϕ→(<t))] = E _q [ p( Y)q( Y) ( φ( )+ _t=1^| Y|+1B φ( Y_<t) ) ] =p[ϕ→(ε)+∑t=1||+1Δϕ→(<t)]=θ, = E _p [ φ( )+ _t=1^| Y|+1B φ( Y_<t) ]=θ, where the last equality follows from Proposition 4. Variance. By Proposition 1, ϕ→⋆(ε)+∑t=1||+1Δϕ→⋆(<t)=θ φ ( )+ _t=1^| Y|+1B φ ( Y_<t)=θ almost surely under p. Hence, MVarq(θ^IS(ϕ→⋆)) MVar_q ( θ_IS( φ ) ) =Varq(p()q()θ) =Var_q ( p( Y)q( Y)θ ) =θ2(q[p()2q()2]−1) =θ^2 ( E _q [ p( Y)^2q( Y)^2 ]-1 ) =θ2(∑∈Σ∗p()2q()−1). =θ^2 ( _ y∈ ^* p( y)^2q( y)-1 ). Finally, ∑∈Σ∗p()2q()≥(∑∈Σ∗p())2=1, _ y∈ ^* p( y)^2q( y)≥ ( _ y∈ ^*p( y) )^2=1, with equality if and only if q=pq=p. Thus, when θ≠0θ≠ 0, the variance is zero if and only if q=pq=p. ∎ The control-variate estimator of Definition 2 is formed with the baseline θ^∗=θ^IS θ_*= θ_IS and θ^(ϕ→) θ( φ) replaced by θ^IS(ϕ→) θ_IS( φ), θ^IS(ϕ→,α)=defθ^IS−α(θ^IS−θ^IS(ϕ→)),α∈ℝ. θ_IS( φ,α) 0.5pt def= θ_IS-α ( θ_IS- θ_IS( φ) ), α . (19) The proof of Proposition 2 uses only the unbiasedness and finite second moments of the two estimators, so it applies verbatim: θ^IS(ϕ→,α) θ_IS( φ,α) is unbiased for every α, and at α⋆α its variance is at most that of θ^IS θ_IS. Appendix C Additional Analysis and Experimental Details This appendix collects the setup behind the experiments of Section 4: the models and prompt banks, the target-construction protocols, the generation configurations, and the supporting figures and tables. One subsection per application, in the order of Section 4. C.1 Analysis of Watermark Detection A language model watermark embeds a statistical signal in generated text while preserving the model’s intended task. Watermark detection then aims to reliably assess whether a text of a given length was generated by a watermarked model. We use the scheme of Kirchenbauer et al. (2023): at each generation step, a pseudorandom green list G()⊆ΣG( [rgb]0,0,0 x) is determined from the prompt and generated prefix [rgb]0,0,0 x, and the logits of tokens in this set are increased by a bias δ. Let γ denote the fraction of the vocabulary assigned to the green list, m the generation length, and σ=defmγ(1−γ)σ 0.5pt def= mγ(1-γ). Detection then counts the number of generated tokens that fall in their respective green lists and standardizes this count, ||G−γmσ,||G=def#t:Yt∈G(<t), | Y|_G-γ mσ, | Y|_G 0.5pt def=\#\t:Y_t∈ G( Y_<t)\, (20) declaring the string watermarked when this exceeds a prescribed threshold τ, equivalently when ||G≥γm+τσ| Y|_G≥γ m+τσ. C.1.1 Estimands and potential definitions Table C.1.1: Expected watermark score Table C.1.1: Probability of watermark detection f()=(||G−γm)/σf( Y)= (| Y|_G-γ m )/σ f()=||G≥γm+τσf( Y)= 1\| Y|_G≥γ m+τσ\ ||G=||G+||G| [rgb]0,0,0 x [rgb]0,0,0 z|_G=| [rgb]0,0,0 x|_G+| [rgb]0,0,0 z|_G not additive in the prefix ϕ→⋆()=(||G+p[||G∣⪰]−γm)/σ φ ( [rgb]0,0,0 x)= (| [rgb]0,0,0 x|_G+ E _p[| [rgb]0,0,0 Z|_G [rgb]0,0,0 Z [rgb]0,0,0 x]-γ m )/σ ϕ→⋆()=Prp[||G≥γm+τσ∣⪰] φ ( [rgb]0,0,0 x)=Pr_p [| Y|_G≥γ m+τσ Y [rgb]0,0,0 x ] ϕ→()=(||G−γm)/σ φ( [rgb]0,0,0 x)= (| [rgb]0,0,0 x|_G-γ m )/σ ϕ→(,q)=Pr[Bin(m−||,q)≥k()] φ( [rgb]0,0,0 x,q)=Pr [Bin(m-| [rgb]0,0,0 x|,q)≥ k( [rgb]0,0,0 x) ] Δϕ→(a∣)=a∈G()/σ φ(a [rgb]0,0,0 x)= 1\a∈ G( [rgb]0,0,0 x)\/σ Δϕ→(a∣,q)=ϕ→(a,q)−ϕ→(,q) φ(a [rgb]0,0,0 x,q)= φ( [rgb]0,0,0 xa,q)- φ( [rgb]0,0,0 x,q) Δϕ→(eos∣)=0 φ( eos [rgb]0,0,0 x)=0 Δϕ→(eos∣,q)=0 φ( eos [rgb]0,0,0 x,q)=0 At a fixed prompt and generation length, two expectations characterize the detector. The expected score p[(||G−γm)/σ] E _p [(| Y|_G-γ m)/σ ] measures how strongly the watermark registers on average; the detection probability Prp[||G≥γm+τσ]Pr_p[| Y|_G≥γ m+τσ] is the detector’s true-positive rate, and is the quantity a deployment actually cares about. The two behave differently under our framework. The score has an additive structure over prefixes (Section C.1.1) so the value of the function at the observed prefix is itself a potential, exactly as for the losses of Section 4.1. The detection indicator is not additive in the prefix, so it admits no such decomposition and must instead be approached by approximating the oracle potential, as in the event analysis of Section 4.3. Here that oracle potential is the conditional probability that the detector will eventually fire. Detection probability: analytical approximation of the oracle potential. For the detection probability a natural analytical approximation of the oracle potential is available. Write k()=γm+τσ−||Gk( [rgb]0,0,0 x)=γ m+τσ-| [rgb]0,0,0 x|_G for the number of green tokens still needed at the prefix [rgb]0,0,0 x. If the remaining m−||m-| [rgb]0,0,0 x| tokens were green independently with a common rate q, the detector would fire with probability Pr[Bin(m−||,q)≥k()]Pr[Bin(m-| [rgb]0,0,0 x|,q)≥ k( [rgb]0,0,0 x)], which is the potential of Section C.1.1. It depends on the single parameter q, which we fit globally for each length m by cross-fitting on 1010 folds of the continuations: the q applied to one fold is the average green-list mass of the other nine, pooled over prompts. Computing this potential requires no additional model evaluations, as it uses only the green-list mass ∑a∈G()p→(a∣) _a∈ G( [rgb]0,0,0 x) p(a [rgb]0,0,0 x) from the already-computed next-token distribution. C.1.2 Experiments and results δ m median [p10,p90][p_10,\,p_90] worst prompt 11 3232 0.4710.471 [0.413, 0.564][0.413,\,0.564] 0.9070.907 6464 0.5130.513 [0.460, 0.597][0.460,\,0.597] 0.8940.894 128128 0.5610.561 [0.508, 0.677][0.508,\,0.677] 0.9310.931 200200 0.6120.612 [0.544, 0.748][0.544,\,0.748] 0.9650.965 22 3232 0.5380.538 [0.471, 0.631][0.471,\,0.631] 0.8360.836 6464 0.5840.584 [0.520, 0.671][0.520,\,0.671] 0.8760.876 128128 0.6440.644 [0.583, 0.748][0.583,\,0.748] 0.9250.925 200200 0.6990.699 [0.629, 0.818][0.629,\,0.818] 0.9620.962 55 3232 0.7220.722 [0.608, 0.807][0.608,\,0.807] 0.9220.922 6464 0.7430.743 [0.665, 0.817][0.665,\,0.817] 0.9600.960 128128 0.7900.790 [0.725, 0.867][0.725,\,0.867] 0.9930.993 200200 0.8370.837 [0.776, 0.910][0.776,\,0.910] 0.9960.996 Table 6: Per-prompt variance of θ^(ϕ→) θ( φ) against Monte Carlo, 512512 prompts. Setup. We use the same setup of Kirchenbauer et al. (2023), with green-list fraction γ=0.5γ=0.5 and logit biases δ∈1,2,5δ∈\1,2,5\. For each configuration, we generate 4,0004,000 continuations per prompt from OPT-1.3B for 512512 prompts from C4 RealNewsLike, using temperature 0.70.7. We evaluate continuation lengths m∈32,64,128,200m∈\32,64,128,200\ and detector thresholds z0∈3,4,5z_0∈\3,4,5\. (a) binomial potential and green-count paths (b) variance ratio across regimes Figure 3: Panel (a) shows how sampled trajectories update the green-token count over 64 generation steps. Red trajectories are those whose final green-token count exceeds the detection threshold, while black trajectories remain below it. The background encodes the potential’s assigned probability of eventual detection at each step and green-count state, with white indicating probability zero. Red trajectories tend to pass through higher-probability regions, whereas black trajectories concentrate in low- or zero-probability regions. Panel (b) shows the variance reduction at different detection probabilities For each of the two estimands, we compare the associated potentials against Monte Carlo. For the watermark-detection probability, we also include the baseline terminal potential, i.e., the corresponding observed-prefix potential for non-additive functions: zero at every nonterminal prefix and equal to the function value at the terminal prefix. Detection probability: the binomial potential removes about 44%44\% of the Monte Carlo variance when detection is uncertain, but little once it saturates. Figure 3(b) shows the variance ratio between the potential estimator and MC for each prompt and configuration (δ,m,z0)(δ,m,z_0). Across the 14,36614,366 pairs for which detection is not deterministic, the median ratio is 0.640.64, corresponding to a 36%36\% reduction in variance. Since variance scales inversely with the number of samples, this is equivalent to reaching the precision of 4,0004,000 Monte Carlo samples with about 2,5602,560 samples, without additional model evaluations. The gain is largest when detection is uncertain: for 0.001≤π^<0.50.001≤ π<0.5, the median ratio is 0.560.56, corresponding to a 44%44\% reduction. It decreases as detection saturates, reaching a median ratio of 0.970.97 for π^≥0.99 π≥ 0.99. Overall, 92.7%92.7\% of pairs improve over Monte Carlo, with similar behavior across δ, m, and z0z_0. The terminal potential, by contrast, remains close to Monte Carlo, with a median ratio of 0.9730.973. Figure 3(a) gives intuition for the variance reduction. From early in generation, trajectories that are eventually detected overlap with regions of positive potential, whereas undetected trajectories lie mostly in the white region where the potential is zero. Thus, the potential already contains information about the final detection outcome long before the trajectory ends. Expected score: every prompt improves, by 1616–53%53\% depending on δ and m. Table 6 gives the per-prompt variance relative to Monte Carlo over the 512512 prompts. The median ratio rises from 0.4710.471 at δ=1,m=32δ=1,m=32 to 0.8370.837 at δ=5,m=200δ=5,m=200: the marginalization helps most when the watermark is weak and the continuation is short, the regime in which a single green indicator carries the most conditional variance. The spread over prompts is narrow, the p10p_10–p90p_90 width never exceeding 0.210.21, and every one of the 512512 prompts lies below Monte Carlo in all twelve cells, the worst single prompt reaching 0.9960.996. We also cross-fit a control-variate coefficient on top of but observe that it buys almost nothing here: across the twelve cells the variance against θ^(ϕ→) θ( φ) stays within 0.9910.991–0.9980.998 and the median prompt within 0.9940.994–1.0001.000. C.2 Losses and Divergences: Details and Further Results C.2.1 Potential definitions Below we give the potential tables for the losses and divergences of Section 4.1; the corresponding gradients are in Section C.3. Continuation sums include termination, with ||+1=||+1=eos Y_| Y|+1= [rgb]0,0,0 z_| [rgb]0,0,0 z|+1= eos. Table 7: KL divergence f()=logp()q()f( Y)= p( Y)q( Y) Test function ϕ→⋆()=logp→()q→()+p[logp()q→()q()p→()|⪰] φ ( [rgb]0,0,0 x)= p( [rgb]0,0,0 x) q( [rgb]0,0,0 x)+ E _p\! [ p( Y)\, q( [rgb]0,0,0 x)q( Y)\, p( [rgb]0,0,0 x)\; |\; Y [rgb]0,0,0 x ] oracle potential ϕ→()=logp→()q→() φ( [rgb]0,0,0 x)= p( [rgb]0,0,0 x) q( [rgb]0,0,0 x) Potential Δϕ→(a∣)=logp→(a∣)q→(a∣) φ(a [rgb]0,0,0 x)= p(a [rgb]0,0,0 x) q(a [rgb]0,0,0 x) Increment (a∈Σa∈ ) Δϕ→(eos∣)=logp→(eos∣)q→(eos∣) φ( eos [rgb]0,0,0 x)= p( eos [rgb]0,0,0 x) q( eos [rgb]0,0,0 x) eos increment Δϕ→()=DKL(p→(⋅∣)∥q→(⋅∣))B φ( [rgb]0,0,0 x)=D_KL ( p(· [rgb]0,0,0 x)\,\|\, q(· [rgb]0,0,0 x) ) Expected increment Table 8: Brier loss f()=∑t∑b∈Σ∪eos(p→(b∣<t)−b=Yt)2f( Y)= _t _b∈ ∪\ eos\ ( p(b Y_<t)- 1\b=Y_t\ )^2 Δϕ→(a∣)=∑b∈Σ∪eos(p→(b∣)−b=a)2=1−2p→(a∣)+∑b∈Σ∪eosp→(b∣)2 φ(a [rgb]0,0,0 x)= _b∈ ∪\ eos\ ( p(b [rgb]0,0,0 x)- 1\b=a\ )^2=1-2\, p(a [rgb]0,0,0 x)+ _b∈ ∪\ eos\ p(b [rgb]0,0,0 x)^2 f()=ϕ→()+∑t=1||+1Δϕ→(t∣<t)f( [rgb]0,0,0 x [rgb]0,0,0 z)= φ( [rgb]0,0,0 x)+ _t=1^| [rgb]0,0,0 z|+1 φ( [rgb]0,0,0 z_t [rgb]0,0,0 x [rgb]0,0,0 z_<t) ϕ→⋆()=ϕ→()+p[∑t=||+1||+1Δϕ→(t∣<t)|⪰] φ ( [rgb]0,0,0 x)= φ( [rgb]0,0,0 x)+ E _p\! [ _t=| [rgb]0,0,0 x|+1^| Y|+1 φ( Y_t Y_<t)\; |\; Y [rgb]0,0,0 x ] ϕ→()=∑t≤||Δϕ→(t∣<t) φ( [rgb]0,0,0 x)= _t≤| [rgb]0,0,0 x| φ( [rgb]0,0,0 x_t [rgb]0,0,0 x_<t) Δϕ→(eos∣)=1−2p→(eos∣)+∑b∈Σ∪eosp→(b∣)2 φ( eos [rgb]0,0,0 x)=1-2\, p( eos [rgb]0,0,0 x)+ _b∈ ∪\ eos\ p(b [rgb]0,0,0 x)^2 Δϕ→()=1−∑b∈Σ∪eosp→(b∣)2B φ( [rgb]0,0,0 x)=1- _b∈ ∪\ eos\ p(b [rgb]0,0,0 x)^2 Table 9: Power-score regret, exponent β f()=∑t[β(p→(Yt∣<t)β−1−q→(Yt∣<t)β−1)−(β−1)∑b∈Σ∪eos(p→(b∣<t)β−q→(b∣<t)β)]f( Y)= _t [β ( p(Y_t Y_<t)^β-1- q(Y_t Y_<t)^β-1 )-(β-1) _b∈ ∪\ eos\ ( p(b Y_<t)^β- q(b Y_<t)^β ) ] Δϕ→(a∣)=β(p→(a∣)β−1−q→(a∣)β−1)−(β−1)∑b∈Σ∪eos(p→(b∣)β−q→(b∣)β) φ(a [rgb]0,0,0 x)=β ( p(a [rgb]0,0,0 x)^β-1- q(a [rgb]0,0,0 x)^β-1 )-(β-1) _b∈ ∪\ eos\ ( p(b [rgb]0,0,0 x)^β- q(b [rgb]0,0,0 x)^β ) f()=ϕ→()+∑t=1||+1Δϕ→(t∣<t)f( [rgb]0,0,0 x [rgb]0,0,0 z)= φ( [rgb]0,0,0 x)+ _t=1^| [rgb]0,0,0 z|+1 φ( [rgb]0,0,0 z_t [rgb]0,0,0 x [rgb]0,0,0 z_<t) ϕ→⋆()=ϕ→()+p[∑t=||+1||+1Δϕ→(t∣<t)|⪰] φ ( [rgb]0,0,0 x)= φ( [rgb]0,0,0 x)+ E _p\! [ _t=| [rgb]0,0,0 x|+1^| Y|+1 φ( Y_t Y_<t)\; |\; Y [rgb]0,0,0 x ] ϕ→()=∑t≤||Δϕ→(t∣<t) φ( [rgb]0,0,0 x)= _t≤| [rgb]0,0,0 x| φ( [rgb]0,0,0 x_t [rgb]0,0,0 x_<t) Δϕ→(eos∣)=β(p→(eos∣)β−1−q→(eos∣)β−1)−(β−1)∑b∈Σ∪eos(p→(b∣)β−q→(b∣)β) φ( eos [rgb]0,0,0 x)=β ( p( eos [rgb]0,0,0 x)^β-1- q( eos [rgb]0,0,0 x)^β-1 )-(β-1) _b∈ ∪\ eos\ ( p(b [rgb]0,0,0 x)^β- q(b [rgb]0,0,0 x)^β ) Δϕ→()=∑b∈Σ∪eos[p→(b∣)β−βp→(b∣)q→(b∣)β−1+(β−1)q→(b∣)β]B φ( [rgb]0,0,0 x)= _b∈ ∪\ eos\ [ p(b [rgb]0,0,0 x)^β-β\, p(b [rgb]0,0,0 x)\, q(b [rgb]0,0,0 x)^β-1+(β-1)\, q(b [rgb]0,0,0 x)^β ] Entropy is the same construction with Δϕ→(a∣)=−logp→(a∣) φ(a [rgb]0,0,0 x)=- p(a [rgb]0,0,0 x), shown in Section 4.1. The four power-score regrets use β∈1.25,1.5,1.75,2β∈\1.25,1.5,1.75,2\; as β→1β→ 1 the per-token term, divided by β−1β-1, tends to the KL term of Section C.2.1. C.2.2 Experimental setup We evaluate the losses and divergences on eight configurations, each defined by a sampling policy, a reference model, and a generation horizon (10). These cover seven distinct policy–reference pairs: the GPT-2 sentiment pair is evaluated at two horizons, one of 5050 tokens that reuses the setting of Amini et al. (2025) for reproducibility, and one of 256256 tokens to match the longer window used for all other models. Model repositories are listed in 11. Shared sampling protocol. Each configuration uses 512512 prompts and B responses per prompt, with B specified in 10. Responses are drawn by ancestral sampling at temperature 11, without top-k or top-p filtering, using fp32 precision. Generation stops at eos or the specified horizon. Table 10: Policy–reference banks of the proper-loss experiment. Labels are those used in the figures; p is the sampling policy and q the reference; model names are resolved in 11. All banks have 512512 prompts and B responses per prompt. bank prompts sampling policy p reference q horizon B IMDb-256 IMDb prefixes gpt2-imdb + DPO gpt2-imdb 256 4000 IMDb-50 IMDb prefixes same policy gpt2-imdb 50 4000 Tülu-3 AlpacaEval 2 Tulu-3-8B-DPO Tulu-3-8B-SFT 256 100 Gemma It3 vs It0 AlpacaEval 2 SPPO-Iter3 gemma-2-9b-it 256 100 Gemma It1 vs It0 AlpacaEval 2 SPPO-Iter1 gemma-2-9b-it 256 100 Gemma It2 vs It0 AlpacaEval 2 SPPO-Iter2 gemma-2-9b-it 256 100 Gemma It2 vs It1 AlpacaEval 2 SPPO-Iter2 SPPO-Iter1 256 100 Gemma It3 vs It2 AlpacaEval 2 SPPO-Iter3 SPPO-Iter2 256 100 GPT-2-imbd + DPO: policy training. We train the policy p ourself, following the setting of Amini et al. (2025). The reference q is gpt2-imdb; we obtain the policy p by fine-tuning it for positive sentiment with direct preference optimization (DPO) on 5,0005,000 IMDb training reviews. Training prompts are 22–88-token prefixes. Four sampled responses per prompt form six response pairs, with preferences determined by distilbert-imdb. We train for one epoch (938938 steps), using βDPO=0.1 _DPO=0.1 and a learning rate of 3×10−63× 10^-6. GPT-2 sentiment: evaluation prompts. We select 512512 reviews from the IMDb test split using a fixed seed. For each review, we draw a prefix length uniformly from 2,…,8\2,…,8\ and retain the corresponding GPT-2 token IDs as the prompt. We draw 4,0004,000 responses per prompt at a primary horizon of 256256 tokens and a control horizon of 5050 tokens, using the same policy and reference in both configurations. Chat models: policy–reference comparisons. For Tülu-3, we compare the published DPO policy with its SFT reference. For Gemma-2, we use the three published SPPO iterations, denoting the initial instruction model by iteration 00. Comparing iterations 11, 22, and 33 with iteration 00 measures accumulated drift. We additionally compare iteration 22 with 11 and iteration 33 with 22 to examine local drift between successive updates. Chat models: evaluation prompts. The chat-model configurations use 512512 instructions from the 805805-instruction AlpacaEval 2 set. Instructions are rendered with the corresponding model’s chat template; within each policy–reference pair, both models receive identical token IDs. Each prompt receives 100100 responses, with a horizon of 256256 tokens. Table 11: Model repositories (Hugging Face) for the names used in 10. model repository gpt2-imdb lvwerra/gpt2-imdb gpt2-imdb + DPO trained by us from lvwerra/gpt2-imdb (recipe above) distilbert-imdb (DPO preference judge) lvwerra/distilbert-imdb Tulu-3-8B-SFT allenai/Llama-3.1-Tulu-3-8B-SFT Tulu-3-8B-DPO allenai/Llama-3.1-Tulu-3-8B-DPO gemma-2-9b-it (SPPO Iter0) google/gemma-2-9b-it SPPO-Iter1 UCLA-AGI/Gemma-2-9B-It-SPPO-Iter1 SPPO-Iter2 UCLA-AGI/Gemma-2-9B-It-SPPO-Iter2 SPPO-Iter3 UCLA-AGI/Gemma-2-9B-It-SPPO-Iter3 C.2.3 Additional results Variance reduction across quantities. Figure 4 extends the magnitude–variance comparison of Figure 1(b) to all seven quantities. For each quantity, it relates the variance ratio to the estimated magnitude of the quantity, showing both configuration-level aggregates and individual prompts. Figure 4: Variance reduction as a function of quantity magnitude for all losses and divergences in the proper-loss experiment. Each panel shows the variance ratio Varp(θ^(ϕ→))/Varp(θ^MC)Var_p( θ( φ))/Var_p( θ_MC) against the estimated magnitude of the quantity, with one colored marker per policy–reference bank and grey points for individual prompts. For KL and the power-score regrets, variance reduction becomes stronger as the policy–reference divergence decreases, with this dependence becoming more pronounced for larger β. Entropy and Brier loss instead remain comparatively close to the Monte Carlo baseline across banks. Prompt-level behavior and control variates. 12 summarizes the per-prompt variance ratios of the potential estimator and its cross-fitted control-variate variant relative to Monte Carlo. It also reports separately the prompts on which the potential alone increases variance, showing how fitting the control-variate coefficient changes their ratios and how often it brings them below one. θ^(ϕ→) θ( φ) vs. θ^MC θ_MC θ^(ϕ→,α⋆) θ( φ,α ) vs. θ^MC θ_MC prompts with Varp(θ^(ϕ→))>Varp(θ^MC)Var_p( θ( φ))>Var_p( θ_MC) quantity median [IQR] % <1<1 median [IQR] % <1<1 n med. θ^(ϕ→) θ( φ) med. θ^(ϕ→,α⋆) θ( φ,α ) [95% CI] % <1<1 DKLD_KL 0.204 [0.061, 0.503] 98.5 0.204 [0.061, 0.482] 99.5 63 1.134 0.974 [0.945, 0.990] 65.1 regret β=1.25β=1.25 0.202 [0.063, 0.441] 99.0 0.193 [0.063, 0.436] 99.8 40 1.083 0.951 [0.932, 0.973] 77.5 regret β=1.5β=1.5 0.164 [0.063, 0.400] 99.6 0.160 [0.063, 0.396] 99.9 15 1.058 0.935 [0.891, 1.044] 66.7 regret β=1.75β=1.75 0.142 [0.062, 0.365] 99.9 0.141 [0.061, 0.363] 99.9 4 1.190 1.080 [0.925, 1.102] 25.0 regret β=2β=2 0.131 [0.060, 0.330] 99.9 0.129 [0.060, 0.329] 99.9 3 1.125 1.040 [1.009, 1.043] 0.0 entropy 0.633 [0.509, 0.796] 99.8 0.630 [0.498, 0.792] 99.5 10 1.047 1.090 [1.050, 1.096] 0.0 Brier loss 0.694 [0.545, 0.903] 99.8 0.698 [0.543, 0.903] 99.5 9 1.040 1.059 [1.024, 1.080] 0.0 Table 12: Per-prompt variance ratios for the losses and divergences experiments, pooling all eight model pairs. We compare the potential estimator θ^(ϕ→) θ( φ) and its cross-fitted control-variate variant θ^(ϕ→,α⋆) θ( φ,α ) against Monte Carlo. The rightmost columns restrict attention to prompts for which θ^(ϕ→) θ( φ) has higher variance than Monte Carlo and show how often the control variate reduces the ratio below one. C.2.4 Reproducing Amini et al. (2025) We first reproduce the KL-estimation experiment of Amini et al. (2025) using our sentiment-tuned policy and compare against the estimators considered in their paper. Their term-wise Rao–Blackwellized estimator is recovered as θ^(ϕ→) θ( φ), corresponding to a unit coefficient on the associated zero-mean correction. We reproduce its variance reduction over Monte Carlo (13). Fitting the control-variate coefficient reduces variance further: at M=1M=1, the standard deviation decreases from 2.052.05 to 1.901.90. Thus, the unit coefficient is not in general variance-optimal. std. dev. estimator mean M=1M=1 M=5M=5 M=10M=10 θ^MC θ_MC 6.69 3.32 1.49 1.05 HT 6.69 3.32 1.49 1.05 Schulman CV, α=1α=1 6.47 35.67 15.95 11.28 Schulman CV, fitted α 6.36 11.44 5.11 3.61 θ^(ϕ→) θ( φ) 6.69 2.05 0.92 0.65 θ^(ϕ→,α⋆) θ( φ,α ) 6.69 1.901.90 0.850.85 0.600.60 Table 13: Reproduction of the KL-estimation experiment of Amini et al. (2025): mean estimate and empirical standard deviation (nats) at sample sizes M=1,5,10M=1,5,10, over 512512 prompts with 4,0004,000 samples each. The mean does not depend on M. Lowest standard deviation per column in bold. Counterexample to the variance guarantee of Amini et al. (2025). Theorem 2 of Amini et al. (2025) claims that their term-wise Rao–Blackwellized estimator never has larger variance than Monte Carlo. Consider language models p and q supported on length-two strings, with p(ax)=p(bx)=12,q(ax)=q(bx)=14,q(bz)=12,p(ax)=p(bx)= 12, q(ax)=q(bx)= 14, q(bz)= 12, and zero probability elsewhere. Under p, logp(Y)q(Y)=log2 p(Y)q(Y)= 2 for both possible strings, so Varp(θ^MC)=0Var_p( θ_MC)=0. Their estimator, corresponding to θ^(ϕ→) θ( φ) in our notation, instead equals 12log43 12 43 or 12log43+log3 12 43+ 3 with equal probability. It remains unbiased, but Varp(θ^(ϕ→))=14(log3)2>0=Varp(θ^MC).Var_p\! ( θ( φ) )= 14( 3)^2>0=Var_p( θ_MC). This also follows directly from Proposition 1: since θ^MC θ_MC is constant, Covp(θ^MC,θ^MC−θ^(ϕ→))=0<12Varp(θ^MC−θ^(ϕ→))Cov_p( θ_MC, θ_MC- θ( φ))=0< 12Var_p( θ_MC- θ( φ)), so the condition for variance reduction fails. Accordingly, Proposition 2 gives α⋆=0α =0, whereas the estimator of Amini et al. (2025) corresponds to α=1α=1. C.3 Gradients of Losses and Divergences C.3.1 Potential definitions Let ∈ℝdw ^d denote the parameters of p, and write ∇=def∇ 0.5pt def= _w, holding q fixed. The original test function f is scalar-valued, whereas its gradient with respect to w is d-dimensional. This gives the vector-valued test function ∇()=def∇f()+f()∇logp()∈ℝd. f_∇( Y) 0.5pt def=∇ f( Y)+f( Y)\,∇ p( Y) ^d. (21) The first term differentiates the function value for a fixed string, while the second accounts for the change in its sampling probability. Under standard regularity conditions, ∇p[f()]=p[∇()].∇\, E _p[f( Y)]= E _p[ f_∇( Y)]. (22) Gradient estimation therefore fits our framework by applying the potential estimator to ∇ f_∇ coordinatewise. Given a scalar potential ϕ→ φ for f, we define the corresponding vector-valued gradient potential ϕ→∇()=def∇ϕ→()+ϕ→()∇logp→()∈ℝd. φ_∇( [rgb]0,0,0 x) 0.5pt def=∇ φ( [rgb]0,0,0 x)+ φ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x) ^d. (23) Replacing ϕ→ φ by ϕ→⋆ φ gives the corresponding oracle gradient potential. Because the next-token probabilities in B also depend on w, the expected gradient increment satisfies Δϕ→∇()=∇Δϕ→()+Δϕ→()∇logp→().B φ_∇( [rgb]0,0,0 x)= φ( [rgb]0,0,0 x)+B φ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x). (24) We extend the operators Δ and B to vector-valued potentials coordinatewise. Section C.3.1 gives the resulting constructions for the losses in Sections 4.1, C.2.1, C.2.1 and C.2.1. For KL, Δϕ→()B φ( [rgb]0,0,0 x) is the KL divergence between the next-token distributions, recovering the local gradient expression of Amini et al. (2025). Table 14: Gradient of a loss, ∇p[f()]∇\, E _p[f( Y)] ∇()=∇f()+f()∇logp() f_∇( Y)=∇ f( Y)+f( Y)\,∇ p( Y) Test function Δϕ→∇(a∣)=∇Δϕ→(a∣)+Δϕ→(a∣)∇logp→(a)+ϕ→()∇logp→(a∣) aligned φ_∇(a [rgb]0,0,0 x)&=∇ φ(a [rgb]0,0,0 x)+ φ(a [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 xa)\\ & + φ( [rgb]0,0,0 x)\,∇ p(a [rgb]0,0,0 x) aligned Increment (a∈Σa∈ ) ∇()=∇ϕ→()+ϕ→()∇logp→()+∇[f()−ϕ→()]+[f()−ϕ→()]∇logp()+ϕ→()∇logp→(∣) aligned f_∇( [rgb]0,0,0 x [rgb]0,0,0 z)&=∇ φ( [rgb]0,0,0 x)+ φ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x)\\ & +∇ [f( [rgb]0,0,0 x [rgb]0,0,0 z)- φ( [rgb]0,0,0 x) ]+ [f( [rgb]0,0,0 x [rgb]0,0,0 z)- φ( [rgb]0,0,0 x) ]\,∇ p( [rgb]0,0,0 x [rgb]0,0,0 z)\\ & + φ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 z [rgb]0,0,0 x) aligned Function structure ϕ→∇⋆()=∇ϕ→⋆()+ϕ→⋆()∇logp→() φ _∇( [rgb]0,0,0 x)=∇ φ ( [rgb]0,0,0 x)+ φ ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x) oracle potential ϕ→∇()=∇ϕ→()+ϕ→()∇logp→() φ_∇( [rgb]0,0,0 x)=∇ φ( [rgb]0,0,0 x)+ φ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x) Potential Δϕ→∇(eos∣)=∇Δϕ→(eos∣)+Δϕ→(eos∣)∇logp()+ϕ→()∇logp→(eos∣) φ_∇( eos [rgb]0,0,0 x)=∇ φ( eos [rgb]0,0,0 x)+ φ( eos [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x)+ φ( [rgb]0,0,0 x)\,∇ p( eos [rgb]0,0,0 x) eos increment Δϕ→∇()=∇Δϕ→()+Δϕ→()∇logp→()B φ_∇( [rgb]0,0,0 x)= φ( [rgb]0,0,0 x)+B φ( [rgb]0,0,0 x)\,∇ p( [rgb]0,0,0 x) Expected increment C.3.2 Control-variate coefficient for vector-valued functions For the gradient we need to fit a global control-variate coefficient when the function value lies in ℝdR^d. We describe the construction for a general vector-valued test function :Σ∗→ℝd f ^* ^d with a potential ϕ→:Σ∗→ℝd φ ^* ^d, both read coordinatewise. Let ^MC,^(ϕ→)∈ℝd θ_MC,\, θ( φ) ^d denote the Monte Carlo and potential estimators for one trajectory, and define the zero-mean correction =def^MC−^(ϕ→). D 0.5pt def= θ_MC- θ( φ). (25) As in Proposition 2, we use a single scalar coefficient, here shared across all d coordinates: ^(ϕ→,α)=def^MC−α. θ( φ,α) 0.5pt def= θ_MC-α\, D. (26) Thus, α=0α=0 recovers Monte Carlo and α=1α=1 the potential estimator. We choose α to minimize the total variance trVarp(^(ϕ→,α)) _p ( θ( φ,α) ), giving α⋆=trCovp(^MC,)trVarp()=p[⟨^MC−p[^MC],−p[]⟩]p[‖−p[]‖22].α = _p ( θ_MC, D ) _p( D)= E _p\! [ θ_MC- E _p[ θ_MC],\; D- E _p[ D] ] E _p\! [\| D- E _p[ D]\|_2^2 ]. (27) Given sampled trajectories indexed by prompt i and response j, we estimate this globally as α^=∑i,j⟨^MC(ij)−^¯MC,(ij)−¯⟩∑i,j‖(ij)−¯‖22, α= _i,j θ_MC^(ij)- θ_MC,\; D^(ij)- D _i,j\| D^(ij)- D\|_2^2, (28) where the bars denote means over all trajectories. The coefficient is therefore shared across prompts and coordinates. For the gradient experiments we take =∇ f= f_∇ and ϕ→=ϕ→∇ φ= φ_∇ from Section C.3.1, with d the number of policy parameters. C.3.3 Experiments and results Marginalization reduces gradient variance across all 135 defined quantity–checkpoint comparisons. We repeat the comparison for the gradient of each quantity with respect to the policy parameters along a DPO training trajectory of the GPT-2 sentiment model (Figure 5). We evaluate 2020 checkpoints using 512512 prompts and 2020 responses per prompt, and measure the variance of a gradient estimator θ by trVarp(^) _p( θ). Across all 135135 quantity–checkpoint comparisons, the potential estimator ^(ϕ→∇) θ( φ_∇) reduces variance relative to the centered Monte Carlo gradient ^MC θ_MC, with variance ratios trVarp(^(ϕ→∇))/trVarp(^MC) _p( θ( φ_∇))/ _p( θ_MC) between 0.110.11 and 0.980.98. The gains are largest early in training. The gains are largest early in training and generally decrease as training progresses. For example, the KL gradient retains 35%35\% of the Monte Carlo variance at step 5050 and 92%92\% at step 938938, while the β=2β=2 regret retains 11%11\% and 65%65\%, respectively. Among the power-score regrets, larger β consistently yields larger reductions at every checkpoint. Fitting α⋆α further improves gradients for entropy and Brier score. For some gradients, the control-variate estimator ^(ϕ→∇,α⋆) θ( φ_∇,α ) provides an additional gain over ^(ϕ→∇) θ( φ_∇). The unit-coefficient estimator ^(ϕ→∇) θ( φ_∇) alone barely improves the entropy and Brier gradients, with variance ratios of 0.970.97–0.980.98. Their fitted coefficients are substantially larger, α⋆≈5α ≈ 5 for entropy and α⋆≈11α ≈ 11 for Brier, reducing the ratios to 0.910.91 and 0.860.86. For the power-score regrets, by contrast, α⋆α remains close to one (0.660.66–0.990.99), and fitting it provides at most about a further 3%3\% reduction. (a) gradient variance ratio (b) cross-fitted control variate (c) fitted coefficient α⋆α Figure 5: Gradient estimation along a DPO trajectory of the GPT-2 sentiment policy (checkpoints every 5050 steps; the variance of an estimator θ is trVarp(^) _p( θ) over the 512×20512× 20 responses). (a) Variance of ^(ϕ→∇) θ( φ_∇) relative to the centered Monte Carlo gradient ^MC θ_MC, one series per quantity (regrets darker for larger β; undefined at step 00). (b) Variance of the cross-fitted control-variate estimator ^(ϕ→∇,α⋆) θ( φ_∇,α ) relative to ^(ϕ→∇) θ( φ_∇). (c) The cross-fitted coefficient, mean over 88 prompt folds with the range as band; α⋆=1α =1 recovers ^(ϕ→∇) θ( φ_∇). Bands in (a) and (b) are 95%95\% prompt-cluster bootstrap intervals; one training run and one seed. C.4 n-gram Counts and Frequencies C.4.1 Definition of frequency estimator For n-gram statistics one is often interested in the frequency of a target rather than solely its expected count. Let f()=#occurrences of w in f( Y)=\#\occurrences of w in Y\ be the count, the test function of Section 4.2, and Nn()N_n( Y) the number of word n-gram slots in Y, i.e., the number of positions in the decoded text of Y at which n consecutive words occur (for n=1n=1, the number of words). Here a word is a maximal run of ASCII letters, and consecutive words must be separated by a single space, so punctuation ends an n-gram. The frequency of a substring w made of n such words is θFR=def[f()]/[Nn()] _FR 0.5pt def= E [f( Y)]/ E [N_n( Y)]. Since every slot holds exactly one n-gram, the counts of all word n-grams of order n sum to Nn()N_n( Y), so, as a function of w, θFR _FR is the marginal distribution over word n-grams of order n. We estimate it with the ratio estimators θ^FR=defθ^MC1M∑m=1MNn((m)),θ^FR(ϕ→)=defθ^(ϕ→)1M∑m=1MNn((m)), θ_FR 0.5pt def= θ_MC 1M _m=1^MN_n( Y^(m)), θ_FR( φ) 0.5pt def= θ( φ) 1M _m=1^MN_n( Y^(m)), (29) where θ^MC θ_MC and θ^(ϕ→) θ( φ) are the Monte Carlo and potential estimators of [f()] E [f( Y)] built from Section 4.2, and the denominator is the Monte Carlo estimate of [Nn()] E [N_n( Y)] over the same M generations. Ratio estimators are biased at finite sample size, but the bias vanishes asymptotically under standard regularity conditions (Cochran, 1977). C.4.2 Definition of cost-adjusted variance ratio For each target w, we compare the potential estimator against Monte Carlo at equal sample size and equal compute. For either the expected count or frequency, we define ρ^(w)=Var^(potential estimator at M=4,000)1000Var^(Monte Carlo estimator at M=4×106),ρ^c(w)=(Lw+1)ρ^(w). ρ(w)= Var (potential estimator at M=4,000 )1000\, Var (Monte Carlo estimator at M=4× 10^6 ), ρ_c(w)=(L_w+1)\, ρ(w). (30) The factor 10001000 rescales the Monte Carlo variance from 4×1064× 10^6 to 4,0004,000 samples. This rescaling is exact for expected counts, whose estimator is a sample mean. For frequencies, which are ratios of sample means, the same 1/M1/M variance scaling holds asymptotically by the multivariate delta method (Rieder, 2012, §1.3). The per-target token-evaluation cost proxy is cw=1+(Lw−1)+1=Lw+1c_w=1+(L_w-1)+1=L_w+1: one generation position supplies the first-token probability, Lw−1L_w-1 positions supply the remaining target-token probabilities, and one position supplies the closing-boundary probability. A single-token target therefore has cost multiplier 22. Under this proxy and the 1/M1/M variance scaling above, ρ^c(w) ρ_c(w) estimates the fraction of Monte Carlo compute needed to attain the same variance. For frequency ratios, this comparison does not establish equal mean squared error because finite-sample bias is not included. C.4.3 Observed-prefix potential The ZendZ^end decomposition gives the observed-prefix potential ϕ→()=f() φ( [rgb]0,0,0 x)=f( [rgb]0,0,0 x) of Section C.4.3. Its expected increment is nonzero only at prefixes that already match the substring up to its final token, so the estimator averages over that token alone. Section C.4.3 reports the median variance of the resulting estimator relative to Monte Carlo at equal sample size (M=4,000M=4,000) for the unigram targets of GPT-2, banded by their occurrences in the 4,0004,000 samples (very rare 11–99, rare 1010–9999, moderate 100100–999999, frequent ≥1,000≥ 1,000). Table 15: Substring counts, observed-prefix potential f()=#occurrences of in f( Y)=\#\occurrences of w in Y\ f()=∑tZ,tendf( Y)= _tZ^end_ w,t, Z,tend= ends at tZ^end_ w,t= 1\ w ends at t\ ϕ→⋆()=p[∑t=1||Z,tend|⪰] φ ( [rgb]0,0,0 x)= E _p\! [ _t=1^| Y|Z^end_ w,t\; |\; Y [rgb]0,0,0 x ] ϕ→()=f()=∑t≤||Z,tend φ( [rgb]0,0,0 x)=f( [rgb]0,0,0 x)= _t≤| [rgb]0,0,0 x|Z^end_ w,t Δϕ→(a∣)= ends at ||+1 in a φ(a [rgb]0,0,0 x)= 1\ w ends at | [rgb]0,0,0 x|+1 in [rgb]0,0,0 xa\ Δϕ→(eos∣)=0 φ( eos [rgb]0,0,0 x)=0 Δϕ→()=Prp[ ends at ||+1∣⪰]B φ( [rgb]0,0,0 x)=Pr_p [ w ends at | [rgb]0,0,0 x|+1 Y [rgb]0,0,0 x ] Table 16: Variance of the observed-prefix estimator relative to Monte Carlo at equal sample size (M=4,000M=4,000): median over the unigram targets of GPT-2, banded by their occurrence count in the 4,0004,000 samples (columns). model 11–99 1010–9999 100100–999999 ≥1,000≥ 1,000 GPT-2 (854854) 0.9970.997 0.9960.996 0.9930.993 0.9970.997 C.4.4 Results for expected counts Results for expected counts. Figure 6 repeats Figure 2 for the expected count θw=[Cw] _w= E [C_w] instead of the frequency; the results again show the largest gains for rare targets, with some common targets costing more than Monte Carlo at equal variance. Half the retained targets reach equal variance on at most 0.80%0.80\% of the Monte Carlo compute. Figure 6a shows the cost-adjusted variance ratio ρ^c ρ_c for the 4,1494,149 of 7,6477,647 targets whose ratio can be estimated with relative standard error at most 22 (3,7003,700 for frequencies), together with approximate per-target 95%95\% log-ratio confidence intervals. Under the per-target cost proxy, the median ratio is 0.79%0.79\% of the Monte Carlo budget, and approximately 76.8%76.8\% of the retained targets require at most 10%10\%. Replacing each ratio by the upper bound of its 95%95\% interval gives a median of 1.78%1.78\% of the Monte Carlo budget; approximately 77.2%77.2\% of the retained targets have upper bounds at most 20%20\%, and 95.4%95.4\% remain below parity with Monte Carlo. The gains are largest for rare targets. We restrict to the same 2,8312,831 targets as for frequencies, those whose count estimates have nominal RSE at most 20%20\% under a Poisson approximation, equivalently at least 2525 expected occurrences in the independent 2×1062× 10^6-sample frequency bank. The compute ratio increases strongly with the expected count per document (Spearman rs=0.87r_s=0.87), with its median rising from 2.3×10−32.3× 10^-3 in the rarest populated decade to 1.241.24 among targets with one to ten expected occurrences per document (Figure 6b). For the latter group, the median exceeds parity with Monte Carlo under the cost proxy. At the same sample size (M=4,000M=4,000), the potential estimator yields substantially more precise count estimates, with 79.7%79.7\% of targets reaching RSE ≤20%≤ 20\% and 96.5%96.5\% reaching RSE ≤50%≤ 50\%, compared with 11.4%11.4\% and 19.0%19.0\% for Monte Carlo (Figure 6c). The gap grows with token length: for targets of length 44 or more, Monte Carlo almost never reaches RSE ≤50%≤ 50\%, whereas the potential estimator does so for at least 90%90\% of targets in every length group. For Lw≥7L_w≥ 7, Monte Carlo does not observe any target even in the 4×1064× 10^6-sample reference bank, while the potential estimator gives a finite estimate for all of them. Every share agrees with the corresponding share for frequencies to within 0.140.14 percentage points. (a) ρ^c(w) ρ_c(w) with log-ratio 95% CI (b) ρ^c(w) ρ_c(w) against expected count (c) resolved by token length Figure 6: Figure 2 for expected counts. Cost-adjusted variance ratio ρ^c(w) ρ_c(w) for n-gram count estimation under the per-target cost proxy Lw+1L_w+1. (a) Empirical cumulative distribution over 4,1494,149 targets with bootstrap RSE at most 22, with approximate 95%95\% log-ratio confidence bounds. (b) Ratio versus estimated expected count per document for the 2,8312,831 targets with at least 2525 expected occurrences in the independent frequency bank. (c) Fraction of targets whose count estimate reaches each RSE threshold, grouped by token length, at the same sample size M=4,000M=4,000 for both estimators. Dashed parity lines in (a–b) mark equal variance at equal cost under the proxy; in (c), line styles distinguish RSE thresholds. C.4.5 Additional diagnostics (a) targets and token lengths per order n n N min max mean med 1 869 1 4 1.110 1 2 1584 2 5 2.128 2 3 1568 3 6 3.174 3 4 1391 4 8 4.219 4 5 1206 5 9 5.266 5 6 1029 6 10 6.320 6 all 7647 – – – – (b) share of targets at each LwL_w, per order (c) same, log scale: tails thin alike across n Figure 7: n-gram target inventory and token-length distribution. The n-gram task provides a controlled way to study targets of different lengths and, consequently, different marginalization costs, which scale with the number of GPT-2 tokens LwL_w in the target. (a) Number of targets at each order n, together with the minimum, maximum, mean, and median LwL_w of their canonical realizations. The inventory contains 7,6477,647 targets in total. (b) Distribution of token lengths within each order. Increasing n shifts the inventory toward longer targets, while GPT-2 tokenization also induces within-order variation in LwL_w. (c) The same distributions on a log scale, showing that each order additionally contains a progressively thinning tail of targets requiring more tokens, and therefore more costly marginalization. Figure 8: Cost-adjusted variance ratio against target rarity before the occurrence filter. Both panels show the frequency-estimator variance ratio, against estimated n-gram frequency (left) and expected count per document (right), before applying the occurrence filter to the independently estimated x-axis quantity. In the main figure, we retain only targets whose independent frequency estimate has RSE at most 20%20\%, corresponding to at least 2525 expected occurrences in the estimation sample, to restrict the analysis to more reliable estimates. Here, we show all 3,8353,835 targets with positive estimated frequency and a finite positive variance ratio; this remains a subset of the 7,6477,647-target inventory. Table 17: Variance reduction is strongest for rare targets, even within fixed n-gram order and token length. The 2,8312,831 targets of Figure 2b are grouped by their independently estimated frequency θ^FR θ_FR. Each entry reports the median cost-adjusted variance ratio 103ρ^c10^3 ρ_c within that frequency bin; lower values indicate larger gains over Monte Carlo. (a) The trend across all targets. (b–c) The same comparison within fixed n-gram order and canonical token length LwL_w, showing that the rarity trend persists after controlling for target length and marginalization cost. The bottom rows give the number of targets in each column and Spearman’s correlation rsr_s between θ^FR θ_FR and ρ^c ρ_c. Cells with fewer than 88 targets are omitted. (a) all targets (b) n-gram order (c) token length LwL_w frequency θ^FR θ_FR # targets 103ρ^c10^3 ρ_c 11 22 33 44 11 22 33 44 [10−8,10−7)[10^-8,10^-7) 354 1.5 2.0 1.5 1.5 1.3 – 1.6 1.5 1.3 [10−7,10−6)[10^-7,10^-6) 701 5.8 6.5 6.6 5.3 5.1 – 6.6 5.4 5.1 [10−6,10−5)[10^-6,10^-5) 725 28 35 28 25 17 42 26 25 17 [10−5,10−4)[10^-5,10^-4) 684 130 160 110 85 – 160 110 85 – [10−4,10−3)[10^-4,10^-3) 299 310 330 260 – – 330 260 – – [10−3,10−1)[10^-3,10^-1) 68 710 740 530 – – 740 530 – – targets 2,831 – 865 1,295 537 115 794 1,307 587 124 Spearman rsr_s – 0.85 0.72 0.80 0.73 0.69 0.67 0.78 0.74 0.69 The variance ratios of Figure 2 are point estimates, so we also ask how many targets are separated from parity once their own estimation error is taken into account. 18 reports this across cumulative relative-standard-error subsets. Table 18: A majority of targets in each precision subset retain a variance advantage under the cost proxy. Among targets with positive estimated frequency and a finite positive variance ratio, each row retains those with bootstrap RSE below the indicated cutoff s. We report the number of retained targets and the percentage satisfying ρ^c<1/(1+1.96s) ρ_c<1/(1+1.96s), a sufficient condition for the approximate normal interval ρ^c[1±1.96RSE] ρ_c[1± 1.96\,RSE] to lie below parity. These are conservative shares for the normal interval, rather than the log-ratio intervals in Figures 2 and 6. Expected count Frequency Max. RSE # targets 95% CI below 11 # targets 95% CI below 11 5% 84 61.9% 41 70.7% 10% 317 84.2% 240 89.6% 20% 1,048 92.2% 953 95.2% 50% 2,661 93.0% 2,583 94.7% C.5 Event Analysis of Dorman et al. (2026) C.5.1 The Automated Readability Index We extend the experiments of (Dorman et al., 2026) that estimate the Automated Readability Index. For a string Y with c()c( Y) characters, w()w( Y) words, and s()s( Y) sentences, ARI()=4.71c()w()+0.5w()s()−21.43,ARI( Y)=4.71\, c( Y)w( Y)+0.5\, w( Y)s( Y)-21.43, (31) a score calibrated so that its value roughly matches the US school grade needed to read the text. The events DjD_j are bins of this score, and the test functions are the indicators ARI()∈Dj 1\ARI( Y)∈ D_j\. Experimental setup. The experiments follow the authors’ setup. The model TinyStories-8M continues a fixed prompt for 100100 tokens, and the target is the probability that the completion’s ARIARI falls in each of 8080 bins of width 0.290.29 spanning [−8,15][-8,15]. The model is tilted by exp(−λARI) (- ) at 1919 values of λ; each tilt is sampled by transition-path sampling (200200 chains of 32,40032,400 retained states) and complemented by 200,000200,000 direct samples, and the resulting 6.68×1066.68× 10^6 states are combined with MBAR weights, which we retain unchanged. Folds are whole chains, stratified by tilting schedule. C.5.2 Potential definitions The events are the J=80J=80 readability bins DjD_j of Section C.5.1. For each bin we learn a potential ϕ→j(,)=gj((),) φ_j( [rgb]0,0,0 x; η)=g_j( ξ( [rgb]0,0,0 x); η) of the form of Eq. 13, approximating the oracle Prp[ARI()∈Dj∣⪰]Pr_p[ARI( Y)∈ D_j Y [rgb]0,0,0 x]. The prefix features () ξ( [rgb]0,0,0 x) collect the character, word, and sentence counts of [rgb]0,0,0 x together with its position in the generation window. We consider two models, a logistic model that predicts bin probabilities directly and a mixture model that forecasts the counts underlying the ARI; both are fitted by the cross-fitting scheme introduced in Section 4.3, with potential parameters fitted on the training folds using the globally estimated MBAR weights. Logistic potential. The logistic model predicts bin probabilities from the prefix features directly, ϕ→jlog(,)=p^jexp(j⊤())∑j′=1Jp^j′exp(j′⊤()),=(1,…,J), φ^log_j( [rgb]0,0,0 x; η)= p_j \! ( β_j ξ( [rgb]0,0,0 x) ) _j =1^J p_j \! ( β_j ξ( [rgb]0,0,0 x) ), η=( β_1,…, β_J), (32) where p^j p_j is the marginal probability of bin DjD_j estimated from the training data. The feature vector () ξ( [rgb]0,0,0 x) here comprises 11 features: an intercept, the normalized prefix position and its square, the inverse remaining length and its inverse square root, the character, word, and sentence counts, and the interactions between each count and the normalized position, so that the counts can play different roles at different stages of generation. We fit η by minimizing MBAR-weighted cross-entropy. Mixture-model potential. At each prefix, we linearly extrapolate the character and word counts in () ξ( [rgb]0,0,0 x) to the terminal length. We separately model the number k of additional sentence endings, since a single further sentence break can change the ARI substantially. Each value of k gives a predicted ARI μk() _k( [rgb]0,0,0 x), around which we place a Gaussian with standard deviation σk() _k( [rgb]0,0,0 x). Combining these components with their predicted probabilities πk() _k( [rgb]0,0,0 x) gives a Gaussian mixture law for the terminal ARI()ARI( Y) of the completed string, and the potential is the mass this law assigns to the bin Dj=[lj,uj)D_j=[l_j,u_j), available in closed form through the standard normal distribution function Φ : qmix(⋅∣) q_mix(\,· [rgb]0,0,0 x) =∑kπk()(μk(),σk2()), = _k _k( [rgb]0,0,0 x)\,N\! ( _k( [rgb]0,0,0 x), _k^2( [rgb]0,0,0 x) ), (33) ϕ→jmix(,) φ^mix_j( [rgb]0,0,0 x; η) =qmix(Dj∣)=∑kπk()[Φ(uj−μk()σk())−Φ(lj−μk()σk())]. =q_mix(D_j [rgb]0,0,0 x)= _k _k( [rgb]0,0,0 x) [ \! ( u_j- _k( [rgb]0,0,0 x) _k( [rgb]0,0,0 x) )- \! ( l_j- _k( [rgb]0,0,0 x) _k( [rgb]0,0,0 x) ) ]. Here η collects the extrapolation rates, mixture weights, and Gaussian widths. We fit these parameters on the combined direct and transition-path samples, giving each retained sample equal weight. To recover target-model bin probabilities, we multiply the fitted mixture density by the MBAR weight as a function of terminal ARI, integrate the resulting weighted density over each bin using fine-grid quadrature, and normalize once across the full mixture. We cap the number of additional sentence endings at 1414, a heuristic upper bound for the 100-token completions used in our experiments. C.5.3 Potential selection algorithm For each potential ϕ→ φ, we use the control-variate estimator of Eq. 19. The coefficient α may either be fixed a priori or estimated from data to minimize variance. In particular, α=1α=1 yields the potential estimator, while α=0α=0 recovers the baseline estimator regardless of the potential. Algorithm 1 gives the nested cross-fitting procedure used in Section 4.3 to select both the potential and, when learned, its coefficient using held-out folds. Algorithm 1 Potential selection by nested cross-fitting 1: Sampled strings D with stored next-token distributions; target function f. 2: P candidate potentials, each fixed or learned, with fitting procedures for learned parameters. 3: For each candidate, a fixed coefficient α (default: 11) or a coefficient-fitting procedure. 4: Outer and inner fold counts K and L. 5: Estimate θ θ of p[f(Y)]E_p[f(Y)]. 6: Split D into K equal-sized outer folds 1,…,KD_1,…,D_K. 7: for k=1,…,Kk=1,…,K do 8: ←∖kT _k ⊳ Outer training data 9: Split T into L equal-sized inner folds 1,…,LV_1,…,V_L. 10: for j=1,…,Pj=1,…,P do 11: for ℓ=1,…,L =1,…,L do 12: ←∖ℓA _ 13: Fit candidate j’s potential on A, if needed, to obtain ϕφ. 14: Set α to its fixed value, or fit it on A. 15: Evaluate θ^(yi,ϕ,α) θ(y_i;φ,α) for every yi∈ℓy_i _ . 16: vjℓ←v_j ← sample variance of these estimates. 17: sj←1L∑ℓ=1Lvjℓs_j← 1L _ =1^Lv_j ⊳ Validation score 18: j⋆←argminj=1,…,Psjj ← _j=1,…,Ps_j 19: Refit the selected potential on T, if needed, to obtain ϕk _k. 20: Set αk _k to the selected candidate’s fixed value, or refit it on T. 21: θ^k←1|k|∑yi∈kθ^(yi,ϕk,αk) θ_k← 1|D_k| _y_i _k θ(y_i; _k, _k) ⊳ Held-out estimate 22: return θ^=1K∑k=1Kθ^k θ= 1K _k=1^K θ_k K=5K=5 outer folds1,…,B y_1,…, y_B1122334455Selectargminj≤pVar^CV(θ^(ϕ→j)) _j≤ p Var_ 0.7$ $( θ( φ_j))11223344Nested CVRefitj⋆j 0.6$ $1122334411ϕ→j⋆(,^1) φ_j 0.6$ $( [rgb]0,0,0 x; η_ 0.8$ 1$)22ϕ→j⋆(,^2) φ_j 0.6$ $( [rgb]0,0,0 x; η_ 0.8$ 2$)33ϕ→j⋆(,^3) φ_j 0.6$ $( [rgb]0,0,0 x; η_ 0.8$ 3$)44ϕ→j⋆(,^4) φ_j 0.6$ $( [rgb]0,0,0 x; η_ 0.8$ 4$)55ϕ→j⋆(,^5) φ_j 0.6$ $( [rgb]0,0,0 x; η_ 0.8$ 5$) Figure 9: Nested cross-fitting, our in-sample selection among p candidate potentials (Algorithm 1). Each of K=5K=5 outer folds is held out in turn (red); on the remaining folds the candidates are compared by inner cross-validation, the winner j⋆j is refitted on all remaining data and applied to the held-out fold. C.5.4 Compute costs and additional results Compute. All potential fitting, scoring, and selection run on CPUs, with no GPU computation beyond sampling. The logistic potential has 880880 coefficients (1111 features across 8080 bins), while the mixture model stores empirical rates, mixture weights, and residual widths. We use five outer folds, with the four remaining folds serving as the inner folds. Each inner fit uses three folds: exchanging the outer held-out fold and the inner validation fold leaves the training set unchanged. Thus, the 5×4=205× 4=20 inner evaluations require only (52)=10 52=10 distinct fits, which we reuse across outer splits. Together with 55 outer refits and 11 full-data fit, this gives 1616 parameter sets per model. Fitting these takes 1414 minutes for logistic regression and 44 minutes for the mixture model on a single 3232-core node. Scoring both potentials using the stored next-token distributions requires ten parallel jobs, each using 3232 cores for 5050–9595 minutes. Additional results. Figure 9 sketches the selection procedure and Figure 10 shows the variance of the selected estimator over the bins on which MBAR’s own variance is well measured, Figure 11 the same view for each single candidate, Figure 12 extends the comparison to all 6464 populated bins with the accompanying diagnostics, and 19 lists the chain-level variance, sample-size efficiency, and MBAR-equivalent sample size of every estimator. Figure 10: Variance ratio of the selected potential on the event-analysis experiment: chain-level variance of the cross-fitted estimator relative to MBAR, per ARI bin, over the 3939 bins 2020–5858: the 3838 bins on which MBAR’s own variance is itself measured to within 20%20\%, plus bin 5757 (at 25%25\%), kept so that the range is contiguous; shading gives the bin probability, the dashed line the median (0.880.88) and the band its 95%95\% interval. All 6464 populated bins are in Figure 12. (a) terminal potential (b) terminal potential with α (c) logistic model (d) mixture model Figure 11: The same view as Figure 10 for each single candidate, each refitted without the held-out fold and applied to every bin: chain-level variance relative to MBAR per ARI bin over the same 3939 bins (3838 by the 20%20\% rule plus bin 5757), with the bin probability as shading and the candidate’s own median over these bins with its 95%95\% interval. The terminal potential and its calibrated version remove 8%8\%, the logistic model 11%11\%, and the mixture model 13%13\%, against 12%12\% for the selection, which picks the mixture model in four folds and the logistic model in one. Figure 12: Variance reduction across ARI bins. (a) variance of the cross-fitted estimator relative to MBAR for all 6464 populated bins. Filled markers denote the 3838 bins for which the MBAR variance is estimated with RSE below 20%20\%; over these bins, the median variance ratio is 0.880.88 with a 95%95\% bootstrap interval of [0.84,0.90][0.84,0.90], and over all 6464 populated bins it is 0.880.88 with interval [0.86,0.91][0.86,0.91]. Ratios in the remaining tail bins are much noisier because the MBAR variance itself is poorly estimated. Whiskers show ±1.96± 1.96 bootstrap standard deviations from 1,0001,000 paired chain-bootstrap replicates. (b) RSE of the MBAR variance and mean estimates across bins, showing that uncertainty in the variance estimate grows sharply in the tails. (c) Variance ratios against the RSE of the MBAR variance estimate. Among the 3838 reliable bins, the ratio shows no systematic dependence on baseline uncertainty (Spearman ρ=0.05ρ=0.05, p=0.78p=0.78), indicating that the observed reduction is not driven by noise in the MBAR variance estimate. Estimator Var. ratio vs. MBAR [95% CI] Sample-size eff. [95% CI] MBAR-equiv. samples (×106× 10^6) [95% CI] Nested selection 0.88 [0.86, 0.91] 1.14 [1.10, 1.17] 7.58 [7.34, 7.80] Terminal 0.92 [0.90, 0.94] 1.09 [1.06, 1.11] 7.29 [7.09, 7.40] Terminal, α 0.92 [0.90, 0.95] 1.09 [1.06, 1.11] 7.27 [7.05, 7.40] Logistic 0.90 [0.86, 0.91] 1.12 [1.09, 1.16] 7.46 [7.31, 7.75] Mixture model 0.87 [0.84, 0.90] 1.15 [1.11, 1.19] 7.69 [7.43, 7.92] Table 19: Event analysis: chain-level variance of each estimator relative to MBAR, median over all 6464 populated ARI bins, with 95%95\% intervals from 1,0001,000 paired chain-bootstrap replicas (fitted potentials, selections, and coefficients held fixed). Sample-size efficiency: median over bins of Var(MBAR)/Var(estimator)Var(MBAR)/Var(estimator); MBAR-equivalent samples: that factor times the actual sample size N=6.68×106N=6.68× 10^6 states.