<figure><img src="image_1.png" /></figure>

HAL
open science

# Efficiency of Estimating Functions of Means In Rare-Event Contexts

Marvin K Nakayama, Bruno Tuffin

<figure><img src="image_2.png" /></figure>

## To cite this version:

Marvin K Nakayama, Bruno Tuffin. Efficiency of Estimating Functions of Means In Rare-Event Contexts. Winter Simulation Conference, Dec 2023, San Antonio, United States. hal-04096711

HAL Id: hal-04096711

https://inria.hal.science/hal-04096711

Submitted on 13 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficiency of Estimating Functions of Means In Rare-Event Contexts

Marvin K. Nakayama

Computer Science Department

New Jersey Institute of Technology

Newark, USA

Bruno Tuffin

Inria, Univ Rennes, CNRS, IRISA

Rennes, FRANCE

## Abstract

When estimating a function of means, where not necessarily all of them correspond to rare events, we provide conditions under which having efficient estimators of each individual mean leads to an efficient estimator of the function of the means. We illustrate this setting through several examples, and numerical results complement the theory.

## 1 Introduction

Many rare-event simulation problems consider an estimand expressed in terms of several quantities, where the behaviors of some but not necessarily all of their estimators are critically influenced by the rarity of the event of interest. For example, for a stable GI/GI/1 queue, the expected time for the queue-length process to first hit a large threshold can be written as a ratio. The denominator is the probability (i.e., expectation of an indicator) of hitting the large threshold within a busy period, which is a small rare-event probability that is difficult to estimate. But the numerator, which is the expectation of the minimum of the busy cycle length and the time to hit the large threshold, turns out to be an easily estimated non-rare mean. Thus, the numerator may be efficiently estimated via crude Monte Carlo, but the denominator requires applying a variance-reduction technique [1, Chapters V and VI] to efficiently estimate.

In this paper, we consider estimating a function of means, some (but not necessarily all) are intricately related to rare events. We provide conditions that ensure when a good estimator for each mean leads to an efficient estimator of the function of means.

The rest of the paper unfolds as follows. Section 2 reviews the concept of relative error and different notions of efficiency for simulation estimators for problems involving rare events. Section 3 describes the setting of the estimand as a function of means, and presents our main theorem. Throughout these sections, we provide examples to motivate the ideas and notation. Section 4 gives numerical results complementing the theory, and concluding remarks appear in Section 5.
## 2 Relative Error

Let $\alpha$ be a performance measure or estimand for a stochastic model. Suppose that we estimate $\alpha$ via an estimator $\hat{\alpha}_n$ based on a sample size of $n$ when using some Monte Carlo (MC) method, such as crude MC (simple random sampling) or a variance-reduction technique. Also, suppose that $\hat{\alpha}_n$ obeys a central limit theorem (CLT)

$$ \sqrt{n}[\hat{\alpha}_n - \alpha] \Rightarrow N(0, \tau^2) \quad \text{as } n \to \infty, \qquad (1) $$

where $\Rightarrow$ denotes convergence in distribution [2, Section 25], $N(a, b^2)$ represents a normal random variable with mean $a$ and variance $b^2$, and $\tau^2$ is the asymptotic variance of the CLT (1). If we have a consistent estimator $\hat{\tau}_n^2$ of $\tau^2$ in the sense that $\hat{\tau}_n^2 \Rightarrow \tau^2$ as $n \to \infty$, then for fixed large $n$, we can construct an approximate 95% confidence interval (CI) for $\alpha$ based on (1) as $[\hat{\alpha}_n \pm 1.96\hat{\tau}_n/\sqrt{n}]$, whose relative half-width is roughly $1.96\tau/[\alpha\sqrt{n}]$. This motivates assessing the quality of the estimator $\hat{\alpha}_n$ through its relative error (RE; [7]), defined for a fixed $n$ and $\alpha \neq 0$ as

$$ \mathrm{RE}[\hat{\alpha}_n] = \frac{\tau}{|\alpha|}. \qquad (2) $$

Thus, obtaining a 95% confidence interval with relative half-width $w$ (say, $w = 0.1$) requires a sample size $n$ that is roughly $n \approx (1.96\tau)^2/(w\alpha)^2$, which is $(1.96/w)^2$ times the squared RE.

**Example 1** (Estimating an expectation via a sample mean) When $\alpha = E[Z]$ is a mean of a random variable $Z$ that can be simulated, we can draw $n$ independent and identically distributed (i.i.d.) observations $Z_1, Z_2, \dots, Z_n$ of $Z$, and compute the estimator as the sample average $\hat{\alpha}_n = (1/n) \sum_{k=1}^n Z_k$ of $\alpha$. Then the CLT (1) holds for $\tau^2 = \mathrm{Var}[Z] \in (0, \infty)$ [2, Theorem 27.1]. This set-up also allows for applying importance sampling (IS; Section V.1 of [1]), as follows. Suppose that $\alpha = E_0[Q]$ for some random variable $Q$, where $E_0[\cdot]$ is the expectation operator under the original distribution $F_0$ of $Q$. Now let $F_1$ be another distribution such that the probability measure corresponding to $F_0$ is absolutely continuous with respect to the measure for $F_1$. Then for $E_1[\cdot]$ as the expectation operator under $F_1$, we may write $\alpha = \int Q dF_0 = \int Q \frac{dF_0}{dF_1} dF_1 = E_1[QL]$ with $L = dF_0/dF_1$ as the likelihood ratio, so $Z = QL$ and $E[\cdot]$ corresponds to $E_1[\cdot]$.

When $\alpha$ is related to a rare event, the relative error of an estimator can become quite large as the event of interest becomes rarer. For example, suppose that we want to estimate the probability that a random variable $W_{\langle r \rangle}$ lies in some set $A_{\langle r \rangle} \subset \Re$, where the distribution of $W_{\langle r \rangle}$ and the set $A_{\langle r \rangle}$ may depend on some “rarity parameter” $r$. In the context of Example 1, we take $Z \equiv Z_{\langle r \rangle} = I(W_{\langle r \rangle} \in A_{\langle r \rangle})$, where $I(\cdot)$ denotes the indicator function. Thus, we consider a sequence of models indexed by $r$, and examine what happens as $r \to \infty$, where we assume the event $\{W_{\langle r \rangle} \in A_{\langle r \rangle}\}$ becomes “rarer” as $r$ grows in the sense that $\alpha_{\langle r \rangle} \equiv E[Z_{\langle r \rangle}] = P(W_{\langle r \rangle} \in A_{\langle r \rangle}) \to 0$ as $r \to \infty$. For example, in a stable GI/GI/1 queue, $r$ might represent a large buffer threshold, and $W_{\langle r \rangle}$ is the maximum queue length before first emptying out; letting $A_{\langle r \rangle} = \{r, r+1, \dots\}$, we have that $\alpha_{\langle r \rangle} = P(W_{\langle r \rangle} \in A_{\langle r \rangle})$ is the probability that the queue length hits level $r$ during a busy cycle, where $\alpha_{\langle r \rangle} \to 0$ as $r \to \infty$ [14, Theorem 1]. For each fixed $r$, the CLT in (1) holds for Example 1 with $Z_{\langle r \rangle} = I(W_{\langle r \rangle} \in A_{\langle r \rangle})$ and asymptotic variance $\tau^2 = \tau_{\langle r \rangle}^2 = \alpha_{\langle r \rangle}(1 - \alpha_{\langle r \rangle})$. Let $\hat{\alpha}_n \equiv \hat{\alpha}_{\langle r \rangle,n}$ be the resulting average of a fixed
number $n$ of i.i.d. copies of $Z_{\langle r\rangle}$, and its relative error in (2) then behaves as

$$ \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\tau_{\langle r \rangle}}{\alpha_{\langle r \rangle}} = \frac{\sqrt{\alpha_{\langle r \rangle}(1 - \alpha_{\langle r \rangle})}}{\alpha_{\langle r \rangle}} = \frac{\sqrt{1 - \alpha_{\langle r \rangle}}}{\sqrt{\alpha_{\langle r \rangle}}} \to \infty \quad \text{as } r \to \infty. $$

Thus, the relative error of $\hat{\alpha}_{\langle r \rangle, n}$ blows up as the event of interest becomes rarer, which means that the required sample size to obtain a fixed relative-width confidence interval grows as roughly $1/\alpha_{\langle r \rangle}$ as $r \to \infty$. This illustrates the primary difficulty with rare-event simulation.

