Chapter 5 Monte Carlo Optimization & Rare-Event Simulation

5.1 Sensitivity Analysis (Score Function Method)

For a parameterized expectation $\alpha(\theta) = E_\theta[h(X)] = \int h(x) f(x; \theta) \, dx$, we want to estimate the gradient $\nabla_\theta \alpha(\theta)$. Using the Score Function / Likelihood Ratio Method: \[\nabla_\theta \alpha(\theta) = E_\theta\left[ h(X) \nabla_\theta \log f(X; \theta) \right]\] The term $S(X; \theta) = \nabla_\theta \log f(X; \theta)$ is the score function.

Suppose you want to know how the average height of a tree species $\alpha(\theta)$ changes as you adjust the soil moisture level $\theta$. The naive way requires planting trees at many different levels of moisture and measuring the difference. The Score Function method is a mathematical shortcut: instead of changing the moisture and replanting, you measure the height of the current trees $h(X)$ and multiply it by a “score” $S(X; \theta)$ that tracks how sensitive each tree’s likelihood of growth is to moisture. Averaging this product gives the exact gradient (sensitivity) without needing to run simulations at new parameters.

Step-by-Step Proof:
1. Write as Integral: By definition of the expectation: \[\alpha(\theta) = \int h(x) f(x; \theta) \, dx\]
2. Differentiate Under the Integral Sign: We take the gradient with respect to $\theta$ on both sides: \[\nabla_\theta \alpha(\theta) = \nabla_\theta \int h(x) f(x; \theta) \, dx = \int h(x) \nabla_\theta f(x; \theta) \, dx\]
3. Divide and Multiply by the Density: Assuming $f(x; \theta) > 0$ on the support, we multiply and divide the integrand by $f(x; \theta)$: \[\nabla_\theta \alpha(\theta) = \int h(x) \frac{\nabla_\theta f(x; \theta)}{f(x; \theta)} f(x; \theta) \, dx\]
4. Apply logarithmic derivative identity: Since $\nabla_\theta \log f(x; \theta) = \frac{\nabla_\theta f(x; \theta)}{f(x; \theta)}$, we substitute: \[\nabla_\theta \alpha(\theta) = \int h(x) \nabla_\theta \log f(x; \theta) f(x; \theta) \, dx\]
5. Rewrite as Expectation: This integral is exactly the expectation of $h(X) \nabla_\theta \log f(X; \theta)$ under $X \sim f(x; \theta)$: \[\nabla_\theta \alpha(\theta) = E_\theta\left[ h(X) \nabla_\theta \log f(X; \theta) \right] \quad \blacksquare\]
Practical Applications: Used in reinforcement learning (Policy Gradient methods like REINFORCE), finance (calculating “Greeks” or option price sensitivities), and engineering design optimization.
Manual Walkthrough (By-Hand Simulation): Let $X \sim N(\theta, 1)$. We want to estimate the gradient of $\alpha(\theta) = E_\theta[X^3]$ at $\theta = 1.0$. The analytical value is $\frac{d}{d\theta} E_\theta[X^3] = \frac{d}{d\theta}(\theta^3 + 3\theta) = 3\theta^2 + 3$, which at $\theta = 1.0$ is $6.0$. The score function is: \[S(x; \theta) = \frac{d}{d\theta} \log f(x; \theta) = \frac{d}{d\theta} \left[ -\frac{1}{2}\log(2\pi) - \frac{(x - \theta)^2}{2} \right] = x - \theta\] Suppose we simulate a single sample $x = 2.0$ from $N(1.0, 1)$. The single-sample score function estimate of the gradient is: \[\text{grad}_{\text{est}} = h(x) S(x; \theta) = x^3 (x - \theta) = 2.0^3 (2.0 - 1.0) = 8.0 \times 1.0 = 8.0\]

The following flowchart visualizes the process of estimating parameter sensitivity using the Score Function method.

## 
## Attaching package: 'DiagrammeR'

## The following object is masked _by_ '.GlobalEnv':
## 
##     render_graph

Figure 5.1: Gradient Estimation via Score Function Method

5.1.1 Score Function Code Verification

We estimate the gradient of $E_\theta[X^3]$ at $\theta = 1.0$ using 10,000 samples and compare it with the analytical gradient ($6.0$):

n_samples <- 10000
theta_val <- 1.0

