stat.MLOct 6, 2026

Steering Diffusion Models to Rare Events with Sequential Monte Carlo

Authors: Aavash Subedi, Tim Reichelt, Christopher Williams, Philip Stier, Yee Whye Teh, Saifuddin Syed

Organizations: Department of Statistics University of Oxford · Department of Physics University of Oxford · Department of Statistics University of British Columbia

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Diffusion models are increasingly used as surrogates for expensive simulators in weather prediction, molecular dynamics, and materials design. In these models, computing the probability p0[E]p_0[E] of an event EE is difficult, especially when the event of interest is rare. A stable estimate using Monte Carlo becomes computationally intractable, requiring a growing sample size ∝ ⁣1/p0[E]\propto\!1/p_0[E] to compensate for an increasing rarity. In this paper, we present Diffusion Importance Sampling of Rare Events or DireSMC, a sequential Monte Carlo scheme that guides a population of weighted samples towards the rare event, giving access not only to samples but also to a calibrated estimate of its probability. We set up our guidance using an analytical relaxation of the event set, allowing the method to easily extend to a wide range of user-defined rare events. We validate our method on a toy problem with analytical solutions and on a score-based climate emulator, where we obtain accurate rare-event probabilities on a range of rarities from 10−310^{-3} to 10−510^{-5}, achieving net speed-ups of 9×9\times to 1413×1413\times over Monte Carlo.

Figures & tables

Figure 1: We steer particles ( ∙\bullet ) from prefp_{\mathrm{ref}} towards π0\pi_{0} , which concentrates mass in the rare-event region of interest, using a guidance drift {\color[rgb]{1,0,0}\nabla_{x}r_{k}} . Particle size reflects the importance weights accumulated along the trajectory, which are used to resample ( ×\times ) and form a rare event estimate, p0[E]p_{0}[E] .
Figure 2: Left: CV plotted against NN for a p0[x1>10]p_{0}[x_{1}>10] problem at a fixed KK . Rejuvenation tracks the lowest CV across all NN . Middle: we compare our estimator against Manshausen et al. (2026) across rarities at a compute-normalised NN . We consistently track a lower CV. Right: measured bias scaling with rarity using exact Doob’s steering and analytical scores, highlighting the bias of Manshausen et al. (2026) , DireSMC tracks the floor, see Tables 23 and App. J.1.2 .
Figure 3: Surface-temperature anomaly (K) for July at ΔT=+1.5\Delta T=+1.5 K over the Pacific Northwest (PNW) target region. Left: average tas field over 10710^{7} samples. Remaining panels: the highest-weight crossing particle produced by the guided sampler at different exceedance thresholds τ\tau and regional average temperature ϕ(⋅)\phi(\cdot) of the sample. p⋆p^{\star} are MC estimates; p^\hat{p} is obtained from the guided estimator; full results in Table 15 . The variance ratio (Var. Ratio) at matched compute is defined as Var(p^unguided)/Var(p^guided)\mathrm{Var}(\hat{p}_{\text{unguided}})/\mathrm{Var}(\hat{p}_{\text{guided}}) .
p⋆p^{\star}p^/p⋆\hat{p}/p^{\star}CVd⋆d^{\star}d^/d⋆\hat{d}/d^{\star}Speed Up ↑\uparrow
3.14×10−33.14\times 10^{-3}0.970.97 ±\pm0.020.020.070.0731311.041.04 ±\pm0.020.029×9\times
2.12×10−42.12\times 10^{-4}0.990.99 ±\pm0.030.030.110.112122121.021.02 ±\pm0.030.0347×47\times
1.39×10−51.39\times 10^{-5}1.041.04 ±\pm0.020.020.080.08138413841.061.06 ±\pm0.060.061413×1413\times
Table 1: Probability estimates of three compound hot and dry extreme events over the PNW region. The dependence ratio dd is defined in the main text, d^\hat{d} denotes that estimated by our method. Net speed-up is the crude-MC cost of matching the CV of the guided run divided by the measured per-sample cost of guidance (Table 21 ); thresholds τ\tau are given in App. H.1 . Mean ±\pm standard error over 16 seeds.
Figure 4: Estimated p^0[E]\hat{p}_{0}[E] versus global temperature anomaly ΔT\Delta T for our method versus unguided MC at equal compute. Points are the mean over 16 runs, bars are ±1\pm 1 std. Bars reaching the axis floor indicate levels where some MC runs record no event; for ΔT≤1.0\Delta T\leq 1.0 each MC run records at most one event. See Tables 19 and 20 .
Appendix figures & tables41 assets

Supplementary material from the paper’s appendix.

Appendix

NNp^SNIS\hat{p}_{\text{SNIS}}Relative BiasZ^0,SNIS\hat{Z}_{0,\text{SNIS}}
1004.6035×10−14.6035\times 10^{-1}+6,593%+6,593\%4.873×10−14.873\times 10^{-1}
1,0003.5854×10−13.5854\times 10^{-1}+5,113%+5,113\%3.770×10−13.770\times 10^{-1}
10,0003.1680×10−13.1680\times 10^{-1}+4,506%+4,506\%3.367×10−13.367\times 10^{-1}
25,0001.9804×10−11.9804\times 10^{-1}+2,779%+2,779\%2.102×10−12.102\times 10^{-1}
1006.9250×10−36.9250\times 10^{-3}+0.68%+0.68\%(analytical Z0,δZ_{0,\delta} )
Appendix
Table 2: Self-normalised importance sampling (SNIS) estimates vs. sample size ( NN ).
 ckc_{k}vkv_{k}ζk2=ck2/vk\zeta^{2}_{k}=c^{2}_{k}/v_{k}