Previous works (e.g., [13] and Chapter VI of [1]) for specific stochastic models devise variance-reduction techniques that produce estimators of such rare-event probabilities $\alpha_{\langle r \rangle}$ for which the relative error grows very slowly, remains bounded or vanishes as $r \to \infty$. Specifically, in Example 1, suppose that for each fixed $r$, we have an estimator $\hat{\alpha}_{\langle r \rangle, n}$ of $\alpha_{\langle r \rangle}$ satisfying the CLT (1) with asymptotic variance $\tau_{\langle r \rangle}^2 = O(\alpha_{\langle r \rangle}^2)$ as $r \to \infty$, where the notation $b_1(r) = O(b_2(r))$ as $r \to \infty$ for functions $b_1$ and $b_2$ means that there exists constants $c_0 > 0$ and $r_0 > 0$ such that $|b_1(r)| \le c_0|b_2(r)|$ for all $r > r_0$. Then the relative error of the estimator $\hat{\alpha}_{\langle r \rangle, n}$ satisfies

$$ \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\sqrt{O(\alpha_{\langle r \rangle}^2)}}{\alpha_{\langle r \rangle}} = \frac{O(\alpha_{\langle r \rangle})}{\alpha_{\langle r \rangle}} = O(1) \quad \text{as } r \to \infty, \quad (3) $$

so the estimator has bounded relative error (BRE). Hence, the necessary sample size to obtain a confidence interval of relative half-width $w$ remains bounded as $r \to \infty$. If instead $\tau_{\langle r \rangle}^2 = o(\alpha_{\langle r \rangle}^2)$ as $r \to \infty$, where the notation $b_1(r) = o(b_2(r))$ means that $\lim_{r \to \infty} b_1(r)/b_2(r) = 0$, then we have

$$ \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\sqrt{o(\alpha_{\langle r \rangle}^2)}}{\alpha_{\langle r \rangle}} = \frac{o(\alpha_{\langle r \rangle})}{\alpha_{\langle r \rangle}} = o(1) \quad \text{as } r \to \infty, \quad (4) $$

so the estimator then has vanishing relative error (VRE).

We next consider an efficiency property that is slightly weaker than BRE. To do this, define

$$ \mathrm{RE}_{\zeta}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\tau_{\langle r \rangle}}{|\alpha_{\langle r \rangle}|^{\zeta}} \quad (5) $$

for each constant $\zeta > 0$, so (5) with $\zeta = 1$ reduces to just RE in (2). Then when $\alpha_{\langle r \rangle} \to 0$ as $r \to \infty$, we say that the estimator $\hat{\alpha}_{\langle r \rangle, n}$ of $\alpha_{\langle r \rangle}$ is logarithmically efficient (LE) if

$$ \limsup_{r \to \infty} \mathrm{RE}_{1-\delta}[\hat{\alpha}_{\langle r \rangle, n}] = 0 \quad (6) $$

for all $\delta \in (0, 1)$. Often a desirable goal when designing an estimator of a small probability based on large-deviations theory, the property (6) is equivalent to $\liminf_{r \to \infty} \frac{|\ln \tau_{\langle r \rangle}|}{|\ln(\alpha_{\langle r \rangle})|} \ge 1$, which shows where LE gets its name. While the literature on rare-event simulation has primarily focused on LE estimators for small probabilities, some estimands of interest have $\alpha_{\langle r \rangle} \to \infty$ as $r \to \infty$, e.g., a
stochastic process's expected hitting time to a rarely visited set of states. When instead considering an estimand for which $\alpha_{\langle r \rangle} \to \pm\infty$ as $r \to \infty$, we redefine LE to mean that for all $\delta > 0$,

$$ \limsup_{r \to \infty} \mathrm{RE}_{1+\delta}[\hat{\alpha}_{\langle r \rangle, n}] = 0, \qquad (7) $$

which is equivalent to $\limsup_{r \to \infty} \frac{|\ln \tau_{\langle r \rangle}|}{|\ln(\alpha_{\langle r \rangle})|} \le 1$. We do not define LE for the case when $\lim_{r \to \infty} \alpha_{\langle r \rangle} = \chi_0 \in (-\infty, \infty)$ with $\chi_0 \ne 0$. Indeed, the idea of LE originates from studies of rare-event probabilities through large-deviations theory, where a goal in simulation is to construct an estimator whose standard deviation and mean both converge to 0 at the same exponential rate as $r \to \infty$. While (7) extends this conceptual framework to handle $\alpha_{\langle r \rangle} \to \pm\infty$, an appropriate definition for LE is not obvious when $\alpha_{\langle r \rangle}$ converges to a finite $\chi_0 \ne 0$, so we do not consider this case. We could also take into account the CPU time to compute $\hat{\alpha}_{\langle r \rangle, n}$ via the work-normalized RE (WNRE) [7], but we focus here only on RE.

In what follows, we will sometimes omit the subscript $\langle r \rangle$ on variables to simplify the notation.

## 3 Functions of Means

We will consider estimands $\alpha$ having a specific form:

$$ \alpha = g(\boldsymbol{\theta}) \qquad (8) $$

for some known function $g: \Re^d \to \Re$ and $\boldsymbol{\theta} = (\theta_1, \theta_2, \dots, \theta_d) \in \Re^d$ is a vector of unknown parameters, for some $d \ge 1$. Both $g$ and $d$ do not depend on the rarity parameter $r$, but the components $\theta_i$ of $\boldsymbol{\theta}$ may. We assume that there is a $d$-dimensional random vector $\mathbf{X} = (X_1, X_2, \dots, X_d)$ such that its mean

$$ E[\mathbf{X}] = \boldsymbol{\theta} \qquad (9) $$

and that we can generate i.i.d. copies of $\mathbf{X}$.

The random vector $\mathbf{X}$ has some joint distribution. The $d$ components of $\mathbf{X}$ may be dependent, as is often the case when they are generated from the same simulation, or they may be independent, generated from different simulations applying different MC methods. For example, $X_1$ could be generated using crude MC, and $X_2 = QL$ could be generated (independently of $X_1$) using importance sampling, as in Example 1. This approach is called *measure-specific importance sampling* (MSIS) [5].

Many estimands of interest have the form in (8).

**Example 2 (Expected hitting time)** Consider a nondelayed regenerative process $Y = (Y(t): t \ge 0)$ evolving on a state space $S \subset \Re^v$ for some $v \ge 1$ [6]. We split the sample path of $Y$ into *regenerative cycles*, where the cycles are i.i.d., and let $T_0 > 0$ be the time of the first regeneration. For example, if $Y$ represents the queue-length process in a stable GI/GI/1 queue, then regenerations occur when a customer arrives to an empty system. If $Y$ is an irreducible and positive-recurrent continuous-time Markov chain (CTMC), returns to any fixed state $s_0 \in S$ constitute regenerations. Also consider a
set $A \subset S$ (with $s_0 \notin A$ in the case of a CTMC), and let $T_A = \inf\{t \ge 0: Y(t) \in A\}$ be the hitting time to A. We are interested in estimating $\alpha = E[T_A]$. As shown in e.g., [5] and [4], we can express $\alpha$ as a ratio $\alpha = \eta/\gamma$, with $\eta = E[\min(T_A, T_0)]$ and $\gamma = E[I(T_A < T_0)] = P(T_A < T_0)$, where both $\min(T_A, T_0)$ and $I(T_A < T_0)$ are measurable with respect to the sigma-field of the process Y up to time $T_0$. Then for the random vector $\mathbf{X} = (X_1, X_2)$ in (9) has dimension $d = 2$, with $X_1 = \min(T_A, T_0)$ and $X_2 = I(T_A < T_0)$, (8) has $\boldsymbol{\theta} = (\eta, \gamma)$, and $g(x_1, x_2) = x_1/x_2$.

