Rare event importance sampling and multilevel splitting

Let $(P_{\lambda },A_{\lambda }{)}_{\lambda \in \mathrm{\Lambda }}$ be a family of probability measures and measurable events with $p_{\lambda }=P_{\lambda }(A_{\lambda })>0$ and $p_{\lambda }\to 0$ along the stated rare-event limit. An estimator ${\widehat{p}}_{\lambda }$ is unbiased with finite variance when $\mathbb{E}[{\widehat{p}}_{\lambda }]=p_{\lambda }$ and $\mathrm{Var}({\widehat{p}}_{\lambda })<\mathrm{\infty }$. Its relative error is $\mathrm{RE}({\widehat{p}}_{\lambda })=\frac{\sqrt{\mathrm{Var}({\widehat{p}}_{\lambda })}}{p_{\lambda }}\mathrm{.}$ It has bounded relative error if these relative errors remain bounded for all sufficiently rare indices. For such an unbiased estimator family, it is logarithmically efficient if ${\widehat{p}}_{\lambda }\ge 0$, its second moment is finite and positive, and ${\lim }_{\lambda }\frac{\log \mathbb{E}[{\widehat{p}}_{\lambda }^2]}{\log (p_{\lambda }^2)}=1.$ Among nonnegative unbiased estimator families, logarithmic efficiency is a weaker asymptotic property than bounded relative error; the definitions impose no implication between an algorithm's efficiency and its proposal-fitting success.

For a specified cost model, let $W_{\lambda }$ be the expected work of one estimator replication. Record both its variance and its variance-work product $W_{\lambda }\mathrm{Var}({\widehat{p}}_{\lambda })$; for relative precision, record $W_{\lambda }\mathrm{Var}({\widehat{p}}_{\lambda })\mathrm{/}{p}_{\lambda }^2$. These quantities depend on the cost model. Comparing algorithms requires the same declared work convention.

For a single event $A$ under a probability measure $P$, put $p=P(A)>0$. If $P$ has density $f$ with respect to $\mu$, the zero-variance event proposal is the probability density $q^{\ast }(x)=\frac{{\mathbf{1}}_A(x)f(x)}{p}\mathrm{.}$ It is a formal benchmark because its normalizing constant is the unknown event probability $p$. Ordinary importance-sampling estimators and their weighted-second-moment condition are defined in Ordinary importance sampling and Ordinary importance sampling has finite variance and a CLT under a weighted second moment.

Let $\varphi (x)=(2\pi {)}^{-1\mathrm{/}2}{e}^{-x^2\mathrm{/}2}$ be the standard-normal density and $\bar{\mathrm{\Phi }}(x)={\int }_x^{\mathrm{\infty }}\varphi (t)dt$. For every $x>0$, $\frac{x}{x^2+1}\varphi (x)\le \bar{\mathrm{\Phi }}(x)\le \frac{\varphi (x)}{x}\mathrm{.}$ Consequently $\bar{\mathrm{\Phi }}(x)\sim \varphi (x)\mathrm{/}x$ and $\log \bar{\mathrm{\Phi }}(x)\sim -x^2\mathrm{/}2$ as $x\to \mathrm{\infty }$.

Let $(\mathsf{X},\mathcal{X},\mu )$ be a measure space, let $A\in \mathcal{X}$, and let the target probability measure $P$ have a pointwise finite measurable density $f:\mathsf{X}\to [0,\mathrm{\infty })$ with respect to $\mu$, so $P(B)={\int }_{B}fd\mu$. Let $\{q_{\theta }:\theta \in \mathrm{\Theta }\}$ be pointwise finite probability densities with respect to $\mu$, and write $Q_{\theta }(dx)=q_{\theta }(x)\mu (dx)$, where $\mathrm{\Theta }\subset {\mathbb{R}}^d$ is a finite nonempty set equipped with lexicographic order. For each $\theta \in \mathrm{\Theta }$, supply an exact admissibility checker $\mathcal{C}(A,f,q_{\theta })$ that returns admissible if and only if $P{\mathrm{\mid }}_A\ll Q_{\theta }$ and the contribution

is integrable under $Q_{\theta }$, where $P{\mathrm{\mid }}_A(B):=P(B\cap A)$ for measurable $B$. Supply an initial parameter ${\theta }_0\in \mathrm{\Theta }$, a pilot size $m\ge 1$, and a production size $n\ge 1$.