Euler–Maruyama (VP/VE)σtk2Δtk\sigma^{2}_{t_{k}}\Delta t_{k}σtk2Δtk\sigma^{2}_{t_{k}}\Delta t_{k}σtk2Δtk\sigma^{2}_{t_{k}}\Delta t_{k}
exponential, VP ( ϵ^\hat{\epsilon} frozen)2σˉtk−1σˉtk(eΔℓk−1)2\bar{\sigma}_{t_{k-1}}\bar{\sigma}_{t_{k}}(e^{\Delta\ell_{k}}-1)σˉtk−12(e2Δℓk−1)\bar{\sigma}^{2}_{t_{k-1}}(e^{2\Delta\ell_{k}}-1)4 σˉtk2tanh⁡(Δℓk/2)4\,\bar{\sigma}^{2}_{t_{k}}\tanh(\Delta\ell_{k}/2)
exponential, VE ( x^\hat{x} frozen)(1−κk) σˉtk2(1-\kappa_{k})\,\bar{\sigma}^{2}_{t_{k}}σˉtk−12(1−κk)\bar{\sigma}^{2}_{t_{k-1}}(1-\kappa_{k})σˉtk2(e2Δℓk−1)\bar{\sigma}^{2}_{t_{k}}(e^{2\Delta\ell_{k}}-1)
exponential, VE ( ϵ^\hat{\epsilon} frozen)2 σˉtk−1σˉtk(eΔℓk−1)2\,\bar{\sigma}_{t_{k-1}}\bar{\sigma}_{t_{k}}\bigl(e^{\Delta\ell_{k}}-1\bigr)σˉtk−12(e2Δℓk−1)\bar{\sigma}^{2}_{t_{k-1}}\bigl(e^{2\Delta\ell_{k}}-1\bigr)4 σˉtk2tanh⁡(Δℓk/2)4\,\bar{\sigma}^{2}_{t_{k}}\tanh(\Delta\ell_{k}/2)
Appendix
Table 3: ζk2=ck2/vk\zeta_{k}^{2}=c_{k}^{2}/v_{k} for different integrators and noise schedules.
Figure 5: Soft reward and its guidance for a threshold event E7={x:ϕ(x)≥7}E_{7}=\{x:\phi(x)\geq 7\} shown for three smoothing levels δ\delta . Left. The soft indicator relaxes the hard indicator. Right. The guidance ∇xr0,δ(x)\nabla_{x}r_{0,\delta}(x) points toward the event set (shaded) and decays inside it, a smaller δ\delta gives a stronger guidance.
Figure 6: Base distribution p0(x)p_{0}(x) of the correlated 2D GMM (log-density, so all modes are visible) with the analytic x1x_{1} and x2x_{2} marginals. The two right-hand modes ( i=5,6i=5,6 ) form the tilted “X” of opposite-sign correlation that gives the rare set a bi-modal conditional.
iiwiw_{i}μi\mu_{i}Σi\Sigma_{i}
10.100.10(−6,−2)(-6,-2)[0.360.300.300.81]\big[\begin{smallmatrix}0.36&0.30\\ 0.30&0.81\end{smallmatrix}\big]
20.180.18(−4, 4)(-4,\ 4)[0.81−0.25−0.250.49]\big[\begin{smallmatrix}0.81&-0.25\\ -0.25&0.49\end{smallmatrix}\big]
30.140.14(−1,−5)(-1,-5)[0.490.200.200.81]\big[\begin{smallmatrix}0.49&0.20\\ 0.20&0.81\end{smallmatrix}\big]
40.220.22(−10, 3)(-10,\ 3)[0.040.040.040.09]\big[\begin{smallmatrix}0.04&0.04\\ 0.04&0.09\end{smallmatrix}\big]
50.100.10( 4, 3)(\ 4,\ 3)[2.401.851.852.00]\big[\begin{smallmatrix}2.40&1.85\\ 1.85&2.00\end{smallmatrix}\big]
60.160.16( 4,−3)(\ 4,-3)[2.40−1.85−1.852.00]\big[\begin{smallmatrix}2.40&-1.85\\ -1.85&2.00\end{smallmatrix}\big]
Appendix
Table 4: Parameters of the correlated 2D GMM. The right-hand modes ( i=5,6,7i=5,6,7 ) carry the rare-event tail; modes 55 and 66 have correlation ρ≈+0.85\rho\approx+0.85 and −0.85-0.85 respectively.
τ\taupdata[E]p_{\mathrm{data}}[E]pmodelp_{\rm model} (MC)pmodel/pdata[E]p_{\rm model}/p_{\mathrm{data}}[E]rel. s.e.path integralratio
76.88 ⁣× ⁣10−36.88\!\times\!10^{-3}7.13 ⁣× ⁣10−37.13\!\times\!10^{-3}1.0370.08%1.0751.037
81.28 ⁣× ⁣10−31.28\!\times\!10^{-3}1.38 ⁣× ⁣10−31.38\!\times\!10^{-3}1.0810.18%1.1561.069
91.62 ⁣× ⁣10−41.62\!\times\!10^{-4}1.87 ⁣× ⁣10−41.87\!\times\!10^{-4}1.1530.50%1.2361.072
101.40 ⁣× ⁣10−51.40\!\times\!10^{-5}1.72 ⁣× ⁣10−51.72\!\times\!10^{-5}1.2291.64%1.2671.031
Appendix
Table 5: Tail bias of the learned 2D-GMM score. pmodelp_{\rm model} (MC) is obtained from a large MC reference pool using unguided draws. Path integral measures the ratio obtained by integrating the pointwise score error along the probability-flow ODE (Eq. ( 245 )) and reweighing the exact conditional. The two induce different laws but agree to within 7%7\% across three decades of rarity. pdata[E]p_{\mathrm{data}}[E] is the analytical probability.
KKEM biasSEEDS-1 bias∣EM∣/∣SEEDS∣|\mathrm{EM}|/|\mathrm{SEEDS}|CV ↓\downarrowk^↓\hat{k}\downarrowESSE↑\mathrm{ESS}_{E}\uparrow
Uniform in tt
10+64.02% ±\pm 11.31%-15.98% ±\pm 1.88%4.010.338/0.1090.85/0.460.032/0.102
15+48.61% ±\pm 8.39%-16.23% ±\pm 1.28%3.000.277/0.0750.70/0.450.063/0.189
20+54.90% ±\pm 16.50%-15.13% ±\pm 0.97%3.630.522/0.0560.68/0.430.099/0.248
30+23.41% ±\pm 1.65%-13.62% ±\pm 0.75%1.720.065/0.0430.56/0.400.176/0.375
50+16.75% ±\pm 1.38%-8.36% ±\pm 0.63%2.000.058/0.0340.55/0.450.315/0.465
Appendix
Table 6: Signed discretisation bias versus step count KK for the two integrators, under exact Doob guidance, at different grids, uniform in tt (top block) and uniform in the half-log SNR ℓt=log⁡(αt/σˉt)\ell_{t}=\log(\alpha_{t}/\bar{\sigma}_{t}) (bottom block). Positive values are overestimates. Paired columns are EM/SEEDS-1 throughout. k^\hat{k} is the Pareto Smoothed Importance Sampling shape estimate ( Vehtari et al., 2024 ) on the weights, and ESSE\mathrm{ESS}_{E} is the ESS of the final particles, measured after thresholding to τ≥10\tau\geq 10 .
Figure 7: Discretisation bias for the three integrators, Euler-Maruyama, SEEDS-1 and DDIM. We unconditionally draw 1e8 samples from the base model at different discretisation KK and evaluate the bias at different rarities. The figure on the right shows the number of discretisation required to get a <10%<10\% bias for different rarities, showing that the SDE requires less steps than the ODE integrator. Top row uses the uniform log-SNR Δℓ\Delta\ell schedule and the bottom is uniform in time Δt\Delta t .
Figure 8: Measured bias versus KK with 1000 particles, analytical scores and Doob’s transition kernels. See App G.2.2 . We show the bias of our estimate falls with a O(1/K)O(1/K) leading order term. We compare both the EM integrator and the SEEDs-1 integrator in the noise prediction form. We compare two different noise schedules: uniform in time and uniform in log-SNR ℓ\ell . Error-bars represent standard error of the mean. Experiment is conducted the τ=10\tau=10 threshold exceed problem. See Table 6 .
Figure 9: Measured bias versus KK with 10,000 particles and the deployed guidance and no rejuvenation for the learned score. We compare the discretisation bias at different rarities, we compare against a 10810^{8} unconditional MC draw of the exponential integrator at KK as our reference. Hence, this measures how quickly the estimates converge to the reference with KK .
τ\tauscorep⋆p^{\star}p^/p⋆\hat{p}/p^{\star}CV ↓\downarrowKS x1{}_{x_{1}}KS x2{}_{x_{2}}EMD
7analytical6.88 ⁣× ⁣10−36.88\!\times\!10^{-3}0.980±0.0120.980\pm 0.0120.0600.881.331.52
8analytical1.28 ⁣× ⁣10−31.28\!\times\!10^{-3}0.966±0.0120.966\pm 0.0120.0611.131.101.06
9analytical1.62 ⁣× ⁣10−41.62\!\times\!10^{-4}0.974±0.0160.974\pm 0.0160.0831.121.000.97
10analytical1.40 ⁣× ⁣10−51.40\!\times\!10^{-5}0.954±0.0130.954\pm 0.0130.0690.941.051.46
7learned7.13 ⁣× ⁣10−37.13\!\times\!10^{-3}0.978±0.0120.978\pm 0.0120.0601.111.201.14
8learned1.38 ⁣× ⁣10−31.38\!\times\!10^{-3}0.985±0.0130.985\pm 0.0130.0640.971.000.73
Appendix
Table 7: Guided generation on the GMM with learned and analytical scores. Each arm is scored against the law it samples from, p⋆p^{\star} ; either pdata[E]p_{\mathrm{data}}[E] (analytical) or an MC estimate of p0[E]p_{0}[E] (learned). Sample quality metrics are evaluated against a reference from the same law, presented as a ratio against a noise-floor. Estimates p^/p⋆\hat{p}/p^{\star} are biased to the MADM kernel’s floor (Fig. 13 ). We evaluate the two-sample Kolmogorov-Smirnov (KS) test statistic, the Earth Mover’s Distance (EMD) and the coefficient of variation CV2=Var(pE^)/pE2\text{CV}^{2}=\mathrm{Var}(\hat{p_{E}})/p_{E}^{2} .
Figure 10: Conditional-sample coverage on the 2D GMM for the analytical (a) and learned (b) scores. Blue is the path-weighted guided cloud ( 20482048 samples; pooled over 88 seeds for the learned score); red is a MC draw from the conditional p(x∣x1≥τ)p(x\mid x_{1}\geq\tau) . Left: the joint cloud, with the threshold dashed. Centre, right: the x1x_{1} and x2x_{2} marginals on a log density scale. Corresponding quantitative KS/EMD metrics in Table 7 .
Figure 11: Effective sample size of the path weights along diffusion time, without resampling (analytic score, left; learned score, right), for τ∈{6,7,8,9,10}\tau\in\{6,7,8,9,10\} . Settings as Table 7 . Time runs right to left, from noise ( t=1t=1 ) to data ( t=0t=0 ).
τ\tauarmCV ↓\downarrowKS x1{}_{x_{1}}KS x2{}_{x_{2}}EMD
7no resampling0.0670.881.521.50
resampling0.0621.050.990.83
resampling ++ rejuv. ( M=3M=3 )0.0600.881.331.52
8no resampling0.0831.141.351.17
resampling0.0981.120.950.77
resampling ++ rejuv. ( M=3M=3 )0.0611.131.101.06
Appendix
Table 8: Selection and rejuvenation on the 2D GMM under the analytical score, 2424 seeds, all other settings as in Table 7 . Resampling and rejuvenation are ablated at different threshold exceedance rarities. Weighted shape-metrics are computed against a conditional reference. Rejuvenation helps reduce the CV across all rarities; without it resampling alone can adversely prune lineages.
eventarmp^/p⋆\hat{p}/p^{\star}CV ↓\downarrowKS 0KS 1EMD
x0>7x_{0}>7no resampling0.9520.0700.741.001.09
resampling0.9880.0930.941.251.47
resampling ++ rejuv. ( M=3M=3 )0.9850.0611.011.171.28
gated ℓ≲0.05\ell\lesssim 0.050.9830.0600.891.151.11
gated ℓ≲−0.21\ell\lesssim-0.210.9780.0601.111.201.14
gated ℓ≲−0.46\ell\lesssim-0.460.9920.0820.901.331.37
Appendix
Table 9: Resampling and rejuvenation ablation for the learned score network. Distributional metrics are presented as ratios to the reference pool’s split-half noise floor. Gate levels are quoted as the half-log SNR ℓ\ell below which the rejuvenation move triggers; the adopted gate ℓ≲−0.21\ell\lesssim-0.21 is t≥0.30t\geq 0.30 . The no-resampling and resampling-only arms run at λ≡1\lambda\equiv 1 . Rejuvenation tracks the lowest CV.
Figure 12: (Left) Magnitude of the antisymmetric part of the Jacobian of the score network, measured as a Frobenius norm for 1024 particles for the different threshold exceedance problems of Table 9 . (Right) Same value as the left multiplied by ρ\rho for τ>10\tau>10 threshold exceedance problem, for different calibrated step sizes ρ\rho , denoted by their average acceptance rate. The shaded band marks the range of gates tested and the dashed line the adopted gate.
Figure 13: Bias introduced by MADM for the two integrators, using analytical scores and Doob’s guidance with M=3M=3 for the first two panels. Left & Middle Injected bias affects rarer samples more, with all biases decreasing as K→∞K\rightarrow\infty . The two integrators have a different error profile. Right Bias plotted for the τ=10\tau=10 threshold exceedance problem. We find that the error grows with MM . Error bars are due to seed variability.
Figure 14: Illustrative example of box boundaries. Plotted against a backdrop of the unconditional density of the 2D-GMM. Quantitative results of rare-event sampling are given in Table 10 below.
eventvariantp⋆p^{\star}p^/p⋆\hat{p}/p^{\star}CV ↓\downarrowKS 0KS 1EMD
[6,∞)×[1.5,∞)[6,\infty)\times[1.5,\infty)no resampling9.97 ⁣× ⁣10−39.97\!\times\!10^{-3}0.9550.0660.950.981.01
resampling1.0040.1071.370.870.98
resampling ++ rejuv. ( M=3M=3 )1.0000.0700.890.910.96
gated ℓ≲−0.21\ell\lesssim-0.211.0000.0700.940.891.00
[6,∞)×[1.5,∞)[6,\infty)\times[1.5,\infty) ∪ [6,∞)×(−∞,−1.5]\cup\,[6,\infty)\times(-\infty,-1.5]no resampling2.60 ⁣× ⁣10−22.60\!\times\!10^{-2}0.9660.0571.151.151.09
resampling0.9790.0550.991.080.95
Appendix
Table 10: Box events on the learned score. Every region is drawn on the model density in Fig. 14 . p⋆p^{\star} is the K=1000K=1000 chain’s own tail (Table 5 pool). We compare the estimate p^\hat{p} and the shape metrics of the final particles, presented as a ratio against a noise floor.
schedulebase ↓\downarrow+ resampling ↓\downarrow+ MADM ↓\downarrow( Manshausen et al., 2026 ) ↓\downarrow
analytic score
λ≡1\lambda\equiv 10.150.230.090.29
ton=1.0, p=2t_{\rm on}=1.0,\ p=20.280.320.080.68
ton=0.9, p=2t_{\rm on}=0.9,\ p=2 (deployed)0.350.330.070.83
ton=0.7, p=2t_{\rm on}=0.7,\ p=20.620.400.101.26
ton=0.5, p=2t_{\rm on}=0.5,\ p=21.730.910.131.94
Appendix
Table 11: Schedule ablation at τ=10\tau=10 , K=1000K=1000 , N=1000N=1000 : CV of p^\hat{p} for the tempering schedule λ(t)=1−(t/ton)p\lambda(t)=1-(t/t_{\rm on})^{p} . Bold marks the minimum of each column. We find that the no-rejuvenation variants converged on the same tempering schedule λt≡1\lambda_{t}\equiv 1 compared to the MADM variant. All schedules with rejuvenation yield a lower CV compared to without.
Short nameDescriptionUnits
tas22 m air temperatureK
hurs22 m relative humidity%
prPrecipitation ratemm day -1
sfcWind1010 m wind speedm s -1
Appendix
Table 12: Emulated variables with their corresponding CMIP6 short names, descriptions and units.
  MarginalUpper tailLower tail