**Example 3 (Derivative of expected hitting time)** Suppose for the CTMC Y in Example 2, its infinitesimal generator matrix is parameterized by $\lambda$, which may represent the transition rate of a specific group of transitions. For example, $\lambda$ may represent the failure rate of a component in a reliability setting, or the arrival rate of customers to a queue. In this case, we write $\beta = E[T_A]$ as $\beta(\lambda) = \eta(\lambda)/\gamma(\lambda)$ as these values now depend on $\lambda$. We are interested in computing the derivative $\beta'(\lambda) = [\eta'(\lambda)\gamma(\lambda) - \eta(\lambda)\gamma'(\lambda)]/\gamma^2(\lambda)$, where $\eta'(\lambda) = \frac{d}{d\lambda}\eta(\lambda)$ and $\gamma'(\lambda) = \frac{d}{d\lambda}\gamma(\lambda)$. Using the likelihood ratio derivative method [3], we can express $\eta'(\lambda) = E[\min(T_A, T_0)L']$ and $\eta'(\lambda) = E[I(T_A < T_0)L']$, where $L'$ is the derivative with respect to $\lambda$ of the likelihood ratio $L$ up to $\min(T_A, T_0)$; see [10] for details. Thus, for random vector $\mathbf{X} = (X_1, X_2, X_3, X_4)$ in (9) with dimension $d = 4$, and $X_1 = \min(T_A, T_0)$, $X_2 = I(T_A < T_0)$, $X_3 = \min(T_A, T_0)L'$, and $X_4 = I(T_A < T_0)L'$, we can write $\beta'(\lambda)$ as $\alpha$ in (8) using $\boldsymbol{\theta} = (\eta(\lambda), \gamma(\lambda), \eta'(\lambda), \gamma'(\lambda))$, and $g(x_1, x_2, x_3, x_4) = [x_3x_2 - x_4x_1]/x_2^2$.

**Example 4 (Conditional expectation)** Let $(Z_1, Z_2)$ be a random vector, and let $A \subset \Re$ be a set. Then the conditional expectation $\alpha = E[Z_1 | Z_2 \in A] = E[Z_1I(Z_2 \in A)]/E[I(Z_2 \in A)]$ may be expressed in the form in (8) by considering the random vector $\mathbf{X} = (X_1, X_2)$ in (9) with dimension $d = 2$, $X_1 = Z_1I(Z_2 \in A)$, $X_2 = I(Z_2 \in A)$, $\boldsymbol{\theta} = (E[X_1], E[X_2])$, and $g(x_1, x_2) = x_1/x_2$.

As in Section 2, we consider a sequence of problems indexed with by rarity parameter $r$, and consider the limit as $r \to \infty$. In (8), we assume that the value of the unknown parameter $\boldsymbol{\theta} = \boldsymbol{\theta}_{\langle r \rangle} = (\theta_{\langle r \rangle,1}, \theta_{\langle r \rangle,2}, \dots, \theta_{\langle r \rangle,d})$ may depend on $r$, but the function $g$ does not. Then the estimand $\alpha = \alpha_r$ in (8) also depends on $r$, and we want to study when $\alpha_{\langle r \rangle}$ can be estimated with BRE as in (3), with VRE as in (4), or with LE as in (6) or (7). As noted before in Section 2, we will sometimes omit the subscript $\langle r \rangle$ to simplify notation.

To estimate $\alpha_{\langle r \rangle} = g(\boldsymbol{\theta}_{\langle r \rangle})$ in (8), we first generate $n \ge 1$ i.i.d. observations $\mathbf{X}_{\langle r \rangle, k}$, $k = 1, 2, \dots, n$, of $\mathbf{X} = \mathbf{X}_{\langle r \rangle} = (X_{\langle r \rangle,1}, X_{\langle r \rangle,2}, \dots, X_{\langle r \rangle,d})$, whose joint distribution may depend on $r$. We then obtain an unbiased estimator $\hat{\boldsymbol{\theta}}_{\langle r \rangle, n} = (\hat{\theta}_{\langle r \rangle, n,1}, \hat{\theta}_{\langle r \rangle, n,2}, \dots, \hat{\theta}_{\langle r \rangle, n,d}) = (1/n) \sum_{k=1}^n \mathbf{X}_{\langle r \rangle, k}$ of $\boldsymbol{\theta} = \boldsymbol{\theta}_{\langle r \rangle} = E[\mathbf{X}_{\langle r \rangle}]$, leading to the plug-in estimator

$$ \hat{\alpha}_{\langle r \rangle, n} = g(\hat{\boldsymbol{\theta}}_{\langle r \rangle, n}) \qquad (10) $$

of $\alpha_{\langle r \rangle}$. This framework allows for the components of $\mathbf{X}_{\langle r \rangle}$ to be sampled using different simulation methods, as in MSIS; e.g., $X_{\langle r \rangle,1}$ may be generated using crude MC, while $X_{\langle r \rangle,2} = Q_{\langle r \rangle}L_{\langle r \rangle}$ is sampled via importance sampling, where $L_{\langle r \rangle}$ is the likelihood ratio, as in Example 2 with dimension $d = 2$, $Q_{\langle r \rangle} = I(T_{A_{\langle r \rangle}} < T_{0,\langle r \rangle})$, and $\theta_{\langle r \rangle,2} = E_0[Q_{\langle r \rangle}] = E_1[Q_{\langle r \rangle}L_{\langle r \rangle}]$.

For the estimator $\hat{\alpha}_{\langle r \rangle, n}$, we will study its relative error as in (2). We want to determine conditions under which the RE properties of an estimator $\hat{\boldsymbol{\theta}}_{\langle r \rangle, n}$ of $\boldsymbol{\theta}_{\langle r \rangle}$ will carry over to the plug-in estimator $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\boldsymbol{\theta}}_{\langle r \rangle, n})$ of $\alpha_{\langle r \rangle}$. Let $\Sigma_{\langle r \rangle} = (\Sigma_{\langle r \rangle, i,j: i, j = 1, 2, \dots, d})$ be the covariance matrix of the random
vector $X_{\langle r \rangle}$, where $\Sigma_{\langle r \rangle, i, j} = \text{Cov}[X_{\langle r \rangle, i}, X_{\langle r \rangle, j}] = E[(X_{\langle r \rangle, i} - E[X_{\langle r \rangle, i}]) (X_{\langle r \rangle, j} - E[X_{\langle r \rangle, j}])]$. Also, let $\nabla g(\boldsymbol{\theta}) = (g_1(\boldsymbol{\theta}), g_2(\boldsymbol{\theta}), \dots, g_d(\boldsymbol{\theta}))$ be the gradient (when it exists) of $g(\boldsymbol{\theta})$, where $g_i(\boldsymbol{\theta}) = \frac{\partial}{\partial \theta_i} g(\boldsymbol{\theta})$, $i = 1, 2, \dots, d$, and for each $\boldsymbol{\psi}, \boldsymbol{\zeta} > 0$, define

$$h_{i, \psi, \zeta}(\boldsymbol{\theta}_{\langle r \rangle}) = \left| \frac{\theta_{\langle r \rangle, i}^{\psi} g_i(\boldsymbol{\theta}_{\langle r \rangle})}{g^{\zeta}(\boldsymbol{\theta}_{\langle r \rangle})} \right|. \quad (11)$$

For analyzing LE, we further impose the following.

**Assumption A1** As $r \to \infty$, $\alpha_{\langle r \rangle} \to \chi_0 \in \{0, \infty, -\infty\}$ and $\theta_{\langle r \rangle, i} \to \chi_i \in [-\infty, \infty]$ for each $i = 1, 2, \dots, d$.

Assumption A1 disallows the estimand $\alpha_{\langle r \rangle}$ to converge to a nonzero finite constant. This is to retain the spirit of the framework that LE was originally defined to consider, as discussed after (7). But for $i = 1, 2, \dots, d$, we permit some of the $\theta_{\langle r \rangle, i}$ to converge to nonzero finite constants $\chi_i$. This will allow us to handle a situation that occurs in Example 2 when estimating the expected time for the queue-length process of a stable GI/GI/1 queue to first hit a high level $r$; this will be explained later in Example 8. For each $\delta \in (0, 1)$ and $i = 0, 1, \dots, d$, define $\delta_i$ based on Assumption A1 as follows:

$$\delta_i = \begin{cases} -\delta & \text{if } \chi_i = 0, \\ 0 & \text{if } \chi_i \in (-\infty, \infty) \text{ with } \chi_i \neq 0, \\ \delta & \text{if } |\chi_i| = \infty. \end{cases} \quad (12)$$

**Theorem 1** Assume that for each $r$, the function $g$ in (8) is differentiable at $\boldsymbol{\theta}_{\langle r \rangle}$ with $\nabla g(\boldsymbol{\theta}_{\langle r \rangle}) \neq \mathbf{0}$ and $\Sigma_{\langle r \rangle}$ is finite and positive definite. Then the following hold as $r \to \infty$:

(i) If $h_{i,1,1}(\boldsymbol{\theta}_{\langle r \rangle})\text{RE}[\hat{\theta}_{\langle r \rangle,n,i}] = O(1)$ as $r \to \infty$ for all $i = 1, 2, \dots, d$, then $\hat{\alpha}_{\langle r \rangle,n}$ has BRE.

(ii) If $h_{i,1,1}(\boldsymbol{\theta}_{\langle r \rangle})\text{RE}[\hat{\theta}_{\langle r \rangle,n,i}] = o(1)$ as $r \to \infty$ for all $i = 1, 2, \dots, d$, then $\hat{\alpha}_{\langle r \rangle,n}$ has VRE.

For analyzing LE, if Assumption A1 also is in force, then the following holds as $r \to \infty$:

(iii) For each $i = 1, 2, \dots, d$, if

$$\limsup_{r \to \infty} h_{i,1+\delta_i,1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) \text{RE}_{1+\delta_i}[\hat{\theta}_{\langle r \rangle,n,i}] = 0 \quad (13)$$

for all $\delta \in (0, 1)$, then $\hat{\alpha}_{\langle r \rangle,n}$ is LE.