# Simulate from target density at theta_val
x_sims <- rnorm(n_samples, mean = theta_val, sd = 1)
# Estimate gradient
sf_gradients <- x_sims^3 * (x_sims - theta_val)
emp_gradient <- mean(sf_gradients)
emp_grad_se <- sd(sf_gradients) / sqrt(n_samples)

The empirical gradient estimate is 5.9001 (standard error: 0.1690), which is within statistical tolerance of the true analytical gradient value $6.0000$.

5.2 Robbins-Monro Stochastic Approximation

To find the root $\theta^*$ of an expectation function $M(\theta) = E[Y(\theta)] = 0$, where we only observe noisy measurements $y_n(\theta_n) = M(\theta_n) + \epsilon_n$, we use the Robbins-Monro algorithm: \[\theta_{n+1} = \theta_n - a_n y_n(\theta_n)\]

Imagine you are trying to balance a scale to find a weight $\theta^*$ where the reading is exactly zero. However, your scale is highly sensitive to wind, so every time you measure, the reading bounces around. The Robbins-Monro algorithm is a rule for adjusting the weight step-by-step: if the scale reads too high, you subtract some weight; if too low, you add weight. Crucially, as time goes on, you make your adjustments smaller and smaller (decaying step size $a_n$). This allows you to explore the scale quickly at first, but prevents wind noise from throwing off your final balance as you get closer to the true weight.