1. Draw an IID pilot $Y_1,\dots ,Y_m\sim q_{{\theta }_0}$.
2. If no $Y_j$ lies in $A$, return pilot-no-hit.
3. For every $\theta \in \mathrm{\Theta }$, compute $s_{\theta }(y)=\{\begin{array}{ll}\log q_{\theta }(y),& y\in A,\\ 0,& y\mathrm{\notin }A,\end{array}S(\theta )={\sum }_{j=1}^m{s}_{\theta }(Y_j),$ with $\log 0=-\mathrm{\infty }$. Pointwise finiteness rules out $+\mathrm{\infty }$, and the piecewise score avoids the undefined product $0\cdot (-\mathrm{\infty })$. If all scores equal $-\mathrm{\infty }$, return no-finite-maximizer. Otherwise choose the lexicographically least maximizer $\widehat{\theta }$ among the finite scores and freeze it.
4. Evaluate $\mathcal{C}(A,f,q_{\widehat{\theta }})$. If it returns inadmissible, return inadmissible-production-proposal. Otherwise draw an IID production sample $X_1,\dots ,X_n\sim q_{\widehat{\theta }}$, independent of the pilot conditional on $\widehat{\theta }$, and return the ordinary importance-sampling estimate ${\widehat{p}}_{\mathrm{C}\mathrm{E}}=\frac{1}{n}{\sum }_{i=1}^n{w}_{\widehat{\theta }}(X_i)$ using the piecewise-defined contribution $w_{\widehat{\theta }}$ above.

The finite candidate set makes the score comparisons and tie rule deterministic. The supplied exact checker is an input oracle; the procedure does not infer support or integrability from a pilot fit. If its check fails, no production estimate or unbiasedness or finite-variance claim is returned.

Let $(\mathrm{\Omega },\mathcal{F},P)$ be a probability space as in Probability measures and probability spaces. A finite nested event sequence of length $L\ge 1$ is a sequence of measurable sets $\mathrm{\Omega }=A_0\supseteq A_1\supseteq \cdots \supseteq A_L=A$ such that $P(A_{\mathrm{\ell }})>0$ for $0\le \mathrm{\ell }\le L$. Its conditional level probabilities and level laws are $p_{\mathrm{\ell }}=\frac{P(A_{\mathrm{\ell }})}{P(A_{\mathrm{\ell }-1})}\in (0,1],{\pi }_{\mathrm{\ell }}(B)=\frac{P(B\cap A_{\mathrm{\ell }})}{P(A_{\mathrm{\ell }})}(B\in \mathcal{F})\mathrm{.}$ The target probability is $p=P(A_L)$. Independence across levels means that the complete sample arrays used at distinct levels are mutually independent; correct marginal laws alone do not imply this condition.

Let $(E,\mathcal{E})$ be a standard Borel space, $M$ a probability kernel on it, and $T,R\in \mathcal{E}$ disjoint target and return sets. Write $S=E\setminus (T\cup R)$ and let $C=\{\mathtt{h}\mathtt{i}\mathtt{t},\mathtt{r}\mathtt{e}\mathtt{t}\mathtt{u}\mathtt{r}\mathtt{n},\mathrm{\partial }\}$ have the discrete sigma-algebra; $\mathrm{\partial }$ is the cemetery outcome meaning the budget expired before either set was reached.

For budget $m=0$, define ${\kappa }_0(x,\cdot )=\{\begin{array}{ll}{\delta }_{\mathtt{h}\mathtt{i}\mathtt{t}},& x\in T,\\ {\delta }_{\mathtt{r}\mathtt{e}\mathtt{t}\mathtt{u}\mathtt{r}\mathtt{n}},& x\in R,\\ {\delta }_{\mathrm{\partial }},& x\in S\mathrm{.}\end{array}$ For $m\ge 0$, define the finite-budget stopped-status kernel recursively by ${\kappa }_{m+1}(x,\cdot )=\{\begin{array}{ll}{\delta }_{\mathtt{h}\mathtt{i}\mathtt{t}},& x\in T,\\ {\delta }_{\mathtt{r}\mathtt{e}\mathtt{t}\mathtt{u}\mathtt{r}\mathtt{n}},& x\in R,\\ {\displaystyle {\int }_E{\kappa }_m(y,\cdot )M(x,dy),}& x\in S\mathrm{.}\end{array}$ Thus a target hit at the initial state counts at budget zero; each additional budget unit permits one more transition. The exact-hitting status kernel, when used, is the limit law of these statuses and treats nontermination as $\mathrm{\partial }$; it is an exact sampler for $P_x({\tau }_T<{\tau }_R)$ only when the requested event is the finite first hit represented by that limit. The budgeted kernel is distinct from an exact-hitting oracle.

Let $\nu$ be a probability measure on $(E,\mathcal{E})$, $M$ a probability kernel on $E$, $G,f:E\to [0,1]$ measurable, $z\in [0,1]$, and $N\ge 1$. Draw $Y_1,\dots ,Y_N$ IID from the probability law $\nu M$. Put $a=N^{-1}{\sum }_{i=1}^{N}G(Y_i)$. If $a>0$, set ${\eta }^N=\frac{{\sum }_{i=1}^{N}G(Y_i){\delta }_{Y_i}}{Na};$ if $a=0$, let ${\eta }^N$ be any probability measure. Define the one-step unnormalized estimate ${\mathrm{\Gamma }}^N(f)=za{\eta }^N(f)$. Then $\mathbb{E}[{\mathrm{\Gamma }}^N(f)]=z\nu M(Gf),\mathrm{Var}({\mathrm{\Gamma }}^N(f))\le \frac{z^2}{4N}\mathrm{.}$ The estimate is the normalizer update $za$ multiplied by the weighted empirical mean of $f$. The zero-weight fallback does not affect ${\mathrm{\Gamma }}^N$.

The cited rare-event SMC papers motivate the particle update; this exact one-step bound is proved locally from the canonical IID sample-mean moments, including the zero-weight branch.

Let $(\mathsf{X},\mathcal{X},\mu )$ be a measure space, let $P$ have density $f\ge 0$ with $\int fd\mu =1$, and let $A\in \mathcal{X}$ have $p=P(A)>0$. For $q^{\ast }(x)={\mathbf{1}}_A(x)f(x)\mathrm{/}p$ and $X\sim Q^{\ast }$ with this density, the ordinary importance-sampling contribution $Z={\mathbf{1}}_A(X)f(X)\mathrm{/}{q}^{\ast }(X)$, defined on the $Q^{\ast }$-support, equals $p$ almost surely. Thus the estimator from $n\ge 1$ IID $Q^{\ast }$-draws is unbiased and has variance zero. This proposal depends on the unknown $p$, so it is a benchmark and not a generally available algorithm.

Assume AC. Let $P=N(0,1)$, $A_a=[a,\mathrm{\infty })$, $p_a=P(A_a)=\bar{\mathrm{\Phi }}(a)$, and $Q_a=N(a,1)$, for $a>0$. Let $Z_a={\mathbf{1}}_{A_a}(X)dP\mathrm{/}d{Q}_a(X)$ for $X\sim Q_a$, and let ${\widehat{p}}_{a,n}$ be the average of $n$ IID copies. Then ${\mathbb{E}}_{Q_a}[Z_a]=p_a,{\mathbb{E}}_{Q_a}[Z_a^2]=e^{a^2}\bar{\mathrm{\Phi }}(2a)\mathrm{.}$ As $a\to \mathrm{\infty }$, ${\widehat{p}}_{a,n}$ is logarithmically efficient, while its one-sample relative error grows asymptotically as $(\pi \mathrm{/}2{)}^{1\mathrm{/}4}\sqrt{a}$. In particular this tilt is not bounded-relative-error.

Use Independent-pilot cross-entropy proposal fitting and production. Let $\mathcal{G}=\sigma (Y_1,\dots ,Y_m)$ and let $H\in \mathcal{G}$ be the event that the pilot succeeds and the algorithm's exact admissibility checker accepts the selected parameter $\widehat{\theta }$. On $H$, the checker guarantees $P{\mathrm{\mid }}_A\ll Q_{\widehat{\theta }}$ and integrability under $Q_{\widehat{\theta }}$ of the piecewise contribution $w_{\widehat{\theta }}$ defined in the algorithm. Define ${\tilde{p}}_{\mathrm{C}\mathrm{E}}$ to equal the returned production estimate on $H$ and zero off $H$. Then

In particular, the production estimate is conditionally unbiased almost surely on the successful, admissible pilot event. This gives no probability-of-success, proposal-optimization, finite-variance, logarithmic-efficiency, or total-work guarantee.

For the nested event levels of Nested rare-event levels and conditional laws, $P(A_L)={\prod }_{\mathrm{\ell }=1}^L{p}_{\mathrm{\ell }}\mathrm{.}$

Input the nested levels and exact conditional samplers for ${\pi }_{\mathrm{\ell }-1}$ from Nested rare-event levels and conditional laws, a particle count $N\ge 1$, and a declared per-draw work cost $c_{\mathrm{\ell }}>0$ for each level. If any conditional sampler is unavailable, return conditional-sampler-unavailable.

For each level $\mathrm{\ell }=1,\dots ,L$, independently of every other level, draw $X_{\mathrm{\ell },1},\dots ,X_{\mathrm{\ell },N}$ IID from ${\pi }_{\mathrm{\ell }-1}$. Set ${\widehat{p}}_{\mathrm{\ell }}=N^{-1}{\sum }_{i=1}^N{\mathbf{1}}_{A_{\mathrm{\ell }}}(X_{\mathrm{\ell },i}),\widehat{p}={\prod }_{\mathrm{\ell }=1}^L{\widehat{p}}_{\mathrm{\ell }}\mathrm{.}$ Return $\widehat{p}$, the level counts, total draw work $W=N{\sum }_{\mathrm{\ell }=1}^L{c}_{\mathrm{\ell }}$, and status success if all level counts are positive or zero-level-estimate otherwise. The zero estimate is not a proof that $P(A_L)=0$.

The procedure is an ideal independent-level baseline. It does not assert that a branching or Markov mutation method produces independent exact conditional draws.

Input the nested levels of Nested rare-event levels and conditional laws on a standard Borel state space $(\mathrm{\Omega },\mathcal{F})$, a particle count $N\ge 1$, and probability kernels $K_{\mathrm{\ell }}$ for $1\le \mathrm{\ell }<L$. Each mutation kernel must satisfy $K_{\mathrm{\ell }}(x,A_{\mathrm{\ell }})=1(x\in A_{\mathrm{\ell }}),{\pi }_{\mathrm{\ell }}{K}_{\mathrm{\ell }}={\pi }_{\mathrm{\ell }}\mathrm{.}$ If a required kernel is unavailable or violates its declared conditions, return mutation-kernel-unavailable.

Apply the canonical full-resampling sequential importance-resampling procedure of Sequential importance resampling for Feynman–Kac path targets to the following finite Feynman–Kac model: initial law $P$, initial potential $G_0\equiv 1$, identity transition $M_1$, event potential $G_1={\mathbf{1}}_{A_1}$, and, for $2\le \mathrm{\ell }\le L$, transition $M_{\mathrm{\ell }}=K_{\mathrm{\ell }-1}$ and potential $G_{\mathrm{\ell }}={\mathbf{1}}_{A_{\mathrm{\ell }}}$. Use multinomial full resampling at every selection and proposal kernels equal to the target transitions, so each Radon–Nikodym correction is one.

Return the SIR normalizer estimate $\widehat{p}$ and the full particle output. If any incremental average event weight is zero, report status zero-particle-prefix and $\widehat{p}=0$, as in the canonical SIR procedure. This status does not assert that the target probability is zero. The exact conditional level sampler is not inferred from the existence of an invariant kernel; kernel accessibility is an input requirement.

Assume AC. Under the hypotheses of Finite-budget stopped-event kernel with a cemetery outcome, every ${\kappa }_m$ is a measurable probability kernel from $E$ to $C$. For an initial probability law $\mu$, let $p_m={\int }_E{\kappa }_m(x,\{\mathtt{h}\mathtt{i}\mathtt{t}\})\mu (dx)\mathrm{.}$ Then $p_m=P_{\mu }({\tau }_T<{\tau }_R,{\tau }_T\le m)$. In the one-step Feynman–Kac model with initial law $\mu$, transition ${\kappa }_m$, initial potential $1$, and terminal potential ${\mathbf{1}}_{\{\mathtt{h}\mathtt{i}\mathtt{t}\}}$, the normalizer is exactly $p_m$. The limit ${\lim }_{m\to \mathrm{\infty }}{p}_m=P_{\mu }({\tau }_T<{\tau }_R)$ is the exact-hitting normalizer; an exact simulator with finite expected work requires a separate termination-cost bound.

Assume AC. In the stopped-event setup of Finite-budget stopped-event kernel with a cemetery outcome, with $(E,\mathcal{E})$ standard Borel and initial probability law $\mu$, let $p=P_{\mu }({\tau }_T<{\tau }_R)$ and let $p_m$ be the finite-budget normalizer of Finite-budget stopped-kernel measurability and event normalizer. Then $0\le p-p_m=P_{\mu }({\tau }_T<{\tau }_R,{\tau }_T>m)\le P_{\mu }({\tau }_T\wedge {\tau }_R>m)\mathrm{.}$ Thus the finite-budget stopped normalizer estimates the truncated event; a rate for its bias requires a tail bound for the stopping time. Separately, assume AC and let both finite-horizon Feynman–Kac models have the same finite horizon $H\in {\mathbb{N}}_0$, the same standard Borel state spaces $({\mathsf{X}}_j,{\mathcal{X}}_j)$ for $0\le j\le H$, the same initial law ${\mu }_0$, and the same potentials $G_j:{\mathsf{X}}_j\to [0,1]$ for $0\le j\le H$. Let their transition kernels be $M_j,{\tilde{M}}_j:{\mathsf{X}}_{j-1}⇝{\mathsf{X}}_j$ for $1\le j\le H$, and suppose

where $\mathrm{\parallel }\nu -{\nu }^{\mathrm{\prime }}{\mathrm{\parallel }}_{\mathrm{T}\mathrm{V}}={\sup }_B\mathrm{\mid }\nu (B)-{\nu }^{\mathrm{\prime }}(B)\mathrm{\mid }$. If $Z_H^M$ and $Z_H^{\tilde{M}}$ are the normalizers in the common Feynman–Kac setup of Feynman–Kac path targets and normalizers, then $\mathrm{\mid }{Z}_H^M-Z_H^{\tilde{M}}\mathrm{\mid }\le {\sum }_{j=1}^H{ϵ}_j$. If a grid-detected crossing event $A_h$ is contained in a continuous-time target event $A$, its probability bias is exactly $P(A\setminus A_h)$; no grid-error rate follows without a separate path regularity or crossing-probability bound.

Run Independent fixed-level splitting on the nested levels of Nested rare-event levels and conditional laws with $N\ge 1$ draws per level. Then $\mathbb{E}[\widehat{p}]=p,\frac{\mathrm{Var}(\widehat{p})}{p^2}={\prod }_{\mathrm{\ell }=1}^L(1+\frac{1-p_{\mathrm{\ell }}}{N{p}_{\mathrm{\ell }}})-1,W=N{\sum }_{\mathrm{\ell }=1}^L{c}_{\mathrm{\ell }}\mathrm{.}$ In absolute error, $\mathrm{Var}(\widehat{p})\le {(1+\frac{1}{4N})}^L-1.$ For fixed $L$, the right side is $O(L\mathrm{/}N)$, uniformly in the rarity of $p$; the relative-variance formula has factors $1\mathrm{/}{p}_{\mathrm{\ell }}$ and gives no uniform relative-efficiency claim.

Assume AC. Use the full-resampling SIR model in Full-resampling fixed-level SIR splitting with invariant mutations with $N\ge 1$, positive level probabilities, and each $K_{\mathrm{\ell }}$ supported on $A_{\mathrm{\ell }}$ and invariant for ${\pi }_{\mathrm{\ell }}=P(\cdot \mid A_{\mathrm{\ell }})$. Then every prefix Feynman–Kac normalizer is $Z_{\mathrm{\ell }}=P(A_{\mathrm{\ell }})$, and the product normalizer estimate $\widehat{p}$ is unbiased: $\mathbb{E}[\widehat{p}]=P(A_L)\mathrm{.}$ This conclusion asserts no finite-$N$ variance bound, accuracy guarantee, or uniform-in-rarity efficiency.

Assume AC. Use the full-resampling fixed-level SIR model of Full-resampling fixed-level SIR splitting with invariant mutations on nested levels $A_0\supseteq \cdots \supseteq A_L$, with $N\ge 1$ particles and positive finite prefixes $Z_p=P(A_p)>0$. Its event potentials satisfy $0\le G_p\le 1$. Then $\mathbb{E}[\widehat{p}]=P(A_L),\mathbb{E}[(\widehat{p}-P(A_L){)}^2]=\mathrm{Var}(\widehat{p})\le \frac{L+1}{4N}\mathrm{.}$ The bound is an absolute mean-square error bound for fixed $L$; it is not a relative-efficiency claim.