In Theorem 1(i) and (ii), if $h_{i,1,1}(\alpha_{\langle r \rangle})$ is bounded as $r \to \infty$ for all $i = 1, 2, \dots, d$, then the overall estimator $\hat{\alpha}_{\langle r \rangle,n}$ has BRE or VRE when each $\hat{\theta}_{\langle r \rangle,n,i}$ has that same property. To understand Theorem 1(iii), consider the special case when $h_{i,1+\delta_i,1+\delta_0}(\alpha_{\langle r \rangle})$ is bounded as $r \to \infty$ for all $i = 1, 2, \dots, d$. Note that if $\chi_i \in \{0, -\infty, \infty\}$, then $\limsup_{r \to \infty} \text{RE}_{1+\delta_i}[\hat{\theta}_{\langle r \rangle,n,i}] = 0$ means that $\hat{\theta}_{\langle r \rangle,n,i}$ is LE, as in (6) or (7); if instead $\chi_i \in (-\infty, \infty)$ with $\chi_i \neq 0$, then because $\delta_i = 0$ in (13), $\limsup_{r \to \infty} \text{RE}_{1+\delta_i}[\hat{\theta}_{\langle r \rangle,n,i}] = \limsup_{r \to \infty} \text{RE}[\hat{\theta}_{\langle r \rangle,n,i}] = 0$ means that $\hat{\theta}_{\langle r \rangle,n,i}$ has VRE. Thus, if $h_{i,1+\delta_i,1+\delta_0}(\alpha_{\langle r \rangle})$ is bounded as $r \to \infty$ for all $i = 1, 2, \dots, d$, and Assumption A1 holds, then the overall estimator $\hat{\alpha}_{\langle r \rangle,n}$ is LE when each $\hat{\theta}_{\langle r \rangle,n,i}$ is LE. But when $\chi_i \in (-\infty, \infty)$ with $\chi_i \neq 0$, Theorem 1 also allows for $\hat{\theta}_{\langle r \rangle,n,i}$ to be just BRE and not VRE if $h_{i,1+\delta_i,1+\delta_0}(\alpha_{\langle r \rangle}) \to 0$ as $r \to \infty$, which we will later see occurs in Example 8.
*Proof of Theorem 1.* For each fixed $r$, the condition that the covariance matrix $\Sigma_{\langle r \rangle}$ is finite and positive definite ensures that $\hat{\theta}_{\langle r \rangle, n}$ obeys a multivariate CLT [2, Theorem 29.5]

$$ \sqrt{n} [\hat{\theta}_{\langle r \rangle, n} - \theta_{\langle r \rangle}] \Rightarrow \Sigma_{\langle r \rangle}^{1/2} \mathcal{N}_d \quad \text{as } n \to \infty, \qquad (14) $$

where $\Sigma_{\langle r \rangle}^{1/2}$ is a $d \times d$ matrix satisfying $\Sigma_{\langle r \rangle}^{1/2} (\Sigma_{\langle r \rangle}^{1/2})^{\top} = \Sigma_{\langle r \rangle}$, with the superscript $\top$ denoting trans- pose, e.g., $\Sigma_{\langle r \rangle}^{1/2}$ is a Cholesky factor, and $\mathcal{N}_d$ is a $d$-dimensional normal random vector with mean vector $\mathbf{0} = (0, 0, \dots, 0)$ and whose covariance matrix is the $d \times d$ identity matrix. Under our assump- tion that $\nabla g(\theta_{\langle r \rangle}) \neq \mathbf{0}$, the delta method [15, p. 124] then yields the CLT

$$ \sqrt{n} [\hat{\alpha}_{\langle r \rangle, n} - \alpha_{\langle r \rangle}] \Rightarrow \tau_{\langle r \rangle} \mathcal{N}_1 \quad \text{as } n \to \infty, $$

with asymptotic variance

$$ \tau_{\langle r \rangle}^2 = \sum_{i=1}^{d} \sum_{j=1}^{d} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}. \qquad (15) $$

The Cauchy-Schwarz inequality implies that $|\Sigma_{\langle r \rangle, i, j}| \le \sqrt{\mathrm{Var}[X_{\langle r \rangle, i}] \mathrm{Var}[X_{\langle r \rangle, j}]}$, and also we have that $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n}] = \sqrt{\mathrm{Var}[X_{\langle r \rangle, i}]/|\theta_{\langle r \rangle, i}|}$, as in (2). Thus, the RE of $\hat{\alpha}_{\langle r \rangle, n}$ then satisfies

$$ \begin{align} \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] &= \left( \sum_{i=1}^{d} \sum_{j=1}^{d} \frac{g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \Sigma_{\langle r \rangle, i, j} \right)^{1/2} \le \left( \sum_{i=1}^{d} \sum_{j=1}^{d} \left| \frac{g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \Sigma_{\langle r \rangle, i, j} \right| \right)^{1/2} \nonumber \\ &\le \left( \sum_{i=1}^{d} \sum_{j=1}^{d} \left| \frac{g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \right| \sqrt{\mathrm{Var}[X_{\langle r \rangle, i}] \mathrm{Var}[X_{\langle r \rangle, j}]} \right)^{1/2} \tag{16} \\ &= \left( \sum_{i=1}^{d} \sum_{j=1}^{d} \left| \frac{\theta_{\langle r \rangle, i} g_i(\theta_{\langle r \rangle}) \theta_{\langle r \rangle, j} g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \right| \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}] \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, j}] \right)^{1/2} \nonumber \\ &= \left( \sum_{i=1}^{d} \sum_{j=1}^{d} h_{i, 1, 1}(\theta_{\langle r \rangle}) \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}] h_{j, 1, 1}(\theta_{\langle r \rangle}) \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, j}] \right)^{1/2} \nonumber \end{align} $$

by (11), so parts (i) and (ii) follow from our assumptions about each $h_{i, 1, 1}(\theta_{\langle r \rangle}) \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}]$ being $O(1)$ or $o(1)$.

For part (iii) under Assumption A1, we simplify the discussion by considering the square of (5)
for $\zeta = 1 + \delta_0$. As in (16), we get

$$ \frac{\tau_{\langle r \rangle}^2}{|\alpha_{\langle r \rangle}|^{2+2\delta_0}} \le \sum_{i=1}^d \sum_{j=1}^d \left| \frac{g_i(\boldsymbol{\theta}_{\langle r \rangle}) g_j(\boldsymbol{\theta}_{\langle r \rangle})}{g^{2+2\delta_0}(\boldsymbol{\theta}_{\langle r \rangle})} \right| \sqrt{\mathrm{Var}[X_{\langle r \rangle, i}] \mathrm{Var}[X_{\langle r \rangle, j}]} \\ = \sum_{i=1}^d \sum_{j=1}^d \left| \frac{\boldsymbol{\theta}_{\langle r \rangle, i}^{1+\delta_i} g_i(\boldsymbol{\theta}_{\langle r \rangle}) \boldsymbol{\theta}_{\langle r \rangle, j}^{1+\delta_j} g_j(\boldsymbol{\theta}_{\langle r \rangle})}{g^{1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) g^{1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle})} \right| \mathrm{RE}_{1+\delta_i}[\hat{\boldsymbol{\theta}}_{\langle r \rangle, n, i}] \mathrm{RE}_{1+\delta_j}[\hat{\boldsymbol{\theta}}_{\langle r \rangle, n, j}] \\ = \sum_{i=1}^d \sum_{j=1}^d h_{i,1+\delta_i,1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) \mathrm{RE}_{1+\delta_i}[\hat{\boldsymbol{\theta}}_{\langle r \rangle, n, i}] h_{j,1+\delta_j,1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) \mathrm{RE}_{1+\delta_j}[\hat{\boldsymbol{\theta}}_{\langle r \rangle, n, j}] $$

by (11), so part (iii) holds by (6) and (7) under our condition (13). $\square$

Even though $\hat{\boldsymbol{\theta}}_{\langle r \rangle, n}$ is an unbiased estimator of $\boldsymbol{\theta}_{\langle r \rangle}$, we may have that $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\boldsymbol{\theta}}_{\langle r \rangle, n})$ is a biased estimator of $\alpha_{\langle r \rangle} = g(\boldsymbol{\theta}_{\langle r \rangle})$ when $g$ is nonlinear, as in Examples 2–4, as well as the ones below. However, for fixed $r$, we still obtain the CLT in (14) as $n \to \infty$ with the true estimand $\alpha_{\langle r \rangle}$ as the centering constant. Thus, for fixed $r$, the bias makes a negligible contribution to the mean-square error of $\hat{\alpha}_{\langle r \rangle, n}$ as $n \to \infty$. Indeed, under appropriate technical conditions on $g$ and the moments of the estimator of $\alpha_{\langle r \rangle}$,