VariableUnitMeanSD10−310^{-3}10−410^{-4}10−510^{-5}10−310^{-3}10−410^{-4}10−510^{-5}
tasK+2.577+2.5771.0231.0235.8135.8136.4986.4987.0967.096−0.539-0.539−1.174-1.174−1.727-1.727
prmm/day−0.250-0.2500.4500.4501.6981.6982.3112.3112.8422.842−1.334-1.334−1.517-1.517−1.650-1.650
hurs%−3.275-3.2754.2204.22010.51910.51913.48513.48516.05416.054−15.838-15.838−18.387-18.387−20.554-20.554
sfcWindm/s+0.193+0.1930.1860.1860.7730.7730.8950.8951.0041.004−0.387-0.387−0.508-0.508−0.613-0.613
Appendix
Table 13: Tail thresholds for the four emulator variables, from N=9,011,200N=9{,}011{,}200 unguided Monte Carlo draws over RR at ΔT=+1.5\Delta T=+1.5 K in July. Upper-tail entries satisfy P(ϕ≥τ)=pP(\phi\geq\tau)=p , lower-tail entries P(ϕ≤τ)=pP(\phi\leq\tau)=p . Binomial relative standard errors are 1.1%1.1\% , 3.3%3.3\% and 10.5%10.5\% at p=10−3,10−4,10−5p=10^{-3},10^{-4},10^{-5} respectively. The guided runs use the same thresholds.
Figure 15: Discretisation bias for the integrator for the τ1\tau_{1} PNW heatwave event using the SEEDS-1 integrator, with 1024 samples at each point and without any rejuvenation. Dotted line represents the reference pool. We find that for SEEDS-1 and Euler-Maruyama integrators using K=200K=200 saturates the discretisation bias. We could not compare against a large KK unconditional reference due to prohibitive costs. Bars represent ±1\pm 1 standard error of the mean over 8 shared seeds across both integrators.
 SymbolValue
