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
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] of an event E 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] 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−3 to 10−5, achieving net speed-ups of 9× to 1413× over Monte Carlo.
Figures & tables
Figure 1: We steer particles ( ∙ ) from pref towards π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 ( × ) and form a rare event estimate, p0[E] .
Figure 2: Left: CV plotted against N for a p0[x1>10] problem at a fixed K . Rejuvenation tracks the lowest CV across all N . Middle: we compare our estimator against Manshausen et al. (2026) across rarities at a compute-normalised N . 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 K over the Pacific Northwest (PNW) target region. Left: average tas field over 107 samples. Remaining panels: the highest-weight crossing particle produced by the guided sampler at different exceedance thresholds τ and regional average temperature ϕ(⋅) of the sample. p⋆ are MC estimates; 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) .
p⋆
p^/p⋆
CV
d⋆
d^/d⋆
Speed Up ↑
3.14×10−3
0.97±
0.02
0.07
31
1.04±
0.02
9×
2.12×10−4
0.99±
0.03
0.11
212
1.02±
0.03
47×
1.39×10−5
1.04±
0.02
0.08
1384
1.06±
0.06
1413×
Table 1: Probability estimates of three compound hot and dry extreme events over the PNW region. The dependence ratio d is defined in the main text, 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 τ are given in App. H.1 . Mean ± standard error over 16 seeds.
Figure 4: Estimated p^0[E] versus global temperature anomaly ΔT for our method versus unguided MC at equal compute. Points are the mean over 16 runs, bars are ±1 std. Bars reaching the axis floor indicate levels where some MC runs record no event; for ΔT≤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
N
p^SNIS
Relative Bias
Z^0,SNIS
100
4.6035×10−1
+6,593%
4.873×10−1
1,000
3.5854×10−1
+5,113%
3.770×10−1
10,000
3.1680×10−1
+4,506%
3.367×10−1
25,000
1.9804×10−1
+2,779%
2.102×10−1
100
6.9250×10−3
+0.68%
(analytical Z0,δ )
Appendix
Table 2: Self-normalised importance sampling (SNIS) estimates vs. sample size ( N ).
ck
vk
ζk2=ck2/vk
Euler–Maruyama (VP/VE)
σtk2Δtk
σtk2Δtk
σtk2Δtk
exponential, VP ( ϵ^ frozen)
2σˉtk−1σˉtk(eΔℓk−1)
σˉtk−12(e2Δℓk−1)
4σˉtk2tanh(Δℓk/2)
exponential, VE ( x^ frozen)
(1−κk)σˉtk2
σˉtk−12(1−κk)
σˉtk2(e2Δℓk−1)
exponential, VE ( ϵ^ frozen)
2σˉtk−1σˉtk(eΔℓk−1)
σˉtk−12(e2Δℓk−1)
4σˉtk2tanh(Δℓk/2)
Appendix
Table 3: ζk2=ck2/vk for different integrators and noise schedules.
Figure 5: Soft reward and its guidance for a threshold event E7={x:ϕ(x)≥7} shown for three smoothing levels δ . Left. The soft indicator relaxes the hard indicator. Right. The guidance ∇xr0,δ(x) points toward the event set (shaded) and decays inside it, a smaller δ gives a stronger guidance.
Figure 6: Base distribution p0(x) of the correlated 2D GMM (log-density, so all modes are visible) with the analytic x1 and x2 marginals. The two right-hand modes ( i=5,6 ) form the tilted “X” of opposite-sign correlation that gives the rare set a bi-modal conditional.
i
wi
μi
Σi
1
0.10
(−6,−2)
[0.360.300.300.81]
2
0.18
(−4,4)
[0.81−0.25−0.250.49]
3
0.14
(−1,−5)
[0.490.200.200.81]
4
0.22
(−10,3)
[0.040.040.040.09]
5
0.10
(4,3)
[2.401.851.852.00]
6
0.16
(4,−3)
[2.40−1.85−1.852.00]
Appendix
Table 4: Parameters of the correlated 2D GMM. The right-hand modes ( i=5,6,7 ) carry the rare-event tail; modes 5 and 6 have correlation ρ≈+0.85 and −0.85 respectively.
τ
pdata[E]
pmodel (MC)
pmodel/pdata[E]
rel. s.e.
path integral
ratio
7
6.88×10−3
7.13×10−3
1.037
0.08%
1.075
1.037
8
1.28×10−3
1.38×10−3
1.081
0.18%
1.156
1.069
9
1.62×10−4
1.87×10−4
1.153
0.50%
1.236
1.072
10
1.40×10−5
1.72×10−5
1.229
1.64%
1.267
1.031
Appendix
Table 5: Tail bias of the learned 2D-GMM score. pmodel (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% across three decades of rarity. pdata[E] is the analytical probability.
K
EM bias
SEEDS-1 bias
∣EM∣/∣SEEDS∣
CV ↓
k^↓
ESSE↑
Uniform in t
10
+64.02% ± 11.31%
-15.98% ± 1.88%
4.01
0.338/0.109
0.85/0.46
0.032/0.102
15
+48.61% ± 8.39%
-16.23% ± 1.28%
3.00
0.277/0.075
0.70/0.45
0.063/0.189
20
+54.90% ± 16.50%
-15.13% ± 0.97%
3.63
0.522/0.056
0.68/0.43
0.099/0.248
30
+23.41% ± 1.65%
-13.62% ± 0.75%
1.72
0.065/0.043
0.56/0.40
0.176/0.375
50
+16.75% ± 1.38%
-8.36% ± 0.63%
2.00
0.058/0.034
0.55/0.45
0.315/0.465
Appendix
Table 6: Signed discretisation bias versus step count K for the two integrators, under exact Doob guidance, at different grids, uniform in t (top block) and uniform in the half-log SNR ℓt=log(αt/σˉt) (bottom block). Positive values are overestimates. Paired columns are EM/SEEDS-1 throughout. k^ is the Pareto Smoothed Importance Sampling shape estimate ( Vehtari et al., 2024 ) on the weights, and ESSE is the ESS of the final particles, measured after thresholding to τ≥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 K and evaluate the bias at different rarities. The figure on the right shows the number of discretisation required to get a <10% bias for different rarities, showing that the SDE requires less steps than the ODE integrator. Top row uses the uniform log-SNR Δℓ schedule and the bottom is uniform in time Δt .
Figure 8: Measured bias versus K 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) 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 ℓ . Error-bars represent standard error of the mean. Experiment is conducted the τ=10 threshold exceed problem. See Table 6 .
Figure 9: Measured bias versus K 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 108 unconditional MC draw of the exponential integrator at K as our reference. Hence, this measures how quickly the estimates converge to the reference with K .
τ
score
p⋆
p^/p⋆
CV ↓
KS x1
KS x2
EMD
7
analytical
6.88×10−3
0.980±0.012
0.060
0.88
1.33
1.52
8
analytical
1.28×10−3
0.966±0.012
0.061
1.13
1.10
1.06
9
analytical
1.62×10−4
0.974±0.016
0.083
1.12
1.00
0.97
10
analytical
1.40×10−5
0.954±0.013
0.069
0.94
1.05
1.46
7
learned
7.13×10−3
0.978±0.012
0.060
1.11
1.20
1.14
8
learned
1.38×10−3
0.985±0.013
0.064
0.97
1.00
0.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⋆ ; either pdata[E] (analytical) or an MC estimate of p0[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⋆ 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 .
Figure 10: Conditional-sample coverage on the 2D GMM for the analytical (a) and learned (b) scores. Blue is the path-weighted guided cloud ( 2048 samples; pooled over 8 seeds for the learned score); red is a MC draw from the conditional p(x∣x1≥τ) . Left: the joint cloud, with the threshold dashed. Centre, right: the x1 and x2 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} . Settings as Table 7 . Time runs right to left, from noise ( t=1 ) to data ( t=0 ).
τ
arm
CV ↓
KS x1
KS x2
EMD
7
no resampling
0.067
0.88
1.52
1.50
resampling
0.062
1.05
0.99
0.83
resampling + rejuv. ( M=3 )
0.060
0.88
1.33
1.52
8
no resampling
0.083
1.14
1.35
1.17
resampling
0.098
1.12
0.95
0.77
resampling + rejuv. ( M=3 )
0.061
1.13
1.10
1.06
Appendix
Table 8: Selection and rejuvenation on the 2D GMM under the analytical score, 24 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.
event
arm
p^/p⋆
CV ↓
KS 0
KS 1
EMD
x0>7
no resampling
0.952
0.070
0.74
1.00
1.09
resampling
0.988
0.093
0.94
1.25
1.47
resampling + rejuv. ( M=3 )
0.985
0.061
1.01
1.17
1.28
gated ℓ≲0.05
0.983
0.060
0.89
1.15
1.11
gated ℓ≲−0.21
0.978
0.060
1.11
1.20
1.14
gated ℓ≲−0.46
0.992
0.082
0.90
1.33
1.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 ℓ below which the rejuvenation move triggers; the adopted gate ℓ≲−0.21 is t≥0.30 . The no-resampling and resampling-only arms run at λ≡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 ρ for τ>10 threshold exceedance problem, for different calibrated step sizes ρ , 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=3 for the first two panels. Left & Middle Injected bias affects rarer samples more, with all biases decreasing as K→∞ . The two integrators have a different error profile. Right Bias plotted for the τ=10 threshold exceedance problem. We find that the error grows with M . 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.
event
variant
p⋆
p^/p⋆
CV ↓
KS 0
KS 1
EMD
[6,∞)×[1.5,∞)
no resampling
9.97×10−3
0.955
0.066
0.95
0.98
1.01
resampling
1.004
0.107
1.37
0.87
0.98
resampling + rejuv. ( M=3 )
1.000
0.070
0.89
0.91
0.96
gated ℓ≲−0.21
1.000
0.070
0.94
0.89
1.00
[6,∞)×[1.5,∞)∪[6,∞)×(−∞,−1.5]
no resampling
2.60×10−2
0.966
0.057
1.15
1.15
1.09
resampling
0.979
0.055
0.99
1.08
0.95
Appendix
Table 10: Box events on the learned score. Every region is drawn on the model density in Fig. 14 . p⋆ is the K=1000 chain’s own tail (Table 5 pool). We compare the estimate p^ and the shape metrics of the final particles, presented as a ratio against a noise floor.
schedule
base ↓
+ resampling ↓
+ MADM ↓
( Manshausen et al., 2026 ) ↓
analytic score
λ≡1
0.15
0.23
0.09
0.29
ton=1.0,p=2
0.28
0.32
0.08
0.68
ton=0.9,p=2 (deployed)
0.35
0.33
0.07
0.83
ton=0.7,p=2
0.62
0.40
0.10
1.26
ton=0.5,p=2
1.73
0.91
0.13
1.94
Appendix
Table 11: Schedule ablation at τ=10 , K=1000 , N=1000 : CV of p^ for the tempering schedule λ(t)=1−(t/ton)p . Bold marks the minimum of each column. We find that the no-rejuvenation variants converged on the same tempering schedule λt≡1 compared to the MADM variant. All schedules with rejuvenation yield a lower CV compared to without.
Short name
Description
Units
tas
2 m air temperature
K
hurs
2 m relative humidity
%
pr
Precipitation rate
mm day -1
sfcWind
10 m wind speed
m s -1
Appendix
Table 12: Emulated variables with their corresponding CMIP6 short names, descriptions and units.
Marginal
Upper tail
Lower tail
Variable
Unit
Mean
SD
10−3
10−4
10−5
10−3
10−4
10−5
tas
K
+2.577
1.023
5.813
6.498
7.096
−0.539
−1.174
−1.727
pr
mm/day
−0.250
0.450
1.698
2.311
2.842
−1.334
−1.517
−1.650
hurs
%
−3.275
4.220
10.519
13.485
16.054
−15.838
−18.387
−20.554
sfcWind
m/s
+0.193
0.186
0.773
0.895
1.004
−0.387
−0.508
−0.613
Appendix
Table 13: Tail thresholds for the four emulator variables, from N=9,011,200 unguided Monte Carlo draws over R at ΔT=+1.5 K in July. Upper-tail entries satisfy P(ϕ≥τ)=p , lower-tail entries P(ϕ≤τ)=p . Binomial relative standard errors are 1.1% , 3.3% and 10.5% at p=10−3,10−4,10−5 respectively. The guided runs use the same thresholds.
Figure 15: Discretisation bias for the integrator for the τ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=200 saturates the discretisation bias. We could not compare against a large K unconditional reference due to prohibitive costs. Bars represent ±1 standard error of the mean over 8 shared seeds across both integrators.
Symbol
Value
Schedules
Tempering onset
σˉon
40
Tempering peak
λmax
1.25 at σˉ=2
Tempering relaxation
σˉend
0.3
Smoothing width, transport
δhi
0.5SD
Smoothing width, floor
δlo
0.01SD
Appendix
Table 14: Deployed schedule and sampler configuration for the climate emulator experiments. σˉ is the marginal noise standard deviation. Selected once on the marginal 10−3 PNW hot and dry events and applied unchanged to other experiments. ‡ δ is quoted in units of the observable’s unconditional standard deviation.
Hot: ϕtas(x)≥τ
Dry: ϕhurs(x)≤τ
p^/p⋆
p^/p⋆
p⋆
Unguided
Guided
VR
p⋆
Unguided
Guided
VR
τ1
1.00e-3
0.90±0.45
0.98±0.08
34×
1.00e-3
1.03±0.44
1.00±0.07
37×
τ2
1.00e-4
1.30±1.67
1.01±0.08
411×
1.00e-4
0.52±0.93
1.02±0.22†
18×
τ3
1.01e-5
2.57±7.03
1.05±0.13
2722×
1.01e-5
2.57±7.03
1.06±0.12
3567×
Appendix
Table 15: Marginal extreme-event probabilities over the PNW target region in July at ΔT=+1.5 K, for hot ( ϕtas≥τ ) and dry ( ϕhurs≤τ ) events. For each variable, τ1,τ2,τ3 are chosen so that the event is successively rarer (exact values in App. H ). Ground-truth p⋆ are computed with ∼107 MC samples. Unguided Monte Carlo and our guided method use the same computational budget: each unguided estimate uses 4813 draws ( 4.7× the 1024 guided particles, Table 21 ), over 16 disjoint replicates. Estimated probabilities, p^ are reported as a ratio p^/p∗ with mean ±1 standard deviation over 16 seeds. The variance ratio (VR) is Var(p^unguided)/Var(p^guided) . The number of unguided replicates with no event is 1 (hot) and 0 (dry) of 16 at τ1 , 9 and 12 at τ2 , and 14 for both at τ3 , where the other two replicates record a single event each, so the VR at τ3 is a rough estimate. † Includes one collapsed run (App. H.2.1 ). The hot results are shown in Fig. 3 .
p⋆
arm
p^/p⋆
CV ↓
ESS/Nk^
W1
Speed-up ↑
tas
hurs
(95% CI)
Hot
10−3
no resampling
1.024±0.051
0.199
0.060.49
1.03
–
5.2× ( 2.8 – 12× )
adaptive resampling
0.982±0.019
0.078
–
1.54
–
34× ( 20 – 60× )
+ MADM moves
0.954±0.020
0.082
–
1.40
–
27× ( 17 – 44× )
10−4
no resampling
1.020±0.049
0.191
0.030.57
1.07
–
56× ( 32 – 94× )
adaptive resampling
1.010±0.021
0.082
–
0.95
–
311× ( 129 – 616× )
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 ± one standard error over 16 seeds. W1 is a weighted distance between the ϕ(⋅) of generated samples and draws from a reference, scaled against a null floor. The relative error of p⋆ is 1.1/3.3/10.5% at the three marginal rarities and 0.6/2.3/8.9% at the three joint rarities. CV is estimated across p^ estimates; k^ is estimated per seed and we report the median over the 16 seeds. † Includes one collapsed run; without it p^/p⋆=0.968 and CV =0.124 . 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 τ . Guided samples based on the tas variable alone. Same setup as reference Table 15 . p denotes p0[E] and ϕ(⋅) is the regional average temperature.
Figure 17: The same fields as Fig. 16 on the global domain.
marginal p
τtas (K)
τhurs (%)
pjoint
pindep
d⋆
10−1
3.8908
−8.6085
(5.028±0.007)×10−2
1.00×10−2
5.03±0.01
10−2
4.9907
−12.7822
(3.138±0.019)×10−3
1.00×10−4
31.4±0.2
10−3
5.8128
−15.8382
(2.121±0.049)×10−4
1.00×10−6
212.0±4.9
10−4
6.4977
−18.3871
(1.387±0.124)×10−5
1.00×10−8
1384±124
Appendix
Table 17: Monte Carlo estimate and dependence for the compound hot-dry event at matched marginal quantiles. p represents the compound hot-dry event at matched marginal quantiles: at each level p both thresholds are the pool’s own p -quantiles, so pindep=p2 exactly; quoted errors are binomial, and the 10−4 row rests on 125 joint crossers.
Factorised reward
Joint Φ2 reward
marginal p
p^/p⋆
CV ↓
d^/d⋆
p^/p⋆
CV ↓
d^/d⋆
10−2
0.990±0.024
0.097
1.07±0.03
0.965±0.016
0.068
1.04±0.02
10−3
0.952±0.031
0.132
0.98±0.03
0.990±0.028
0.113
1.02±0.03
10−4
1.124±0.087
0.310
1.12±0.08
1.044±0.021
0.080
1.06±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−4 reference itself carries an 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 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−3 ; the box is the region of interest. Guiding on temperature alone reduces the humidity across a region.
ΔT (K)
p^
s.e. ↓
CV ↓
p⋆
p^/p⋆
1.5
9.79×10−4
±0.18×10−4
0.073
1.000×10−3
0.978
2.0
1.51×10−2
±0.03×10−2
0.073
1.485×10−2
1.016
3.0
3.93×10−1
±0.02×10−1
0.025
3.932×10−1
1.000
4.0
9.51×10−1
±0.02×10−1
0.006
9.503×10−1
1.001
Appendix
Table 19: Warming sweep at the frozen event τ=5.81 K over the region computed with 16 seeds per point. p∗ is obtained from a large Monte Carlo references, which exist only at ΔT≥1.5 K (relative errors 1.1% , 3.2% , 0.7% , 0.13% ) and the guided estimates match all four.
DireSMC
Crude MC
ΔT (K)
p^
CV
p^MC
CV
Zero
Var. Ratio
0.0
1.06×10−6
0.402
0
—
16/16
N/A
0.5
1.65×10−6
0.214
0
—
16/16
N/A
1.0
3.51×10−5
0.092
1.30×10−5
4.000
15/16
257
1.2
1.46×10−4
0.140
1.04×10−4
1.789
11/16
83
1.5
9.79×10−4
0.073
1.03×10−3
0.594
0/16
73
Appendix
Table 20: Warming sweep at the frozen event τ=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× 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^) .
Figure 19: Highest weight samples under the no resampling regime. We guide the sampler to the same τtas≥5.81 K for a range of global average temperatures ΔT . See Table 19 for the qualitative analysis.
Threshold Exceedances
Joint hot & dry
p∗
hot
dry
factorised reward
Joint reward
Per-sample cost
4.7×
4.7×
4.7×
7.7×
10−1
0.7×
1.0×
1.3×
1.2×
10−2
8.7×
3.9×
7.0×
8.6×
10−3
34×
40×
56×
47×
10−4
311×
44×
156×
1413×
Appendix
Table 21: Net speed-up over crude Monte Carlo at matched CV for N=1024 guided draws against unguided draws. Computed using CVs obtained from Tables 15 , 1 and 18 . The dry 10−4 entry is low because of one collapsed run (App. H.2.1 ); without it the CV is 0.124 , close to the other 10−4 cells.
Hot ( tas )
Dry ( hurs )
Knob
Deployed
Ablated
p^/p∗
CV
p^/p∗
CV
Deployed recipe
0.982±0.019
0.078
0.996±0.018
0.072
Tempering schedule λ(σˉ)
Peak λmax
1.25
1.0 (no overshoot)
1.010±0.045
0.180
0.984±0.023
0.095
Peak λmax
1.25
1.5
0.984±0.018
0.071
0.984±0.022
0.090
Peak λmax
1.25
2.0
1.025±0.051
0.199
0.942±0.037
0.159
Appendix
Table 22: One-knob ablation of the guidance schedules at p∗=10−3 . We vary a single parameter reference in Table 14 ; the same 16 seeds are used in every row. Collapsed runs are included: one hot run at λmax=2.0 , and two hot and two dry runs at δhi=0.25SD .
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 . Our method tracks a lower bias at finite discretisation K compared to ( Manshausen et al., 2026 ) . Positive % reflects an overestimate and are presented as a percentage against the true rare-event probability.
score function
r
∇r
JVP
ours (Alg. 1 )
K+2RM
K+RM
K+RM
–
odds ratio (Eq. 259 )
3K
2K
2K
2PK
Appendix
Table 24: Network evaluations per particle. For a K step denoising process, with R rejuvenation steps and M MCMC moves of Algorithm 1 . P denotes the number of Hutchinson’s probes in Eq. ( 259 ).
DireSMC
N
Manshausen et al. (2026)
base
+ resampling
+ MADM
100
0.631
0.999
0.611
0.321
250
0.385
0.506
0.394
0.184
500
0.259
0.349
0.303
0.134
750
0.206
0.266
0.168
0.082
1000
0.167
0.213
0.166
0.079
Appendix
Table 25: CV of p^ against particle count N (Fig. 2 , left), at τ=10 , K=1000 , with the learned score, over 24 seeds. Measured against model’s tail mass p∗ (Table 7 ).
DireSMC
τ
p⋆
Manshausen et al. (2026)
base
+ resampling
+ MADM
7
7.13×10−3
0.129
0.070
0.093
0.060
8
1.38×10−3
0.167
0.092
0.093
0.064
9
1.87×10−4
0.221
0.127
0.159
0.090
10
1.72×10−5
0.349
0.213
0.166
0.079
Appendix
Table 26: CV of p^ relative to models threshold exceedance probability p⋆ (Table 7 ). p^ across rarities at equal compute (Fig. 2 , middle), with K=1000 , the learned score and 24 seeds. We use N=1000 and Manshausen et al. (2026) use n=N/3.5=286 (Table 24 ). All variants of DireSMC outperform Manshausen et al. (2026) .
τ
p⋆
K
Manshausen et al. (2026)
DireSMC
7
6.88×10−3
50
1.466±0.001
0.950±0.006
100
1.233±0.001
0.972±0.003
200
1.128±0.001
0.981±0.003
500
1.069±0.001
0.990±0.002
8
1.28×10−3
50
1.701±0.002
0.957±0.005
100
1.336±0.001
0.970±0.003
Appendix
Table 27: Discretisation bias against the step count K (Fig. 2 , right), shown as the mean p^/p⋆± s.e. over 24 seeds. Computed using analytical socres, Doob’s transition kernels and withut resampling at N=1000 . p⋆ is pdata[E] . Manshausen et al. (2026) are given the exact divergence.
We study estimating rare-event probabilities I=P(g(X)>γ) with X∼N(μ,Σ) and general g:Rd→R. We address this problem through importance sampling, and propose a framework that substantially improves efficiency and robustness over baselines such as crude Monte Carlo, adaptive cross-entropy, variational-inference-based methods (including reverse- and forward-KL approaches), as well as Safe-ICE, Subset Simulation, and Sequential Monte Carlo, drawing on ideas from both rare-event estimation and cross-entropy optimization. The key contribution has two parts: first, we separate the problem into coverage, to overcome the cold-start barrier, and fitting, to refine proposals once a meaningful signal is available; second, we constrain the final GMM proposal so that it has finite importance-sampling variance (since coverage alone is not sufficient -- without safeguards, importance sampling may still suffer from infinite variance). Together, these ingredients yield expressive proposals; finite variance does not by itself guarantee practical stability at a fixed sampling budget. Extensive experiments demonstrate substantial variance reduction, strong robustness across diverse benchmarks, and favorable cost--efficiency trade-offs, with the proposed approach often outperforming these baselines, particularly in high-dimensional and multimodal settings where competing methods frequently become unstable or fail. Our code is available at https://github.com/lorek/robust-cfi-is.
Paweł Lorek, Rafał Nowak, Rafał Topolnicki +2
University of Wrocław · Tooploox · TRAILS University of Warsaw +3
As agents are deployed with increased autonomy, even extremely rare events along their stochastic output trajectories can occur and prove catastrophic. Safe deployment therefore does not depend on whether these events can occur, but on how often they might. We study the problem of estimating the probability of rare events that arise from stochastic variation in the agent's own actions. Estimating this type of risk requires searching over the combinatorially vast space of trajectories. Naive Monte Carlo is computationally prohibitive in this regime, and constructing effective importance sampling (IS) proposals requires coordinated changes to a context-dependent chain of conditional distributions. We develop a new IS method that perturbs the original model's weights to construct the proposal. The proposal is itself a differentiably parameterized language model, enabling gradient-based search over weight space. We formulate an objective that combines a differentiable surrogate for event amplification and an adaptive regularization scheme that dynamically balances amplification against estimator stability. We evaluate our approach on ∼120M and ∼2.6B models across three event families spanning 300+ rare events as rare as 10−9, with reference probabilities computed with <10% relative standard error. In our most verifiable settings, we observe that our IS estimator achieves over 800× compute-weighted efficiency gains over naive Monte Carlo for events with probabilities lower than 10−7. Our implementation is available at https://github.com/namkoong-lab/iterative-unalignment.
Hanming Yang, Daksh Mittal, Jing Dong +1
Decision, Risk, and Operations Division, Columbia Business School
The rare-event sampling problem has long been the central limiting factor in molecular dynamics (MD), especially in biomolecular simulation. Recently, diffusion models such as BioEmu have emerged as powerful equilibrium samplers that generate independent samples from complex molecular distributions, eliminating the cost of sampling rare transition events. However, a sampling problem remains when computing observables that rely on states which are rare in equilibrium, for example folding free energies. Here, we introduce enhanced diffusion sampling, enabling efficient exploration of rare-event regions while preserving unbiased thermodynamic estimators. The key idea is to perform quantitatively accurate steering protocols to generate biased ensembles and subsequently recover equilibrium statistics via exact reweighting. We instantiate our framework in three algorithms: UmbrellaDiff (umbrella sampling with diffusion models), MetaDiff (a batchwise analogue for metadynamics), and ΔG-Diff (free-energy differences via tilted ensembles). Across toy systems, protein folding landscapes and folding free energies, our methods achieve fast, accurate, and scalable estimation of equilibrium properties within GPU-minutes to hours per system-closing the rare-event sampling gap that remained after the advent of diffusion-model equilibrium samplers.