$$ E[\hat{\alpha}_{\langle r \rangle, n}] = E[g(\hat{\boldsymbol{\theta}}_{\langle r \rangle, n})] = g(\boldsymbol{\theta}_{\langle r \rangle}) + \frac{1}{n} \sum_{i=1}^{d} \sum_{j=1}^{d} g_i(\boldsymbol{\theta}_{\langle r \rangle}) g_j(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} + o\left(\frac{1}{n}\right), $$

so the squared bias shrinks as $O(1/n^2)$, faster than the estimator's variance, which typically decreases as $c_0/n$ for some constant $c_0$.

## 3.1 Examples

We next provide examples (some expanding on our previous ones) that fit into the framework of Theorem 1.

**Example 5 (Mean time to failure)** We now specialize Example 2 for a highly reliable Markovian system (HRMS), as in [5] and [16]. An HRMS comprises a fixed collection of dependable components subject to random failures and repairs, where the component failure times and repair times are exponentially distributed. The stochastic process $Y$ models the evolution of the HRMS over time, and let $S$ be its state space, where each state $s \in S$ keeps track of the components that are currently up and down, along with any other information (e.g., queueing of failed components at the repair stations) necessary to ensure that $Y$ is a CTMC. The system begins at time 0 in $s_0 \in S$, which is the state with all components operational. The system fails when certain combinations of components are failed, and let $A \subset S$ denote the set of failed states, where we assume that $s_0 \notin A$, and the time to failure is $T_A$. To incorporate the rarity parameter $r$, we assume that the component failure rates are expressed as positive powers of $\epsilon = 1/r$, and repair rates are constants, and we analyze the HRMS as $r \to \infty$. Hence, as $r$ increases, individual components become more reliable, but the set $A$ of failed states does not change. The mean
time to failure (MTTF) $\alpha_{\langle r \rangle} = E_0[T_A] = E_0[T_{A, \langle r \rangle}]$ can be expressed as a ratio $\eta_{\langle r \rangle} / \gamma_{\langle r \rangle}$, where $\eta_{\langle r \rangle} = E_0[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})]$ and $\gamma_{\langle r \rangle} = E_0[I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})] = E_1[I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L]$, where $E_0$ is the expectation operator under the original system dynamics, and $E_1[\cdot]$ denotes expectation under importance sampling, with $L$ as the resulting likelihood ratio. We consider applying MSIS to estimate the ratio, with crude MC for estimating the numerator and importance sampling for handling the denominator. Therefore, in (9), the random vector $\mathbf{X}$ has dimension $d = 2$, and $\mathbf{X} = (X_1, X_2) = (\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle}), I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L)$, where the expectation operator $E[\cdot]$ uses $E_0[\cdot]$ (resp., $E_1[\cdot]$) on the first (resp., second) coordinate of $\mathbf{X}$. Also, (8)–(10) have $\theta_{\langle r \rangle} = (\eta_{\langle r \rangle}, \gamma_{\langle r \rangle})$, and $g(x_1, x_1) = x_1/x_2$. Since (11) has $h_{1,1,1}(\theta_{\langle r \rangle}) = h_{2,1,1}(\theta_{\langle r \rangle}) = 1$, BRE (resp., VRE) for the MTTF ratio estimator follows if the estimators for $\eta_{\langle r \rangle}$ and $\gamma_{\langle r \rangle}$ both have BRE (resp., VRE) by Theorem 1(i) and (ii). When applying balanced failure biasing (described further in Example 10) for the IS for estimating the denominator $\gamma_{\langle r \rangle}$, [16] establishes the BRE of the MTTF MSIS ratio estimator; [9], [13] and [8] provide related results for other forms of IS to estimate $\gamma_{\langle r \rangle}$.

**Example 6** For Example 3 specialized to an HRMS, [10] establishes the BRE of the derivative estimator of the MTTF with respect to a component failure rate when using MSIS with balanced failure biasing under certain assumptions.

**Example 7** (Alternative MTTF estimator) Instead of using the ratio estimator of the MTTF in Example 5, [11] propose a different approach, which can be more efficient in certain situations and also applies more generally, as for Example 2. Specifically, they express the MTTF as

$$ \alpha_{\langle r \rangle} = \frac{E_0[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} > T_{0, \langle r \rangle})] + E_0[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})]}{E_0[I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})]} \\ = \frac{E_0[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} > T_{0, \langle r \rangle})] + E_1[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L]}{E_1[I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L]} $$

Thus, to express this in the framework of (8)–(10), we apply MSIS by defining the random vector $\mathbf{X}_{\langle r \rangle} = (X_{\langle r \rangle,1}, X_{\langle r \rangle,2}, X_{\langle r \rangle,3})$ with dimension $d = 3$ and $X_{\langle r \rangle,1} = \min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} > T_{0, \langle r \rangle})$, $X_{\langle r \rangle,2} = \min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L$, and $X_{\langle r \rangle,3} = I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L$, where $X_{\langle r \rangle,1}$ is generated with crude MC, and $(X_{\langle r \rangle,2}, X_{\langle r \rangle,3})$ is generated with importance sampling. Hence, the expectation operator $E[\cdot]$ in (9) applies $E_0[\cdot]$ (resp., $E_1[\cdot]$) for the first (resp., second and third) coordinate(s), so $\theta = (\theta_1, \theta_2, \theta_3)$ with $\theta_1 = E_0[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} > T_{0, \langle r \rangle})]$, $\theta_2 = E_1[\min(T_{A, \langle r \rangle}, T_{0, \langle r \rangle})I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L]$, and $\theta_3 = E_1[I(T_{A, \langle r \rangle} < T_{0, \langle r \rangle})L]$, and $g(x_1, x_2, x_3) = (x_1 + x_2)/x_3$ in (8). Since (11) has $h_{1,1,1}(\theta_{\langle r \rangle}) = \frac{\theta_{\langle r \rangle,1}}{\theta_{\langle r \rangle,1} + \theta_{\langle r \rangle,2}} \le 1$, $h_{2,1,1}(\theta_{\langle r \rangle}) = \frac{\theta_{\langle r \rangle,2}}{\theta_{\langle r \rangle,1} + \theta_{\langle r \rangle,2}} \le 1$ and $h_{3,1,1}(\theta_{\langle r \rangle}) = 1$, they are all bounded, so by Theorem 1(i), the estimator of $\alpha_{\langle r \rangle}$ has BRE if the estimators of each of the three components have the same property. We can also establish if the estimator of $\alpha_{\langle r \rangle}$ is LE (which is more appropriate in, e.g., the GI/GI/1 setting) through an analysis that will be used in the next example.

**Example 8** (Expected time for GI/GI/1 queue length to reach a large threshold) For the queue-length process of a stable GI/GI/1 queue, we consider the estimand $\alpha_{\langle r \rangle} = E[T_{A, \langle r \rangle}]$, which is the expected hitting time to the set $A_{\langle r \rangle} = \{r, r+1, r+2, \dots\}$, as discussed in Example 2. We can then express $\alpha_{\langle r \rangle} = \eta_{\langle r \rangle} / \gamma_{\langle r \rangle}$, where $\eta_{\langle r \rangle} = E_0[\min(T_{A, \langle r \rangle}, T_0)]$ and $\gamma_{\langle r \rangle} = E_0[I(T_{A, \langle r \rangle} < T_0)] = E_1[I(T_{A, \langle r \rangle} < T_0)L]$, where $E_0[\cdot]$ is the expectation operator under the original system dynamics, and $E_1[\cdot]$ denotes expectation under importance sampling with $L$ as the resulting likelihood ratio. The distribution
of $T_0$ does not depend on $r$, but the distribution of $T_{A_{\langle r \rangle}}$ does. In (9), we apply MSIS using a random vector $\mathbf{X}$ of dimension $d=2$, with $\mathbf{X} = (\min(T_{A_{\langle r \rangle}}, T_0), I(T_{A_{\langle r \rangle}} < T_0)L)$, and the expectation $E[\cdot]$ applies $E_0[\cdot]$ (resp., $E_1$) to the first (resp., second) coordinate. Thus, $\boldsymbol{\theta}_{\langle r \rangle} = (\boldsymbol{\theta}_{\langle r \rangle,1}, \boldsymbol{\theta}_{\langle r \rangle,2}) = (\eta_{\langle r \rangle}, \gamma_{\langle r \rangle})$ and $g(x_1,x_1) = x_1/x_2$. In this case, we have that $0 < \eta_{\langle r \rangle} \le E[T_0] \in (0, \infty)$ for all $r$, and $\eta_{\langle r \rangle} \to \chi_1 \in (0, \infty)$ as $r \to \infty$, so $\delta_1 = 0$ in (12). Also, Theorem 1 of [14] shows that $\gamma_{\langle r \rangle} \to 0$ as $r \to \infty$, so $\delta_2 = -\delta$ in (12) for each $\delta > 0$. Hence, the estimand satisfies $\alpha_{\langle r \rangle} = \eta_{\langle r \rangle}/\gamma_{\langle r \rangle} \to \infty$ as $r \to \infty$, so $\delta_0 = \delta$ in (12). We now want to verify (13) for $i=1,2$. Note that $g_1(\boldsymbol{\theta}_{\langle r \rangle}) = 1/\theta_{\langle r \rangle,2} = 1/\gamma_{\langle r \rangle}$ and $g_2(\boldsymbol{\theta}_{\langle r \rangle}) = -\theta_{\langle r \rangle,1}/\theta_{\langle r \rangle,2} = -\eta_{\langle r \rangle}/\gamma_{\langle r \rangle}^2$. For $i=1$, because $\delta_1 = 0$ and $\delta_0 = \delta$, (11) becomes

