3.1. Variance of the Estimator

First, note that

is unbiased since

As done in [11], the variance of the estimator

is derived by induction and the variance for k thresholds is given by

where

represents the estimator of P in a simulation with k thresholds.

Clearly, the formula holds in straightforward simulation (i.e., when k = 0), since

is a renormalized sum of i.i.d. Bernoulli variables with parameter P.

To go from k to k + 1, assume (12); thus, we have to prove that this formula holds for k + 1 thresholds. First, note that for all X and Y random variables, which are independent given the set B and X σ(B)-measurable, we have

Now let

The random variables X_i0 are i.i.d. with common law Ber(P₁), and conditionally at the event B₁, X_i0 and Z_i0 are independent. Note that each Z_i0 is the estimator of P in a model with k thresholds, T₂ to T_k+1 for the trajectory issued from the success of X_i0. Thus,

and by the induction hypothesis,

So applying (13) with X ∼ Ber(P₁) and Z ∼ Z_i0, we have

Thus, for M thresholds,

Remark 3.1: The induction principle has a concrete interpretation: If in a simulation with M steps, the retrials generated in the first level are not taken into account except one that we call the main trial, we have a simulation with M − 1 steps.

3.2. Optimization of the Parameters

As stated in Section 1, our aim is to minimize the variance for a fixed budget, giving optimal values for N,R₁,…,R_M, P₁,…,P_M+1, and M. Therefore, we have to describe the cost of a given simulation: Each time a particle is launched, it generates an average cost function h. We assume the following:

• The cost h for a particle to reach B_i starting from B_i−1 depends only on P_i (not on the starting level).
• h is decreasing in x (which means that the smaller the transition probability is, the harder the transition is and the higher the cost is).
• h is nonnegative.
• h converges to a constant (in general, small) when x converges to 1.

The (average) cost is then

where N_i is the number of trials that have reached threshold i. Finally,

Example 3.1: We want to study the model of the simple random walk on

starting from zero that we kill as soon as it reaches the level −1 or k (success if we reach k, failure otherwise).

So let X_n be such that

, where {Y_n} is a sequence of random variables valued in {−1,1} with

and define T_k = inf{n ≥ 0 : X_n = −1 or k}.

One can easily check that X_n and X_n² − n are martingales. By Doob's stopping theorem,

, which yields

(i.e., the cost needed to reach the next level is (1/p) − 1 if p is the success probability).

To minimize the variance of

, the optimal values are derived in three steps:

1. The optimal values of N,R₁,…,R_M are derived when we consider that P₁,…,P_M+1 are constant (i.e., the thresholds B_i are fixed).
2. Replacing these optimal values in the variance, we derive the optimal transition probabilities: P₁,…,P_M+1.
3. Replacing these optimal values in the variance, we derive M, the optimal number of thresholds.

Optimal values for N, R₁,…, R_M. Using the method of Lagrange multipliers, we get

Optimal values for P₁,…,P_M+1. Thus, the variance becomes

Proceeding as previously under the constraint P = P₁···P_M+1, we obtain that all of the P_i's satisfy

. If we assume that there exists a unique solution to this equation, we have P_i = g(λ); hence, P = g(λ)^M+1 and g(λ) = P^1/(M+1). Finally,

Optimal value for M. The optimal values for P₁,…,P_M+1 imply that the optimal R_i becomes 1/P_i, i = 1,…,M; thus,

which we want to minimize in M. Note that R_i P_i = 1. Let

whose derivative cancels in

In general, this does not give an integer. We have y₀ = ln P/(M + 1) (i.e., M + 1 = [ln P/y₀] or [ln P/y₀] + 1). Let ln P/y₀ = n + x with 0 < x < 1. Then the following hold:

• If we take M + 1 = n, y = ln P/n.
• If we take M + 1 = n + 1, y = ln P/(n + 1).

The value of the ratio ρ := f (n − 1)/f (n) gives the best choice for M as follows:

• If ρ < 1, M = n − 1.
• If ρ > 1, M = n.

Thus, the optimal number of thresholds is given by M = [ln P/y₀] − 1 or M = [ln P/y₀], where y₀ solves F(y) = 0. Then M minimizes

Example 3.2: For h = 1, we have to solve y = 2(e^y − 1). We get y₁ = 0 and y₂ ≈ −1.5936. y₂ is a minimum and the optimal value of M is

With P = 10^−k, we have

Note that M increases while P decreases, and with this value of M, each R_i and P_i become

Thus, the optimal sampling number and the optimal transition probabilities are independent of the rare event probability.

Moreover, asymptotically, M = n = [ln P/y₀] − 1; thus,

Application 3.1: In approximate counting, remember that the goal is to estimate the number of Knapsack solutions (i.e., the cardinal of Ω defined by

for a given positive real vector a = (a_i)_i=1ⁿ and real number b). We might try to apply the Markov chain Monte Carlo method (MCMC) [9]: Construct a Markov chain

with state space Ω = {x ∈ {0,1}ⁿ : a·x ≤ b} and transitions from each state x = (x₁,…,x_n) ∈ Ω defined by the following:

• With probability ½, let y = x; otherwise
• select i uniformly at random in {1,…,n} and let y′ = (x₁,…,x_i−1, 1 − x_i,x_i+1,…,x_n)
• If ay′ ≤ b, then let y = y′, else let y = x.

The new state is y. This random walk on the hypercube truncated by the hyperplane a·x = b converges to the uniform distribution over Ω. This suggests a procedure for selecting Knapsack solutions almost uniformly at random. Starting in state (0,..,0), simulate

for sufficiently many steps that the distribution over states is “close”

The problem is to bound the number of steps necessary to make the Markov chain

“close” to stationarity. More precisely, we need a bound of the mixing time:

where Δ_x(t) = max_S⊂Ω|P^t(x,S) − Π(S)| and Π is the stationary distribution. In [7], it is shown that

steps suffice.

Keep the vector a fixed but allow b to vary. Use

instead of

to emphasize the dependence on b. Assume without loss of generality that a₁ ≤ ··· ≤ a_n and define

. One can check that

Now write

The ratio ρ_i = |Ω(b_i)|/|Ω(b_i+1)| may be estimated by sampling almost uniformly from Ω(b_i+1) using the Markov chain

and computing the fraction of the samples that lie within Ω(b_i).

Now take a = [1,2,3,4], b = 3, h = 1, R = 5, and C = 2600. We chose the levels as follows: First, define b₁ = 0, b₂ = 1, b₃ = 3, b₄ = 3, and b₅ = b; second, define B₀ = Ω, B₁ = Ω(b₄), B₂ = Ω(b₃), B₃ = Ω(b₂), and B₄ = Ω(b₁). Thus, here, M = n − 1, N = C/n, and n_step = 1020. Obviously, Card(Ω) = 5. We run three different simulations: The first, suggested in [7], consists of estimating the n ratios independently, the crude and splitting ones. We obtain different estimations for Card(Ω):

• Estimation by crude simulation = 4.088
• Estimation by the n ratios independently = 5.44
• Estimation by splitting = 5.019

Even though the levels are not optimal, splitting provides an improvement.

Let us describe briefly the possible solutions of (30). Remember that we want to solve (30); that is, if z = e^y and z ≠ 0,1,