Schedules
Tempering onsetσˉon\bar{\sigma}_{\mathrm{on}}4040
Tempering peakλmax⁡\lambda_{\max}1.251.25 at σˉ=2\bar{\sigma}=2
Tempering relaxationσˉend\bar{\sigma}_{\mathrm{end}}0.30.3
Smoothing width, transportδhi\delta_{\mathrm{hi}}0.5 SD0.5\,\mathrm{SD}
Smoothing width, floorδlo\delta_{\mathrm{lo}}0.01 SD0.01\,\mathrm{SD}
Appendix
Table 14: Deployed schedule and sampler configuration for the climate emulator experiments. σˉ\bar{\sigma} is the marginal noise standard deviation. Selected once on the marginal 10−310^{-3} PNW hot and dry events and applied unchanged to other experiments. ‡ δ\delta is quoted in units of the observable’s unconditional standard deviation.
 Hot: ϕtas(x)≥τ\phi_{\texttt{tas}}(x)\geq\tauDry: ϕhurs(x)≤τ\phi_{\texttt{hurs}}(x)\leq\tau
  p^/p⋆\hat{p}/p^{\star}  p^/p⋆\hat{p}/p^{\star} 
 p⋆p^{\star}UnguidedGuidedVRp⋆p^{\star}UnguidedGuidedVR