$$h_{1,1+\delta_1,1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) = \left| \frac{\theta_{\langle r \rangle,1}^{1+\delta_1} g_1(\boldsymbol{\theta}_{\langle r \rangle})}{g^{1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle})} \right| = \left| \frac{\eta_{\langle r \rangle}/\gamma_{\langle r \rangle}}{(\eta_{\langle r \rangle}/\gamma_{\langle r \rangle})^{1+\delta}} \right| = \frac{1}{\alpha_{\langle r \rangle}^{\delta}} \to 0 \quad \text{as } r \to \infty \quad (17)$$

for each $\delta > 0$. We can also estimate $\theta_{\langle r \rangle,1} = \eta_{\langle r \rangle}$ with BRE using crude MC, so $\text{RE}_{1+\delta_1}[\hat{\theta}_{\langle r \rangle,n,1}] = \text{RE}[\hat{\theta}_{\langle r \rangle,n,1}] = O(1)$ as $r \to \infty$, which combined with (17) verifies (13) for $i=1$. For $i=2$, because $\delta_2 = -\delta$ and $\delta_0 = \delta$, (11) becomes

$$h_{2,1+\delta_2,1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) = \left| \frac{\theta_{\langle r \rangle,2}^{1+\delta_2} g_2(\boldsymbol{\theta}_{\langle r \rangle})}{g^{1+\delta_0}(\boldsymbol{\theta}_{\langle r \rangle})} \right| = \left| \frac{\gamma_{\langle r \rangle}^{1-\delta} (-\eta_{\langle r \rangle}/\gamma_{\langle r \rangle}^2)}{(\eta_{\langle r \rangle}/\gamma_{\langle r \rangle})^{1+\delta}} \right| = \frac{1}{\eta_{\langle r \rangle}^{\delta}}, \quad (18)$$

which, for each $\delta > 0$, is bounded for all $r$ sufficiently large. [14] provides an importance-sampling estimator for $\theta_{\langle r \rangle,2} = \gamma_{\langle r \rangle}$ that is LE, so this estimator satisfies $\limsup_{r \to \infty} \text{RE}_{1+\delta_2}[\hat{\theta}_{\langle r \rangle,n,2}] = \limsup_{r \to \infty} \text{RE}_{1-\delta}[\hat{\theta}_{\langle r \rangle,n,2}] = 0$ as in (6), which combined with (18) confirms (13) for $i=2$. Hence, the MSIS estimator of $\alpha_{\langle r \rangle}$ is LE by Theorem 1(iii).

## 3.2 Special Case when $d=1$

We now examine Theorem 1 in more detail for the special case when the dimension $d=1$, so in (8), the unknown parameter is $\boldsymbol{\theta}_{\langle r \rangle} = \theta_{\langle r \rangle,1} \equiv \theta_{\langle r \rangle}$, the estimand is $\alpha_{\langle r \rangle} = g(\theta_{\langle r \rangle})$, and the random vector $\mathbf{X}_{\langle r \rangle} = X_{\langle r \rangle,1} \equiv X_{\langle r \rangle}$ has mean $\theta_{\langle r \rangle}$. Theorem 1 then requires that the derivative $g'$ of $g$ is nonzero at $\theta_{\langle r \rangle}$, and the asymptotic variance in (15) becomes $\tau_{\langle r \rangle}^2 = [g'(\theta_{\langle r \rangle})]^2 \text{Var}[X_{\langle r \rangle}]$.

We first focus on BRE/VRE properties, which correspond to parts (i) and (ii) of Theorem 1. For $\hat{\theta}_{\langle r \rangle,n}$ as the estimator of $\theta_{\langle r \rangle}$, the relative error of the estimator $g(\hat{\theta}_{\langle r \rangle,n})$ of $g(\theta_{\langle r \rangle})$ is

$$\text{RE}[g(\hat{\theta}_{\langle r \rangle,n})] = \frac{|g'(\theta_{\langle r \rangle})| \sqrt{\text{Var}[X_{\langle r \rangle}]} |g(\theta_{\langle r \rangle})|}{|g(\theta_{\langle r \rangle})|} = \left| \frac{\theta_{\langle r \rangle} g'(\theta_{\langle r \rangle})}{g(\theta_{\langle r \rangle})} \right| \frac{\sqrt{\text{Var}[X_{\langle r \rangle}]} |g(\theta_{\langle r \rangle})|}{|\theta_{\langle r \rangle}|} = h(\theta_{\langle r \rangle}) \text{RE}[\hat{\theta}_{\langle r \rangle,n}], \quad (19)$$

where $h(\theta_{\langle r \rangle}) = |\theta_{\langle r \rangle} g'(\theta_{\langle r \rangle})/g(\theta_{\langle r \rangle})|$ and $\text{RE}[\hat{\theta}_{\langle r \rangle,n}] = \sqrt{\text{Var}[X_{\langle r \rangle}]}/|\theta_{\langle r \rangle}|$. Thus, $g(\hat{\theta}_{\langle r \rangle,n})$ has BRE (resp., VRE) if and only if $h(\theta_{\langle r \rangle}) \text{RE}[\hat{\theta}_{\langle r \rangle,n}]$ remains bounded (resp., vanishes) as $\varepsilon \to 0$. As discussed just before the proof of Theorem 1, a sufficient condition for $g(\hat{\theta}_{\langle r \rangle,n})$ to have BRE (resp., VRE) is that $\hat{\theta}_{\langle r \rangle,n}$ has BRE (resp., VRE) and $h(\theta_{\langle r \rangle})$ remains bounded as $\varepsilon \to 0$. But $g(\hat{\theta}_{\langle r \rangle,n})$ can
still have BRE or VRE even when $\hat{\theta}_{\langle r \rangle, n}$ does not if $h(\theta_{\langle r \rangle}) \to 0$ fast enough as $r \to \infty$. We next consider some specific examples.

First suppose that $g(\theta_{\langle r \rangle}) = c[\theta_{\langle r \rangle}]^v$ for constants $c \neq 0$ and $v \neq 0$. In this case, we get $h(\theta_{\langle r \rangle}) = |\theta_{\langle r \rangle}cv[\theta_{\langle r \rangle}]^{v-1}/(c[\theta_{\langle r \rangle}]^v)| = |v|$, so $h(\theta_{\langle r \rangle})$ is a constant and is trivially bounded, no matter how $\theta_{\langle r \rangle}$ depends on $r$. Thus, it then follows that $g(\hat{\theta}_{\langle r \rangle, n})$ has BRE or VRE if and only if $\hat{\theta}_{\langle r \rangle, n}$ does.

Now instead suppose that $g(\theta_{\langle r \rangle}) = ce^{\theta_{\langle r \rangle}}$ for some constant $c \neq 0$, so $h(\theta_{\langle r \rangle}) = |\theta_{\langle r \rangle}ce^{\theta_{\langle r \rangle}}/(ce^{\theta_{\langle r \rangle}})| = |\theta_{\langle r \rangle}|$. If $\theta_{\langle r \rangle} = c_1r^\ell$ for constants $c_1 \neq 0$ and $\ell > 0$, then $h(\theta_{\langle r \rangle}) = |c_1r^\ell| \to \infty$ as $r \to \infty$. Thus, it is possible to have unbounded $h(\theta_{\langle r \rangle})$. On the other hand, if $\theta_{\langle r \rangle} = c_1r^{-\ell}$ for $\ell > 0$, then $h(\theta_{\langle r \rangle}) = |c_1r^{-\ell}| \to 0$ as $r \to \infty$. In this case, if $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n}] = O(r^\ell)$ (resp., $o(r^\ell)$), then $g(\hat{\theta}_{\langle r \rangle, n})$ has BRE (resp., VRE), even though $\hat{\theta}_{\langle r \rangle, n}$ may have unbounded relative error.

For LE, focus now on part (iii) of Theorem 1. As in (19), with $\delta_0$ and $\delta_1$ defined in (12), we get