First, let z₀ be the solution of 2(z − 1) = ln z. Since h′ ≤ 0, H is negative and a quick survey shows that L is positive on]0,z₀ [and negative on]z₀,1[. As a consequence, the solutions of (37) lie in]z₀,1[, if they exist. Thus, solving (30) is equivalent to studying the intersections between H and L. A quick survey of these functions shows that we have two cases (see Fig. 3):

Note that y = 0 is a solution of (30). In case 1, it corresponds to a maximum, and in case 2, it corresponds to a minimum. The second case is excluded since we made the assumption h(1) > 0.

Behavior of H and L.

Remark 3.2: The solution y = 0 corresponds to the following optimal values:

However P_i = 1 implies that P = 1 and R_i = 1 means that we just perform a crude simulation.

Example 3.3: Here, P = 10⁻¹² and C = 10⁴.

Thus, the control of the variance of

gives a crude confidence interval for P. Indeed, we get

This estimation is, in general, useless. For example, for h = 1, M = 12, and α = 10⁻², the upper bound becomes ≈ 5 × 10⁵/N. To obtain a bound lower than 1, we need N ≥ 5 × 10⁵. To improve it, we will use Chernoff's bounding method instead of the Markov inequality: For all λ > 0,

where

is the log-Laplace of

. Optimization on λ > 0 provides

Similarly,

Let ψ* be the Crämer transform of ψ: ψ*(τ) = sup_λ[λτ − ψ(λ)]. Thus,

So we want to obtain an accurate lower bound of ψ*.

Remark 3.3: Although we would therefore like to take R_i so that R_i P_i = 1, we are constrained to choose R_i to be a positive integer. Hereafter, we suppose that we are in the optimal case, where R_i = 1/P_i is an integer.

4.1. Description of the Model and First Results

We consider a Galton–Watson model (Z_n), where the size of the nth-generation Z_n is the number of particles that have reached the level B_n, with one particle run at the beginning. Then Z₀ = 1 and Z_n satisfies the following recurrence relation:

where X_iⁿ is the number of particles among R_i that have reached the (n + 1)-st level. The (X_iⁿ)_n≥1 are i.i.d. with common law Binomial, with parameters (R_n,P_n+1) and X_i⁰ ∼ Ber(P₁). Take the optimal values of Section 3.2:

Let

, the g.f. of Z₁. Then the g.f. of Z_n is the nth iterate of f. Since

, we get

where g is the g.f. of a Ber(P₀) and f the g.f. of a Bin(R,P₀). Thus, we are interested in the expression of f_M, the Mth functional iterate of f.

Here,

, so we are in the critical case of the branching process that ensures that the algorithm of the simulation stops with probability 1 when M → ∞ (see [6]), since if f⁽³⁾(1) < ∞,

This emphasizes the rarity character when the number M of thresholds increases and the probabilities between the levels decrease.

In our case,

The iterated function f_M has no explicit tractable form and we will derive bounds for f_M(s) around s = 1. To do this, we state a general result on the Laplace transform in critical Galton–Watson models, which we could not find in the literature.

4.2. Bounds of f_n(s) for 0 ≤ s < 1 and m = 1

Remark 4.1: Remember that f_n and its derivatives are convex. Furthermore, for all 0 ≤ s ≤ 1, s ≤ f (s) ≤ f (1) = 1, and by induction, f (s) ≤ f₂(s) ≤ ··· ≤ 1. Finally, we obtain f_n(s) → 1 since f_n(s) ≥ f_n(0).

Proposition 4.1: Let α₁ = f′′(1)/2, C = (max_s∈[0,1] f′′′(s))/6α₁, and γ_n = nα₁ [1 − (C/n)(log n + 1)] − α₁. Then, for s close to 1 and large n,

Proof: Upper bound: Using Taylor's expansion, with f_n(s) ≤ θ_n ≤ f_n(1) = 1,

since f′(1) = 1. Let r_n = 1 − f_n(s); r_n satisfies

Now let α₀ = f′′(0)/2. Define the decreasing sequences (a_n) and (b_n) satisfying

Then

By substituting (66) in

, we get

Lower bound: In fact, we prove by induction that

For n = 1, the left-hand side of (68) is given by Remark 4.1. Then note that

; thus, for n large enough,

. For all 1 − (1/n) ≤ s ≤ 1,

However, by induction, we have

, and so, since f is increasing, taking

,

(i.e.,

, where

. Note that

; more precisely, we have

and we finally obtain the left-hand side of (58) since γ → h_γ is increasing. █

In the particular case of f (s) = (P₀ s + 1 − P₀)^R, we can derive a more precise lower bound:

Proposition 4.2: For s close to 1,

Observe that this is precise at s = 1.

Proof: Let h(s) = 1 − ((1 − s)/(1 + α₁(1 − s))). Since f (1) = h(1) = 1, f′(1) = h′(1) = 1, and f′′(1) = h′′(1) = 2α₁, the sign of f − h trivially depends on the sign of the third derivative of f − h, which is obviously negative here. Then h ≤ f. Since f is increasing, we deduce (74) by induction. █

We plot in Figure 4a the upper bound and the two lower bounds for P = 10⁻¹² and s near 1.

Bounds of fM(s).

4.3. Bounds of f_n(s) for 1 ≤ s and m = 1

Remark 4.2: First, let us note that, by convexity, for all s ≥ 1,

hence, f (s) ≥ s, and by induction on n,

We remark that for s > 1, the iterated function increases rapidly to infinity.

Proposition 4.3: Let γ_n′ = nα₁ [1 + (C/n)(log n + 1)] − α₁. Then, for s close to 1 and large n,

Proof: Proceeding as in Proposition 4.1 leads to the upper bound. Here, f_n →_n→+∞ ∞, which prevents us from making a Taylor expansion around 1. To overcome this difficulty, consider k_n, the inverse function of f_n; it is the nth functional iterate of the g.f. k (inverse function of f) that takes the value 1 in 1, whose derivative is 1 in 1 and second derivative is negative, and k_n →_n→+∞ 1. Thus, making a Taylor development and using the same tools as previously, we get

where β₂ = k′′(s)/2 and s_n := 1 + (1/nα₁). Using the link between k_n and f_n and the upper bound of k_n,

The lower bound of k_n leads to an upper bound of f_n. However, it provides no improvement. █

As done earlier, we can derive a more precise upper bound in the particular case of f (s) = (P₀ s + 1 − P₀)^R:

Proposition 4.4: For s close to 1,

We plot in Figure 4b these three bounds for P = 10⁻¹² and s near 1.

About the geometric distribution. If the law of X is such that the probabilities p_k are in a geometric proportion (

for k = 1,2…. and p₀ = 1 − p₁ − p₂… with b,c > 0 and b ≤ 1 − c), then the associated g.f. is a rational function:

Taking b = (1 − c)² and c = α₁ /(1 + α₁) leads to

So we have compared the nth functional iterate of a Binomial g.f. to the one of a geometric g.f. It suggests comparing the importance splitting models with Binomial and with geometric laws. The geometric laws model is set in the following way: We run particles one after the other. As long as the next level is not reached, we keep on generating particles; then we start again from the level the particle is at (the geometric distribution is the law of the first success).

This link is also stressed by Cosnard and Demongeot in [4]: for m = 1 and σ² = f′′(1) = 2α₁, the asymptotic behavior of f²ⁿ is the same as the geometric distribution with the same variance (i.e., h).

4.4. Optimization of the Crämer Transform

Remember that

Considering the gradient of the functions, we prove that the supremum for λ ≥ 0 is reached near zero. So we can use the upper bounds for f_M obtained in the previous subsection, which leads to lower bounds for ψ*:

where

Finally,

One can easily obtain explicit but complex expressions for F(x) and G(x). We plot in Figure 5 the upper bounds obtained by the variance and by the Laplace transform, for different values of α, the prescribed error of the confidence interval. We take P = 10⁻⁹ and the optimal values obtained above for the parameters. Note that the upper bound given by the Laplace transform is better than the bound given by Chebychev's inequality, with the variance. We obtain

. In the preceding example where P = 10⁻⁹, if we fix α = 0.05 and L close to 0.01, then the corresponding costs needed are 3 × 10⁷ for the variance and 3 × 10⁶ for the Laplace transform.

Upper bounds obtained by the variance and the Laplace transform.