τ1\tau_{1}1.00e-31.00\text{e-}30.90±0.450.90\pm 0.450.98±0.08\mathbf{0.98\pm 0.08}34×34\times1.00e-31.00\text{e-}31.03±0.441.03\pm 0.441.00±0.07\mathbf{1.00\pm 0.07}37×37\times
τ2\tau_{2}1.00e-41.00\text{e-}41.30±1.671.30\pm 1.671.01±0.08\mathbf{1.01\pm 0.08}411×411\times1.00e-41.00\text{e-}40.52±0.930.52\pm 0.931.02±0.22†\mathbf{1.02\pm 0.22}^{\dagger}18×18\times
τ3\tau_{3}1.01e-51.01\text{e-}52.57±7.032.57\pm 7.031.05±0.13\mathbf{1.05\pm 0.13}2722×2722\times1.01e-51.01\text{e-}52.57±7.032.57\pm 7.031.06±0.12\mathbf{1.06\pm 0.12}3567×3567\times
Appendix
Table 15: Marginal extreme-event probabilities over the PNW target region in July at ΔT=+1.5\Delta T=+1.5 K, for hot ( ϕtas≥τ\phi_{\texttt{tas}}\geq\tau ) and dry ( ϕhurs≤τ\phi_{\texttt{hurs}}\leq\tau ) events. For each variable, τ1,τ2,τ3\tau_{1},\tau_{2},\tau_{3} are chosen so that the event is successively rarer (exact values in App. H ). Ground-truth p⋆p^{\star} are computed with ∼107\sim 10^{7} MC samples. Unguided Monte Carlo and our guided method use the same computational budget: each unguided estimate uses 48134813 draws ( 4.7×4.7\times the 10241024 guided particles, Table 21 ), over 16 disjoint replicates. Estimated probabilities, p^\hat{p} are reported as a ratio p^/p∗\hat{p}/p^{*} with mean ±1\pm 1 standard deviation over 16 seeds. The variance ratio (VR) is Var(p^unguided)/Var(p^guided)\mathrm{Var}(\hat{p}_{\text{unguided}})/\mathrm{Var}(\hat{p}_{\text{guided}}) . The number of unguided replicates with no event is 1 (hot) and 0 (dry) of 16 at τ1\tau_{1} , 9 and 12 at τ2\tau_{2} , and 14 for both at τ3\tau_{3} , where the other two replicates record a single event each, so the VR at τ3\tau_{3} is a rough estimate. † Includes one collapsed run (App. H.2.1 ). The hot results are shown in Fig. 3 .
 p⋆p^{\star}armp^/p⋆\hat{p}/p^{\star}CV ↓\downarrowESS/N\mathrm{ESS}/N k^\hat{k}W1W_{1}Speed-up ↑\uparrow
      tashurs(95% CI)
Hot10−310^{-3}no resampling1.024±0.0511.024\pm 0.0510.1990.1990.060.06 0.490.491.031.03–5.2×5.2\times ( 2.82.8 – 12×12\times )
 adaptive resampling0.982±0.0190.982\pm 0.0190.0780.078–1.541.54–34×34\times ( 2020 – 60×60\times )
 + MADM moves0.954±0.0200.954\pm 0.0200.0820.082–1.401.40–27×27\times ( 1717 – 44×44\times )
10−410^{-4}no resampling1.020±0.0491.020\pm 0.0490.1910.1910.030.03 0.570.571.071.07–56×56\times ( 3232 – 94×94\times )
 adaptive resampling1.010±0.0211.010\pm 0.0210.0820.082–0.950.95–311×311\times ( 129129 – 616×616\times )