$$ \begin{aligned} \mathrm{RE}_{1+\delta_0}[g(\hat{\theta}_{\langle r \rangle, n})] &= \frac{|g'(\theta_{\langle r \rangle})|\sqrt{\mathrm{Var}[X_{\langle r \rangle}]} |g(\theta_{\langle r \rangle})|^{1+\delta_0}}{\delta_0} \\ &= \left| \frac{\theta_{\langle r \rangle}^{1+\delta_1} g'(\theta_{\langle r \rangle})}{(g(\theta_{\langle r \rangle}))^{1+\delta_0}} \right| \frac{\sqrt{\mathrm{Var}[X_{\langle r \rangle}]} |\theta_{\langle r \rangle}|^{1+\delta_1}}{\delta_0} \\ &= h_{1+\delta_1, 1+\delta_0}(\theta_{\langle r \rangle})\mathrm{RE}_{1+\delta_1}[\hat{\theta}_{\langle r \rangle, n}], \end{aligned} $$

where $h_{1+\delta_1, 1+\delta_0}(\theta_{\langle r \rangle}) = |\theta_{\langle r \rangle}^{1+\delta_1} g'(\theta_{\langle r \rangle})/(g(\theta_{\langle r \rangle}))^{1+\delta_0}|$. Here too, a sufficient condition for $g(\hat{\theta}_{\langle r \rangle, n})$ to have LE is that $\theta_{\langle r \rangle}$ has LE and $h_{1+\delta_1, 1+\delta_0}(\theta_{\langle r \rangle})$ remains bounded. Also if $h_{1+\delta_1, 1+\delta_0}(\theta_{\langle r \rangle}) \to 0$ sufficiently fast, it may be possible for $g(\hat{\theta}_{\langle r \rangle, n})$ to have BRE or VRE when $\hat{\theta}_{\langle r \rangle, n}$ is only LE. For example, again considering $g(\theta_{\langle r \rangle}) = ce^{\theta_{\langle r \rangle}}$ for some constant $c \neq 0$ and assuming $\theta_{\langle r \rangle} \to 0$ as $r \to \infty$ (so that $\delta_1 = -\delta$), we get $h_{1-\delta, 1+\delta_0}(\theta_{\langle r \rangle}) = |\theta_{\langle r \rangle}|^{1-\delta}/|g(\theta_{\langle r \rangle})|^{1+\delta_0}$. With $\delta_0 = 0$ and recalling that (5) with $\zeta = 1$ reduces to (2), we obtain $\mathrm{RE}[g(\hat{\theta}_{\langle r \rangle, n})] = |\theta_{\langle r \rangle}|^{1-\delta}\mathrm{RE}_{1-\delta}[\hat{\theta}_{\langle r \rangle, n}]/|g(\theta_{\langle r \rangle})| \to 0$ as $r \to \infty$, so $g(\hat{\theta}_{\langle r \rangle, n})$ has VRE.

# 4 Numerical Results

**Example 9** [12] describe a reliability system with two components, which we now consider. The first component fails after an exponential amount of time, and then with probability $1-p$, it automatically resets, and the process repeats again independently of the past. If the first component does not reset, which occurs with probability $p$, the second component takes over, and has a uniformly distributed lifetime. After the second component fails, the system fails, and then the entire system resets after an exponential amount of time. We model the system as a 3-state semi-Markov process (SMP) with state space $S = \{0, 1, 2\}$, whose embedded discrete-time Markov chain (DTMC) has transition probabilities as in Figure 1. The figure also shows the holding time distribution for each successive visit to a state $s \in S$, where $\mathcal{E}(\lambda)$ denotes an exponential distribution with rate $\lambda > 0$, so its mean is $1/\lambda$, and $\mathcal{W}(a, b)$ is a uniform distribution on the interval $(a, b)$. We define a rarity
parameter $\epsilon = 1/r$, and let $w_0, w_1$, and $\lambda_2$ be positive constants, whose values we will vary to study different behaviors of the model as $r \to \infty$.

We aim at computing the mean hitting time of state 2, which $\alpha_{\langle r \rangle}$ now represents. Typically, as explained in [12], if $w_0 + 1 > w_1$, the time spent in state 1 will make a negligible contribution to the hitting time compared to the total time spent in state 0 before moving to state 1, and both the estimator based on a ratio as in Example 5 or the alternative estimator of Example 7 will be efficient. In contrast, if $w_0 + 1 \le w_1$, the estimator in Example 5 encounters issues, motivating the approach in Example 7.

<figure><img src="image_4.png" /><figcaption>Figure 1: The edge labels are the transition probabilities of the embedded DTMC Z of an SMP X, with the holding-time distribution to each visit to a state below that state, where $w_0, w_1$, and $\lambda_2$ are constants.</figcaption></figure>

Adapting the notation in Example 7 to this specific example, we get that the rare set is $A = \{2\}$, $T_{0,\langle r \rangle}$ is the time to first return to state 0, and $T_{A,\langle r \rangle}$ is the system failure time, that is, the hitting time of state 2. Our experiments used $w_0 = 1$ and $w_1 = 3$, and we considered several values of $\epsilon = 1/r$, simulating with MSIS a total number $n = 10^6$ of cycles. For MSIS, we generate a number $n_1 = \lfloor \kappa n \rfloor$ (with $0 < \kappa < 1$) of cycles to estimate $\theta_{\langle r \rangle,1} = E_0[\min(T_{2,\langle r \rangle}, T_{0,\langle r \rangle})I(T_{2,\langle r \rangle} > T_{0,\langle r \rangle})]$ by crude Monte Carlo, while the remaining $n_2 = n - \lfloor \kappa n \rfloor$ cycles are used to estimate $\theta_{\langle r \rangle,2} = E_1[\min(T_{2,\langle r \rangle}, T_{0,\langle r \rangle})I(T_{2,\langle r \rangle} < T_{0,\langle r \rangle})L]$ and $\theta_{\langle r \rangle,3} = E_1[I(T_{2,\langle r \rangle} < T_{0,\langle r \rangle})L]$ by importance sampling, which replaces the transition probability $p = \epsilon = 1/r$ from state 0 to state 1 by a probability 0.8. The proportion $\kappa$ is computed from a presimulation to minimize the work-normalized variance of the resulting estimator (its expression is denoted $\gamma^*$ in Equation (23) of [11]), but keeping at least 10% of the cycles to both crude and IS simulations. This potentially unequal sampling allocation using $\kappa$ differs from the setting in (10), which assumes that each $\theta_i$ is estimated with the same sample size (but possibly different simulation methods), but this does not affect the asymptotic behavior (ignoring leading constants) since we constrain $\kappa \in (0.1, 0.9)$.

Table 1 displays the output for the estimation of $\hat{\theta}_{\langle r \rangle, n_i, i}$ of $\theta_{\langle r \rangle, i}$ ($i \in \{1, 2, 3\}$), with $n_3 = n_2$), an estimate $\widehat{RE}^{(i)}$ of $RE[\hat{\theta}_{\langle r \rangle, n_i, i}]$, and the resulting estimation $\hat{\alpha}_{\langle r \rangle, n}$ of the expected failure time with the corresponding estimated relative error and confidence interval. The results appear to show that the RE of the estimator of each $\theta_i$ does not grow as $\epsilon$ shrinks, and the same is true when estimating $\alpha_{\langle r \rangle} = E[T_{A,\langle r \rangle}]$. This behavior indicates BRE, which agrees with the theory in Theorem 1 and Example 7.

**Example 10** We consider an HRMS (as in Example 5) with $c = 3$ types of component and a redundancy 5 for each type. Each component has an exponentially distributed time to failure with rate
<table><caption>Table 1: Results for the MTTF in the example of Figure 1, from $n = 10^6$ independent cycles.</caption><thead><tr><th>ε = 1/r</th><th>θ̂<sub>(r),n1,1</sub></th><th>R̂E<sup>(1)</sup></th><th>θ̂<sub>(r),n2,2</sub></th><th>R̂E<sup>(2)</sup></th><th>θ̂<sub>(r),n3,3</sub></th><th>R̂E<sup>(3)</sup></th><th>α̂<sub>(r),n</sub></th><th>R̂E[α̂<sub>(r),n</sub>]</th><th>95% conf. int.</th></tr></thead><tbody><tr><td>1.00e-01</td><td>9.03e+00</td><td>1.10e+00</td><td>5.11e+01</td><td>8.06e-01</td><td>1.00e-01</td><td>5.00e-01</td><td>6.01e+02</td><td>7.09e-01</td><td>(5.999e+02, 6.016e+02)</td></tr><tr><td>1.00e-02</td><td>9.92e+01</td><td>1.01e+00</td><td>5.01e+03</td><td>8.15e-01</td><td>1.00e-02</td><td>4.99e-01</td><td>5.10e+05</td><td>6.69e-01</td><td>(5.098e+05, 5.111e+05)</td></tr><tr><td>1.00e-03</td><td>1.00e+03</td><td>1.00e+00</td><td>5.01e+05</td><td>8.16e-01</td><td>1.00e-03</td><td>4.99e-01</td><td>5.01e+08</td><td>6.78e-01</td><td>(5.007e+08, 5.021e+08)</td></tr><tr><td>1.00e-04</td><td>1.00e+04</td><td>1.00e+00</td><td>5.01e+07</td><td>8.16e-01</td><td>1.00e-04</td><td>4.99e-01</td><td>5.01e+11</td><td>6.80e-01</td><td>(4.998e+11, 5.012e+11)</td></tr></tbody></table>