Convergence Conditions: The algorithm converges to the true root $\theta^*$ almost surely if:
1. $\sum_{n=1}^\infty a_n = \infty$ (ensures step sizes are large enough to reach the root from any distance).
2. $\sum_{n=1}^\infty a_n^2 < \infty$ (ensures the steps become small enough to cancel out accumulated noise). An example sequence is $a_n = c/n$.
Newton-Raphson Analogy: To gain intuitive understanding, compare Robbins-Monro with the classical deterministic Newton-Raphson root-finding algorithm:
- In Newton-Raphson, we update $\theta_{n+1} = \theta_n - [M'(\theta_n)]^{-1} M(\theta_n)$.
- In a stochastic setting, we cannot evaluate $M(\theta_n)$ directly (we only observe noisy $y_n(\theta_n)$), and the derivative $M'(\theta_n)$ is unknown.
- Robbins-Monro replaces the derivative term $[M'(\theta_n)]^{-1}$ with a deterministic step size sequence $a_n$ (which slowly decays to suppress observation noise $\epsilon_n$), and replaces $M(\theta_n)$ with the noisy sample $y_n(\theta_n)$.
Practical Applications: Used in online learning, adaptive parameter tuning in neural networks, and stochastic optimization where data arrives sequentially.
Manual Walkthrough (By-Hand Simulation): We want to find the root of $3\theta - 9 = 0$ (root $\theta^* = 3$) using step sizes $a_n = 1/n$, starting from initial guess $\theta_1 = 0$.
- Iteration 1 ($n=1$): Suppose we draw noisy observation $y_1 = 3\theta_1 - 9 + \epsilon_1$. Let $\epsilon_1 = -0.5 \implies y_1 = 3(0) - 9 - 0.5 = -9.5$. Update: \[\theta_2 = \theta_1 - a_1 y_1 = 0 - 1 \times (-9.5) = 9.5\]
- Iteration 2 ($n=2$): Suppose we draw noisy observation $y_2 = 3\theta_2 - 9 + \epsilon_2$. Let $\epsilon_2 = 1.2 \implies y_2 = 3(9.5) - 9 + 1.2 = 20.7$. Update: \[\theta_3 = \theta_2 - a_2 y_2 = 9.5 - 0.5 \times 20.7 = 9.5 - 10.35 = -0.85\]

The following flowchart visualizes the iterative step-by-step root-finding process of the Robbins-Monro stochastic approximation algorithm.

Figure 5.2: Robbins-Monro Stochastic Root Finding Algorithm

5.2.1 Robbins-Monro Code Verification

We find the root of $M(\theta) = 3\theta - 9 = 0$ (analytical root: 3) where observations are corrupted by standard normal noise:

theta_est <- 0  # Initial guess
for (n in 1:1000) {
  y_obs <- 3 * theta_est - 9 + rnorm(1)
  a_n <- 1 / n
  theta_est <- theta_est - a_n * y_obs
}

The estimated root after 1000 iterations is 2.9936 (analytical root: 3.0000).

5.3 Splitting Methods (Rare Event Simulation)

To estimate the probability $\ell = P(X \in A)$ of a rare event $A$, direct Monte Carlo requires an impractically large sample size. Splitting methods (such as fixed factor or fixed effort) decompose the rare event path into a sequence of nested intermediate events $E_1 \supset E_2 \dots \supset E_m = A$: \[\ell = P(E_1) \prod_{k=2}^m P(E_k \mid E_{k-1})\]

Suppose you want to estimate the probability of a balls-in-a-maze toy successfully reaching a tiny exit gate $A$ at the bottom. If you drop balls from the top, only 1 in a million might reach the gate (a rare event). In a splitting method, you define intermediate checkpoint lines (levels $E_1, E_2$). When a ball successfully crosses a checkpoint, you freeze its state and clone (split) it into multiple balls, releasing them to continue the maze. If a ball fails and rolls backwards, you discard (cull) it. By multiplying the fraction of balls that pass each checkpoint, you get the exact success rate without having to run millions of failed top-level trials.

Practical Applications: Used in nuclear safety analysis (modeling radiation leakage), telecommunications (estimating packet loss rates), and financial risk modeling.
Manual Walkthrough (By-Hand Simulation): We trace the Fixed-Effort multi-level splitting algorithm for a 3-step random walk starting at 0 with standard normal increments $Z_t \sim N(0, 1^2)$. We want to estimate the rare probability $P(X_3 \ge 3.0)$ using intermediate levels $L_1 = 1.5$ and $L_2 = 3.0$, with a constant effort of $N = 2$ paths:
- Level 1: Run $N = 2$ independent paths of length 3:
  - Path 1: $(0, 1.0, 1.2, 1.4) \implies$ maximum value is $1.4 < 1.5$. Fail (terminated).
  - Path 2: $(0, 0.8, 1.6, 2.2) \implies$ maximum value is $2.2 \ge 1.5$. Success! Crossing occurs at step $\tau_1 = 2$.
  - Level 1 probability estimate: $p_1 = 1 / 2 = 0.5$.
- Level 2: We resample $N = 2$ starting paths from the successful Path 2 up to crossing time $\tau_1 = 2$: the prefix is $(0, 0.8, 1.6)$. We simulate the remaining step 3 for both paths:
  - Path 2a: starts from $X_2 = 1.6$, draws increment $1.5 \implies X_3 = 3.1 \ge 3.0$. Success!
  - Path 2b: starts from $X_2 = 1.6$, draws increment $-0.2 \implies X_3 = 1.4 < 3.0$. Fail (terminated).
  - Level 2 conditional probability estimate: $p_2 = 1 / 2 = 0.5$.
- Total Probability Estimate: \[\hat{p} = p_1 \times p_2 = 0.5 \times 0.5 = 0.25\]

The following tree diagram visualizes the branching and culling of the trajectories at each level checkpoint for the 2-level splitting process.

Figure 5.3: Fixed-Factor Rare-Event Splitting Path Tree

5.3.1 Rare-Event Splitting Code Verification

We estimate the rare joint probability $p = P(X_1 \ge 1.2, X_2 \ge 2.8, X_3 \ge 4.5, X_4 \ge 6.2, X_5 \ge 8.0)$ for a 5-step random walk $X_t = X_{t-1} + Z_t$ with increments $Z_i \sim N(0, 1.5^2)$ using the time-step checkpoints Fixed-Effort splitting algorithm from helpers.R. The analytical probability of this joint event is approximately $0.00365$. We set the levels to c(1.2, 2.8, 4.5, 6.2, 8.0) and compare the splitting estimate and standard error against standard Monte Carlo:

# Run Fixed-Effort Splitting
nsims_split <- 1000
split_estimates <- replicate(50, run_rare_event_splitting(N = nsims_split, levels = c(1.2, 2.8, 4.5, 6.2, 8.0)))
split_mean <- mean(split_estimates)
split_se <- sd(split_estimates)

# Compare with Standard Monte Carlo with same total sample size (5000 samples per run)
mc_estimates <- replicate(50, {
  mc_draws <- replicate(5000, {
    x <- 0
    path <- numeric(5)
    for (t in 1:5) {
      x <- x + rnorm(1, mean = 0, sd = 1.5)
      path[t] <- x
    }
    all(path >= c(1.2, 2.8, 4.5, 6.2, 8.0))
  })
  mean(mc_draws)
})
mc_mean <- mean(mc_estimates)
mc_se <- sd(mc_estimates)

The splitting estimate of the rare joint probability is 0.00373 (standard error: 0.00057), while the standard Monte Carlo estimate is 0.00357 (standard error: 0.00102), confirming that both converge to the same value (approx. 0.0037) and that multi-level splitting achieves a significantly lower standard error for rare-event estimation.

5.4 Adaptive MCMC (Adaptive Metropolis)

Adaptive MCMC adjusts proposal distributions dynamically during chain execution to target optimal criteria. For random walk proposals, Roberts-Gelman-Gilks proved that the optimal acceptance rate is $0.234$ in high dimensions:

Online Covariance Adaptation: At iteration $t$, we update the log of the proposal standard deviation based on the history: \[\log \sigma_{t} = \log \sigma_{t-1} + \gamma_t (\alpha_t - 0.234)\] where $\gamma_t = t^{-\kappa}$ is a vanishing step size and $\alpha_t$ is the acceptance rate.
Manual Walkthrough (By-Hand Simulation): We trace the online proposal variance adaptation over 2 sequential adaptation windows using step sizes $\gamma_t = \frac{1}{\sqrt{t}}$ starting from an initial proposal standard deviation $\sigma_0 = 1.0$ ($\log \sigma_0 = 0.0$).
- Adaptation Window 1 ($t = 1$, $\gamma_1 = 1.0$): Suppose we run the sampler for a window of $K$ steps and observe an empirical acceptance rate of $\alpha_1 = 0.40$ (which exceeds the target rate of $0.234$). We update the log-scale standard deviation: \[\log \sigma_1 = \log \sigma_0 + \gamma_1 (\alpha_1 - 0.234) = 0.0 + 1.0 \times (0.40 - 0.234) = 0.166\] The new proposal standard deviation is: \[\sigma_1 = e^{0.166} \approx 1.1806\] Interpretation: Since the acceptance rate is too high, the algorithm increases the proposal standard deviation to make larger, more explorative moves.
- Adaptation Window 2 ($t = 2$, $\gamma_2 = 1/\sqrt{2} \approx 0.7071$): Suppose we run the sampler for the next $K$ steps using proposal standard deviation $\sigma_1 = 1.1806$ and observe an empirical acceptance rate of $\alpha_2 = 0.10$ (which is below the target rate of $0.234$). We update the log-scale standard deviation: \[\log \sigma_2 = \log \sigma_1 + \gamma_2 (\alpha_2 - 0.234) = 0.166 + 0.7071 \times (0.10 - 0.234) = 0.166 - 0.0948 = 0.0712\] The new proposal standard deviation is: \[\sigma_2 = e^{0.0712} \approx 1.0738\] Interpretation: Since the acceptance rate is too low, the algorithm decreases the proposal standard deviation to make smaller, more conservative moves.

The following flowchart visualizes the adaptive proposal standard deviation update cycle.

Figure 5.4: Adaptive Metropolis-Hastings Parameter Tuning Loop

5.4.1 Adaptive Metropolis Code Verification

We run the Adaptive Metropolis sampler from helpers.R targeting a standard Normal distribution and check that the acceptance rate converges to the optimal target:

# Target standard normal
normal_target <- function(x) dnorm(x, mean = 0, sd = 1)

# Run adaptive sampler
adapt_chain <- run_adaptive_metropolis(normal_target, init = 0.0, n_iter = 10000)
# Calculate acceptance rate in latter 50% of the chain (post-adaptation phase)
post_burn_chain <- adapt_chain[5001:10000]
transitions <- diff(post_burn_chain)
acc_rate_final <- mean(transitions != 0)

The final acceptance rate in the post-burn-in phase is 0.2386 (optimal target: 0.2340).

5.5 How to Report This in a Paper

When publishing optimization or rare-event simulation results:

Robbins-Monro Sequence: Report the exact step-size formula $a_n$ and the initial parameter vector $\theta_0$.
Intermediate Boundaries: For splitting methods, explicitly define the boundaries of the nested subsets $E_k$ and the replica factors $r_k$.
Adaptive MCMC Settings: Report the adaptation schedule, the target acceptance rate, and verify that the adaptation satisfies diminishing adaptation conditions.

Key Formulas Summary

Algorithm / Technique	Equation / Definition	Key Conditions / Parameters
Score Function Gradient	$\nabla_\theta E_\theta[h(X)] = E_\theta[h(X) \nabla_\theta \log f(X; \theta)]$	Assumes differentiability under the integral sign.
Robbins-Monro Updater	$\theta_{n+1} = \theta_n - a_n y_n(\theta_n)$	Decaying step size sequence $a_n$.
RM Convergence	$\sum_{n=1}^\infty a_n = \infty \quad \text{and} \quad \sum_{n=1}^\infty a_n^2 < \infty$	e.g. $a_n = c/n$.
Multi-Level Splitting	$\ell = P(E_1) \prod_{k=2}^m P(E_k \mid E_{k-1})$	Decomposes rare event $A = E_m$ into nested checkpoints.
Adaptive MCMC proposal	$Y \sim N(X_n, \sigma_n^2 \mathbf{I})$	$\log \sigma_{n+1} = \log \sigma_n + \gamma_n (\alpha_n - \alpha^)$, target $\alpha^ = 0.234$.
Diminishing Adaptation	$\lim_{n \to \infty} \gamma_n = 0$	Ensures transition kernels asymptotically stabilize.

Common Mistakes

Violating Robbins-Monro Convergence Criteria: Using a step-size sequence like $a_n = 1/n^2$ will cause the step size to decay too quickly. If your initial guess $\theta_1$ is far from the true root $\theta^*$, the sum of steps is bounded ($\sum a_n < \infty$), meaning the algorithm will freeze before ever reaching the root. Conversely, a sequence like $a_n = 1/\sqrt{n}$ decays too slowly ($\sum a_n^2 = \infty$), which fails to filter out observation noise, causing the estimator to oscillate indefinitely.

Failing to Track Path Weights in Splitting: During rare-event splitting, cloned (split) paths must be correctly weighted. If you split a path into $R$ replicas, each replica carries a relative probability weight of $1/R$. Forgetting to apply this weight factor will lead to an exponential overestimation of the rare-event probability.

Undetected Non-diminishing MCMC Adaptation: In adaptive MCMC, if the adaptation step size $\gamma_n$ does not decay to 0, the proposal variance will continue to adapt forever. This violates the ergodic theorem (since the chain is non-homogeneous and does not stabilize to a single Markov kernel), leading to samples that do not reflect the true target distribution.

5.6 Exercises

5.6.1 Sensitivity Analysis & Score Function

Derive the score function estimator for a Gamma distribution $\Gamma(k, \theta)$ with respect to the scale parameter $\theta$.
Prove that the expected value of the score function is always zero: $E_\theta[\nabla_\theta \log f(X; \theta)] = 0$.
Construct a manual walkthrough to estimate the sensitivity of $E[X^2]$ for $X \sim \text{Exp}(\lambda)$ at $\lambda = 2.0$ using a single uniform draw $u = 0.5$.

5.6.2 Robbins-Monro Stochastic Approximation

Write an R script to run the Robbins-Monro algorithm to find the root of $M(\theta) = \theta^3 - 8 = 0$ and verify convergence.
Show why the step size sequence $a_n = 1/n^2$ fails to satisfy the first convergence condition of the Robbins-Monro theorem.
Compare the convergence speed of Robbins-Monro using step sizes $a_n = 1/n$ versus $a_n = 1/\sqrt{n}$.

5.6.3 Rare-Event Simulation

Explain the difference between Fixed Factor and Fixed Effort splitting in terms of implementation complexity and variance characteristics.
Formulate an Importance Sampling proposal density to estimate $P(Z \ge 5)$ for $Z \sim N(0, 1)$ and prove that it has bounded relative error.
Build a manual walkthrough showing the splitting transitions of a 2-level splitting process for a simple 3-step random walk.

5.6.4 Adaptive MCMC

Prove why a constant adaptation of proposal variance without a vanishing step-size sequence ($\gamma_t \to 0$) violates the Markov property and can lead to a wrong stationary distribution.
Describe the Roberts-Gelman-Gilks optimal acceptance rate of $0.234$ and explain why it differs from the optimal rate of $0.44$ for a 1D target.

Chapter 5 covered:

Score Function / Likelihood Ratio method: Gradient estimation without re-simulation; $\nabla_\theta E_\theta[h(X)] = E_\theta[h(X) \nabla_\theta \log f(X; \theta)]$; applications in RL (REINFORCE) and finance (Greeks)
Robbins-Monro algorithm: Iterative root-finding under noisy observations; the theoretical ancestor of stochastic gradient descent; convergence requires $\sum a_n = \infty$ and $\sum a_n^2 < \infty$
Multi-level splitting: Decompose $P(A) = \prod_{k=1}^m P(A_k | A_{k-1})$ to estimate rare event probabilities using a series of intermediate checkpoints; far more efficient than crude Monte Carlo for probabilities $< 10^{-6}$
Adaptive MCMC: Automatically tune the proposal distribution during the chain run; caution: naive adaptation can destroy the Markov property — adaptation must satisfy diminishing adaptation conditions

Method	Problem Type	Key Idea
Score Function	Gradient of expectation w.r.t. $ heta$	Differentiate log-density, not the function
Robbins-Monro	Root-finding / optimization under noise	Noisy gradient steps with decreasing learning rate
Multi-level Splitting	Rare event probability $P(A) \ll 10^{-4}$	Decompose into conditional probabilities
Adaptive MCMC	Efficient MCMC in unknown geometry	Auto-tune proposal covariance from past samples

Congratulations — you have completed the Statistical Research Methodology course notes! The five chapters have taken you from research design (Ch. 1), through simulation fundamentals (Ch. 2), Monte Carlo integration and variance reduction (Ch. 3), Bayesian MCMC sampling (Ch. 4), to stochastic optimization and rare-event estimation (Ch. 5). These tools form the computational backbone of modern statistical research.

Exercises

Tier 1 — Conceptual

The Robbins-Monro algorithm is said to “converge almost surely” under its step-size conditions. What does “almost surely” mean in the probability-theoretic sense? Why is this a stronger statement than convergence in probability?
Multi-level splitting requires defining a sequence of nested “intermediate” sets $A_1 \supset A_2 \supset \cdots \supset A_m = A$. What happens if you choose only two levels? What happens if you choose too many?
Adaptive MCMC adapts the proposal covariance using past samples. Why must the adaptation eventually stop (or diminish to zero) for the sampler to remain valid? What theorem is violated if it does not?

Tier 2 — Applied

Implement the Robbins-Monro algorithm in R to find the root of $M( heta) = heta^3 - 2$ (i.e., $ heta^* = 2^{1/3}$) using noisy observations $Y_n = M( heta_n) + \epsilon_n$ with $\epsilon_n \sim N(0, 1)$. Use step sizes $a_n = 1/n$, start at $ heta_0 = 0$, and run for 500 iterations. Plot the iterates $\{ heta_n\}$ and report the final estimate.
Use multi-level splitting to estimate $P(Z > 5)$ where $Z \sim N(0,1)$ (true value $pprox 2.87 imes 10^{-7}$). Use intermediate levels at $z = 2, 3, 4, 5$ and $N = 1{,}000$ particles at each level.

Tier 3 — Challenge

Stochastic Gradient Descent for Logistic Regression: Using the course_dataset from SCNPIR (or the built-in spam dataset from kernlab), implement mini-batch SGD to fit a logistic regression model. Use a learning rate schedule $\eta_t = \eta_0 / (1 + t)^{0.6}$ with $\eta_0 = 0.1$ and mini-batches of size 32. Plot the cross-entropy loss vs. iteration and compare the final coefficients to those from glm().

Solutions to Tier 1 and 2: see Appendix B.

Algorithm / Technique	Equation / Definition	Key Conditions / Parameters
Score Function Gradient	\(\nabla_\theta E_\theta[h(X)] = E_\theta[h(X) \nabla_\theta \log f(X; \theta)]\)	Assumes differentiability under the integral sign.
Robbins-Monro Updater	\(\theta_{n+1} = \theta_n - a_n y_n(\theta_n)\)	Decaying step size sequence \(a_n\).
RM Convergence	\(\sum_{n=1}^\infty a_n = \infty \quad \text{and} \quad \sum_{n=1}^\infty a_n^2 < \infty\)	e.g. \(a_n = c/n\).
Multi-Level Splitting	\(\ell = P(E_1) \prod_{k=2}^m P(E_k \mid E_{k-1})\)	Decomposes rare event \(A = E_m\) into nested checkpoints.
Adaptive MCMC proposal	\(Y \sim N(X_n, \sigma_n^2 \mathbf{I})\)	\(\log \sigma_{n+1} = \log \sigma_n + \gamma_n (\alpha_n - \alpha^)\), target \(\alpha^ = 0.234\).
Diminishing Adaptation	\(\lim_{n \to \infty} \gamma_n = 0\)	Ensures transition kernels asymptotically stabilize.