Appendix
Table 16: Ablation of resampling and rejuvenation on the climate emulator for the hot, dry, and joint hot-and-dry events. Same parameters as in the main text; details in App. H.2.1 . Ratios are mean ±\pm one standard error over 16 seeds. W1W_{1} is a weighted distance between the ϕ(⋅)\phi(\cdot) of generated samples and draws from a reference, scaled against a null floor. The relative error of p⋆p^{\star} is 1.1/3.3/10.5%1.1/3.3/10.5\% at the three marginal rarities and 0.6/2.3/8.9%0.6/2.3/8.9\% at the three joint rarities. CV is estimated across p^\hat{p} estimates; k^\hat{k} is estimated per seed and we report the median over the 16 seeds. † Includes one collapsed run; without it p^/p⋆=0.968\hat{p}/p^{\star}=0.968 and CV =0.124=0.124 . k^>0.7\hat{k}>0.7 represents weights that have no finite second moment ( Vehtari et al., 2024 ) , so the CV estimates, and the speed-ups of those rows, are optimistic. Speed-up is over crude MC at matched CV, including each arm’s measured wall-clock (App. H.4 ); in brackets, its 95% bootstrap CI (BCa, 20,000 resamples of the 16 seeds).
Figure 16: Top-weight crossing fields over the Pacific Northwest (PNW) target region produced by the guided sampler at different exceedance threshold τ\tau . Guided samples based on the tas variable alone. Same setup as reference Table 15 . pp denotes p0[E]p_{0}[E] and ϕ(⋅)\phi(\cdot) is the regional average temperature.
Figure 17: The same fields as Fig. 16 on the global domain.
marginal ppτtas\tau_{\mathrm{tas}} (K)τhurs\tau_{\mathrm{hurs}} (%)pjointp_{\mathrm{joint}}pindepp_{\mathrm{indep}}d⋆d^{\star}
10−110^{-1}3.8908−8.6085-8.6085(5.028±0.007)×10−2(5.028\pm 0.007)\times 10^{-2}1.00×10−21.00\times 10^{-2}5.03±0.015.03\pm 0.01
10−210^{-2}4.9907−12.7822-12.7822(3.138±0.019)×10−3(3.138\pm 0.019)\times 10^{-3}1.00×10−41.00\times 10^{-4}31.4±0.231.4\pm 0.2
10−310^{-3}5.8128−15.8382-15.8382(2.121±0.049)×10−4(2.121\pm 0.049)\times 10^{-4}1.00×10−61.00\times 10^{-6}212.0±4.9212.0\pm 4.9
10−410^{-4}6.4977−18.3871-18.3871(1.387±0.124)×10−5(1.387\pm 0.124)\times 10^{-5}1.00×10−81.00\times 10^{-8}1384±1241384\pm 124
Appendix
Table 17: Monte Carlo estimate and dependence for the compound hot-dry event at matched marginal quantiles. pp represents the compound hot-dry event at matched marginal quantiles: at each level pp both thresholds are the pool’s own pp -quantiles, so pindep=p2p_{\rm indep}=p^{2} exactly; quoted errors are binomial, and the 10−410^{-4} row rests on 125125 joint crossers.
 Factorised rewardJoint Φ2\Phi_{2} reward
marginal ppp^/p⋆\hat{p}/p^{\star}CV ↓\downarrowd^/d⋆\hat{d}/d^{\star}p^/p⋆\hat{p}/p^{\star}CV ↓\downarrowd^/d⋆\hat{d}/d^{\star}
10−210^{-2}0.990±0.024\textbf{0.990}\pm\textbf{0.024}0.0970.0971.07±0.031.07\pm 0.030.965±0.0160.965\pm 0.0160.0681.04±0.02\textbf{1.04}\pm\textbf{0.02}
10−310^{-3}0.952±0.0310.952\pm 0.0310.1320.1320.98±0.030.98\pm 0.030.990±0.028\textbf{0.990}\pm\textbf{0.028}0.1131.02±0.031.02\pm 0.03
10−410^{-4}1.124±0.0871.124\pm 0.0870.3100.3101.12±0.08\textbf{1.12}\pm\textbf{0.08}1.044±0.021\textbf{1.044}\pm\textbf{0.021}0.0801.06±0.06\textbf{1.06}\pm\textbf{0.06}
Appendix
Table 18: Comparison of two reward estimates for compound hot & dry events. The difference between the two rewards is most prominent at the rarest level. Notably the 10−410^{-4} reference itself carries an 8.9%8.9\% standard error. Setup remains identical to Table 1 , MC reference: Table 17 .
Figure 18: Conditional mean fields for the compound event of Table 1 , over the PNW box in July at GMST +1.5+1.5 K. Top: conditioned on the hot constraint alone. Bottom: conditioned jointly on hot & dry. Left: near-surface temperature. Right: relative humidity. Both constraints are at marginal 10−310^{-3} ; the box is the region of interest. Guiding on temperature alone reduces the humidity across a region.
ΔT\Delta T (K)p^\hat{p}s.e. ↓\downarrowCV ↓\downarrowp⋆p^{\star}p^/p⋆\hat{p}/p^{\star}
1.59.79 ⁣× ⁣10−49.79\!\times\!10^{-4}±0.18 ⁣× ⁣10−4\pm 0.18\!\times\!10^{-4}0.0731.000 ⁣× ⁣10−31.000\!\times\!10^{-3}0.9780.978
2.01.51 ⁣× ⁣10−21.51\!\times\!10^{-2}±0.03 ⁣× ⁣10−2\pm 0.03\!\times\!10^{-2}0.0731.485 ⁣× ⁣10−21.485\!\times\!10^{-2}1.0161.016
3.03.93 ⁣× ⁣10−13.93\!\times\!10^{-1}±0.02 ⁣× ⁣10−1\pm 0.02\!\times\!10^{-1}0.0253.932 ⁣× ⁣10−13.932\!\times\!10^{-1}1.0001.000
4.09.51 ⁣× ⁣10−19.51\!\times\!10^{-1}±0.02 ⁣× ⁣10−1\pm 0.02\!\times\!10^{-1}0.0069.503 ⁣× ⁣10−19.503\!\times\!10^{-1}1.0011.001
Appendix
Table 19: Warming sweep at the frozen event τ=5.81\tau=5.81 K over the region computed with 16 seeds per point. p∗p^{*} is obtained from a large Monte Carlo references, which exist only at ΔT≥1.5\Delta T\geq 1.5 K (relative errors 1.1%1.1\% , 3.2%3.2\% , 0.7%0.7\% , 0.13%0.13\% ) and the guided estimates match all four.
 DireSMCCrude MC 