$ε = 1/r$. Any failed component has an exponentially distributed repair time with rate 1, where there are 15 repairpersons, so no failed component ever has to queue before receiving service. Component failure and repair times are all independent. The system is down whenever fewer than two components of any one type are operational. We want to compute here the steady-state unavailability, so $α_{(r)}$ now represents the long run fraction of time that the resulting CTMC is down. As in [5], we can express $α_{(r)} = g(θ_{(r),1}, θ_{(r),2}) = θ_{(r),1}/θ_{(r),2}$, where the numerator is the expected down time in a cycle, and the denominator is the expected cycle length. Applying conditional Monte Carlo [1, Section V.4] by conditioning on the embedded DTMC, we can write $θ_{(r),1}$ as the mean of the sum of the expected sojourn times of each failed state visited within a cycle, and $θ_{(r),2}$ is the mean of the sum of sojourn times of all states visited within a cycle. We consider MSIS to estimate $θ_{(r),1}$ and $θ_{(r),2}$ from a total number $n = 10^6$ of independent cycles of the embedded DTMC, with $n_1 = ⌊κn⌋$ (resp. $n_2 = n - n_1$) cycles used to estimate $θ_{(r),1}$ (resp. $θ_{(r),2}$), with $κ$ obtained (as in Example 9) from a presimulation to minimize the work-normalized variance of the resulting estimator. We estimate $θ_{(r),1}$ through an approach known as *dynamic IS* [5], which within a cycle applies IS only up to the first time the failed set is reached, at which point the rest of the cycle reverts to crude simulation. For the IS, we employ *balanced failure biasing* [16], which, in a state from which both failures and repairs are possible, assigns a total probability $ρ = 0.8$ to all failure transitions, equally allocating $ρ$ across the individual failure transitions; repair transitions get a total probability of $1 - ρ$, which is distributed to individual repair transitions in proportion to their original probabilities. [16] establishes that this dynamic IS procedure yields an estimator of $θ_{(r),1}$ with BRE. We utilize crude simulation to estimate $θ_{(r),2}$. Table 2 illustrates again that BRE is satisfied for all estimators.

<table><caption>Table 2: Results for unavailability computation for the HRMS of Example 10, from $n = 10^6$ independent cycles.</caption><thead><tr><th>ε</th><th>θ̂<sub>(r),n1,1</sub></th><th>R̂E<sup>(1)</sup></th><th>θ̂<sub>(r),n2,2</sub></th><th>R̂E<sup>(2)</sup></th><th>α̂<sub>(r),n</sub></th><th>R̂E[α̂<sub>(r),n</sub>]</th><th>95% conf. int.</th></tr></thead><tbody><tr><td>1.00e-01</td><td>2.90e-03</td><td>8.85e+01</td><td>2.78e+00</td><td>8.85e-01</td><td>1.04e-03</td><td>9.33e+01</td><td>(8.516e-04, 1.233e-03)</td></tr><tr><td>1.00e-02</td><td>1.11e-06</td><td>5.62e+00</td><td>7.74e+00</td><td>7.29e-02</td><td>1.43e-07</td><td>5.93e+00</td><td>(1.414e-07, 1.447e-07)</td></tr><tr><td>1.00e-03</td><td>1.01e-09</td><td>6.98e+00</td><td>6.77e+01</td><td>2.55e-03</td><td>1.49e-11</td><td>7.35e+00</td><td>(1.467e-11, 1.510e-11)</td></tr><tr><td>1.00e-04</td><td>1.00e-12</td><td>7.17e+00</td><td>6.68e+02</td><td>8.45e-05</td><td>1.50e-15</td><td>7.56e+00</td><td>(1.476e-15, 1.521e-15)</td></tr><tr><td>1.00e-05</td><td>9.98e-16</td><td>7.20e+00</td><td>6.67e+03</td><td>3.46e-06</td><td>1.50e-19</td><td>7.59e+00</td><td>(1.475e-19, 1.519e-19)</td></tr><tr><td>1.00e-06</td><td>9.98e-19</td><td>7.20e+00</td><td>6.67e+04</td><td>1.16e-06</td><td>1.50e-23</td><td>7.59e+00</td><td>(1.475e-23, 1.519e-23)</td></tr></tbody></table>

The literature provides many other numerical examples displaying empirical behavior consistent with the theory in Theorem 1. For the MTTF ratio estimator in Example 5, [5] show the BRE
property through numerical results on a model of a computing system. When estimating the mean time to hit a large buffer size *r* in a stable M/M/1 queue, [11] present simulation results using the ratio estimator in Examples 2 and 8 and the alternative estimator described in Example 7, with the latter approach outperforming the former for heavier congestion levels. The M/M/1 estimators have LE but not BRE.

# 5 Concluding Remarks

We considered estimating a function of a vector of means in a rare-event context. This covers a wide range of estimands that arise in practice, as we illustrated through numerous examples. Theorem 1 provides conditions that show when the overall estimator can be estimated efficiently when each of the individual means are.

# References

[1] S. Asmussen and P. Glynn. *Stochastic Simulation: Algorithms and Analysis*. Springer, New York, 2007.

[2] P. Billingsley. *Probability and Measure*. John Wiley and Sons, New York, 3rd edition, 1995.

[3] P. W. Glynn. Likelihood ratio derivative estimators for stochastic systems. *Communications of the ACM*, 33:75–84, 1990.

[4] P. W. Glynn, M. K. Nakayama, and B. Tuffin. On the estimation of the mean time to failure by simulation. In W. K. V. Chan, A.D'Ambrogio, G. Zacharewicz, N. Mustafee, G. Wainer, and E. Page, editors, *Proceedings of the 2017 Winter Simulation Conference*, pages 1844–1855, Piscataway, NJ, 2017. Institute of Electrical and Electronics Engineers.

[5] A. Goyal, P. Shahabuddin, P. Heidelberger, V. Nicola, and P. W. Glynn. A unified framework for simulating Markovian models of highly dependable systems. *IEEE Transactions on Computers*, C-41(1):36–51, 1992.

[6] V. Kalashnikov. *Topics on Regenerative Processes*. CRC Press, Boca Raton, 1994.

[7] P. L'Ecuyer, J. H. Blanchet, B. Tuffin, and P. W. Glynn. Asymptotic robustness of estimators in rare-event simulation. *ACM Transactions on Modeling and Computer Simulation*, 20(1):Article 6, 2010.

[8] P. L'Ecuyer and B. Tuffin. Approximating zero-variance importance sampling in a reliability setting. *Annals of Operations Research*, 189(1):277–297, 2012.

[9] M. K. Nakayama. General conditions for bounded relative error in simulations of highly reliable Markovian systems. *Advances in Applied Probability*, 28:687–727, 1996.

[10] M. K. Nakayama. On derivative estimation of the mean time to failure in simulations of highly reliable Markovian systems. *Operations Research*, 46:285–290, 1998.

[11] M. K. Nakayama and B. Tuffin. Efficient estimation of the mean hitting time to a set of a regenerative system. In N. Mustafee, K.-H.G. Bae, S. Lazarova-Molnar, M. Rabe, C. Szabo, P. Haas, and Y.-J. Son, editors, *Proceedings of the 2019 Winter Simulation Conference*, pages 416–427, Piscataway, New Jersey, 2019. Institute of Electrical and Electronics Engineers.

[12] M. K. Nakayama and B. Tuffin. Density estimators of the cumulative reward up to a hitting time to ararely visited set of a regenerative system. In B. Feng, G. Pedrielli, Y. Peng, S. Shashaani, E. Song, C.G. Corlu, L.H. Lee, E.P. Chew, and T. Roeder and P. Lendermann, editors, *Proceedings of the 2022 Winter Simulation Conference*, pages 121–132, Piscataway, NJ, 2022. Institute of Electrical and Electronics Engineers.

[13] G. Rubino and B. Tuffin, editors. *Rare Event Simulation using Monte Carlo Methods*. John Wiley, Chichester, UK, 2009.
[14] J. S. Sadowsky. Large deviations theory and efficient simulation of excessive backlogs in a gi/gi/m queue. *IEEE Transactions on Automatic Control*, 36:1383–1394, 1991.

[15] R. J. Serfling. *Approximation Theorems of Mathematical Statistics*. John Wiley and Sons, New York, 1980.

[16] P. Shahabuddin. Importance sampling for highly reliable Markovian systems. *Management Science*, 40:333–352, 1994.