ΔT\Delta T (K)p^\hat{p}CVp^MC\hat{p}_{\rm MC}CVZeroVar. Ratio
0.01.06 ⁣× ⁣10−61.06\!\times\!10^{-6}0.40200—16/16N/A
0.51.65 ⁣× ⁣10−61.65\!\times\!10^{-6}0.21400—16/16N/A
1.03.51 ⁣× ⁣10−53.51\!\times\!10^{-5}0.0921.30 ⁣× ⁣10−51.30\!\times\!10^{-5}4.00015/16257257
1.21.46 ⁣× ⁣10−41.46\!\times\!10^{-4}0.1401.04 ⁣× ⁣10−41.04\!\times\!10^{-4}1.78911/168383
1.59.79 ⁣× ⁣10−49.79\!\times\!10^{-4}0.0731.03 ⁣× ⁣10−31.03\!\times\!10^{-3}0.5940/167373
Appendix
Table 20: Warming sweep at the frozen event τ=5.81\tau=5.81 K over the region computed with 16 seeds per point. We compare against an unguided MC estimate at matched clock time allocating it 4.7×4.7\times particles, over 16 replicates. Zero counts the number of MC replicates that record no event. See Table 19 for DireSMC validated against a large MC pool. Var. Ratio = Var(p^MC)/Var(p^)\mathrm{Var}(\hat{p}_{\rm MC})/\mathrm{Var}(\hat{p}) .
Figure 19: Highest weight samples under the no resampling regime. We guide the sampler to the same τtas≥5.81\tau_{\texttt{tas}}\geq 5.81 K for a range of global average temperatures ΔT\Delta T . See Table 19 for the qualitative analysis.
 Threshold ExceedancesJoint hot & dry
p∗p^{*}hotdryfactorised rewardJoint reward
Per-sample cost4.7×4.7\times4.7×4.7\times4.7×4.7\times7.7×7.7\times
10−110^{-1}0.7×0.7\times1.0×1.0\times1.3×1.3\times1.2×1.2\times
10−210^{-2}8.7×8.7\times3.9×3.9\times7.0×7.0\times8.6×8.6\times
10−310^{-3}34×34\times40×40\times56×56\times47×47\times
10−410^{-4}311×311\times44×44\times156×156\times1413×1413\times
Appendix
Table 21: Net speed-up over crude Monte Carlo at matched CV for N=1024N=1024 guided draws against unguided draws. Computed using CVs obtained from Tables 15 , 1 and 18 . The dry 10−410^{-4} entry is low because of one collapsed run (App. H.2.1 ); without it the CV is 0.1240.124 , close to the other 10−410^{-4} cells.
   Hot ( tas )Dry ( hurs )
KnobDeployedAblatedp^/p∗\hat{p}/p^{\ast}CVp^/p∗\hat{p}/p^{\ast}CV
Deployed recipe  0.982±0.0190.982\pm 0.0190.0780.0780.996±0.0180.996\pm 0.0180.0720.072
Tempering schedule λ(σˉ)\lambda(\bar{\sigma})
Peak λmax⁡\lambda_{\max}1.251.251.01.0 (no overshoot)1.010±0.0451.010\pm 0.0450.1800.1800.984±0.0230.984\pm 0.0230.0950.095
Peak λmax⁡\lambda_{\max}1.251.251.51.50.984±0.0180.984\pm 0.0180.0710.0710.984±0.0220.984\pm 0.0220.0900.090
Peak λmax⁡\lambda_{\max}1.251.252.02.01.025±0.0511.025\pm 0.0510.1990.1990.942±0.0370.942\pm 0.0370.1590.159
Appendix
Table 22: One-knob ablation of the guidance schedules at p∗=10−3p^{\ast}=10^{-3} . We vary a single parameter reference in Table 14 ; the same 1616 seeds are used in every row. Collapsed runs are included: one hot run at λmax⁡=2.0\lambda_{\max}=2.0 , and two hot and two dry runs at δhi=0.25 SD\delta_{\mathrm{hi}}=0.25\,\mathrm{SD} .
guidanceKK\mathchoice{\vtop{\halign{#\cr\hfil\displaystyle p\hfil\cr\hfil\raisebox{4.30554pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}{\vtop{\halign{#\cr\hfil\textstyle p\hfil\cr\hfil\raisebox{4.30554pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}{\vtop{\halign{#\cr\hfil\scriptstyle p\hfil\cr\hfil\raisebox{3.01389pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}{\vtop{\halign{#\cr\hfil\scriptscriptstyle p\hfil\cr\hfil\raisebox{2.15277pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}/\mathchoice{\vtop{\halign{#\cr\hfil\displaystyle q\hfil\cr\hfil\raisebox{4.30554pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}{\vtop{\halign{#\cr\hfil\textstyle q\hfil\cr\hfil\raisebox{4.30554pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}{\vtop{\halign{#\cr\hfil\scriptstyle q\hfil\cr\hfil\raisebox{3.01389pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}}{\vtop{\halign{#\cr\hfil\scriptscriptstyle q\hfil\cr\hfil\raisebox{2.15277pt}{\scalebox{0.8}{\scriptstyle\rightharpoonup}}\hfil\cr}}} ( ( Manshausen et al., 2026 ) ) ↓\downarrowours ↓\downarrow
exact Doob50+155.7%±0.3+155.7\%{\scriptstyle\pm 0.3}−2.2%±0.7\mathbf{-2.2\%}{\scriptstyle\pm 0.7}
exact Doob100+65.3%±0.2+65.3\%{\scriptstyle\pm 0.2}−2.6%±0.4\mathbf{-2.6\%}{\scriptstyle\pm 0.4}
exact Doob200+31.7%±0.2+31.7\%{\scriptstyle\pm 0.2}−1.9%±0.3\mathbf{-1.9\%}{\scriptstyle\pm 0.3}
exact Doob500+14.7%±0.1+14.7\%{\scriptstyle\pm 0.1}−1.0%±0.2\mathbf{-1.0\%}{\scriptstyle\pm 0.2}
clipped Doob, ∼10%{\sim}10\% cross50+42.8%±4.2+42.8\%{\scriptstyle\pm 4.2}+2.7%±2.5\mathbf{+2.7\%}{\scriptstyle\pm 2.5}
clipped Doob, ∼10%{\sim}10\% cross100+21.2%±3.4+21.2\%{\scriptstyle\pm 3.4}−2.1%±1.9\mathbf{-2.1\%}{\scriptstyle\pm 1.9}
Appendix
Table 23: Signed percentage bias for the rare-event probability estimated against the analytical ground truth (obtained exactly without discretisation) for the 2D GMM threshold exceedance problem at τ=10\tau=10 . Our method tracks a lower bias at finite discretisation KK compared to ( Manshausen et al., 2026 ) . Positive %\% reflects an overestimate and are presented as a percentage against the true rare-event probability.
 score functionrr∇r\nabla rJVP
ours (Alg. 1 )K+2RMK{+}2RMK+RMK{+}RMK+RMK{+}RM–
odds ratio (Eq. 259 )3K3K2K2K2K2K2PK2PK
Appendix
Table 24: Network evaluations per particle. For a KK step denoising process, with RR rejuvenation steps and MM MCMC moves of Algorithm 1 . PP denotes the number of Hutchinson’s probes in Eq. ( 259 ).
  DireSMC
NNManshausen et al. (2026)base+ resampling+ MADM
1000.6310.9990.6110.321
2500.3850.5060.3940.184
5000.2590.3490.3030.134
7500.2060.2660.1680.082
10000.1670.2130.1660.079
Appendix
Table 25: CV of p^\hat{p} against particle count NN (Fig. 2 , left), at τ=10\tau=10 , K=1000K=1000 , with the learned score, over 24 seeds. Measured against model’s tail mass p∗p^{*} (Table 7 ).
   DireSMC
τ\taup⋆p^{\star}Manshausen et al. (2026)base+ resampling+ MADM
77.13×10−37.13\times 10^{-3}0.1290.0700.0930.060
81.38×10−31.38\times 10^{-3}0.1670.0920.0930.064
91.87×10−41.87\times 10^{-4}0.2210.1270.1590.090
101.72×10−51.72\times 10^{-5}0.3490.2130.1660.079
Appendix
Table 26: CV of p^\hat{p} relative to models threshold exceedance probability p⋆p^{\star} (Table 7 ). p^\hat{p} across rarities at equal compute (Fig. 2 , middle), with K=1000K=1000 , the learned score and 24 seeds. We use N=1000N=1000 and Manshausen et al. (2026) use n=N/3.5=286n=N/3.5=286 (Table 24 ). All variants of DireSMC outperform Manshausen et al. (2026) .
τ\taup⋆p^{\star}KKManshausen et al. (2026)DireSMC
76.88×10−36.88\times 10^{-3}501.466±0.0011.466\pm 0.0010.950±0.006\mathbf{0.950}\pm 0.006
1001.233±0.0011.233\pm 0.0010.972±0.003\mathbf{0.972}\pm 0.003
2001.128±0.0011.128\pm 0.0010.981±0.003\mathbf{0.981}\pm 0.003
5001.069±0.0011.069\pm 0.0010.990±0.002\mathbf{0.990}\pm 0.002
81.28×10−31.28\times 10^{-3}501.701±0.0021.701\pm 0.0020.957±0.005\mathbf{0.957}\pm 0.005
1001.336±0.0011.336\pm 0.0010.970±0.003\mathbf{0.970}\pm 0.003
Appendix
Table 27: Discretisation bias against the step count KK (Fig. 2 , right), shown as the mean p^/p⋆±\hat{p}/p^{\star}\pm s.e. over 24 seeds. Computed using analytical socres, Doob’s transition kernels and withut resampling at N=1000N=1000 . p⋆p^{\star} is pdata[E]p_{\mathrm{data}}[E] . Manshausen et al. (2026) are given the exact divergence.

Explore similar work

CardsList

  1. Robust Importance Sampling for Rare Events via Constrained Gaussian Mixtures

    Oct 5, 2026Paweł Lorek, Rafał Nowak, Rafał Topolnicki +2Rare EventsFinite-Sample

    University of Wrocław · Tooploox · TRAILS University of Warsaw +3

  2. Rare Event Estimation via Iterative Unalignment

    Sep 21, 2026Hanming Yang, Daksh Mittal, Jing Dong +1Rare EventsMonte Carlo

    Decision, Risk, and Operations Division, Columbia Business School

  3. Enhanced Diffusion Sampling: Efficient Rare Event Sampling and Free Energy Calculation with Diffusion Models

    Feb 18, 2026Yu Xie, Ludwig Winkler, Lixin Sun +10Diffusion SamplingMolecular Dynamics

    Microsoft Research AI for Science