• $q$ be a probability distribution on $\mathcal{X}$,
  • $w$ be a nonnegative function from $\mathcal{X}$ to $\mathbf{R}$ which is bounded away from $0$ and $\infty$, and
  • $s$ be a bounded function from $\mathcal{X}$ to $\mathbf{R}_+$ with mean $0$ under $q$

I am interested in unbiased estimation of the following quantity:

\begin{align} \gamma = \int_{\mathcal{X}^N} \left( \prod_{i = 1}^N q (x_i)\right) \cdot \frac{w(x_1)}{\sum_{i = 1}^N w(x_i)} \cdot s(x_1) \, dx_{1:N} \end{align}

given $(N + 1)$ iid samples from $q$, where $N$ is a positive integer.

Now, given that $s$ has mean $0$ under $q$, I can actually rewrite this expression in terms of a covariance, i.e. define

\begin{align} W ( x_1 | x_{2:N} ) = \frac{w(x_1)}{\sum_{i = 1}^N w(x_i)}, \end{align}

then \begin{align} \gamma = \int_{\mathcal{X}^{N - 1}} \left( \prod_{i = 2}^N q (x_i)\right) \cdot \text{Cov}_{q(x_1)} \left( W ( x_1 | x_{2:N} ), s(x_1) \right) \, dx_{2:N} \end{align}

Given samples $x_{0:N}$, one can then form an unbiased estimator of the inner covariance as

\begin{align} \hat{c} (x_0, x_1 | x_{2:N} ) \triangleq \frac{1}{2} \left\{ W ( x_0 | x_{2:N} ) - W ( x_1 | x_{2:N} ) \right\} \cdot \left\{ s(x_0) - s(x_1) \right\}, \end{align}

and so one can form an overall estimator by re-using this estimator for all splits of the data, i.e.

\begin{align} \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ). \end{align}

Unfortunately, naive evaluation of the above sum costs $\mathcal{O} \left( N^2 \right)$, which is more expensive than I would like. My question is thus: can one simplify this sum such that it can be evaluated exactly in time $\mathcal{O} \left( N \right)$?

One reason why I suspect that this may be possible is that $\hat{c}$ can be re-written as

\begin{align} \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ) = \frac{1}{2} \left\{ \frac{w (x_i)}{W - w (x_j)} - \frac{w (x_j)}{W - w (x_i)} \right\} \cdot \left\{ s(x_i) - s(x_j) \right\}, \end{align}

where $W = \sum_{i = 0}^N w (x_i)$, which may be able to simplify some things.

  • $\begingroup$ The term $W(x_1|x_{2:N})\times s(x_1)$ is an unbiased estimator of $\gamma$. $\endgroup$
    – Xi'an
    Nov 26, 2020 at 18:28
  • $\begingroup$ @Xi'an indeed! My aim is not just to construct an unbiased estimator, but to construct one which i) uses the knowledge that $s$ is mean-0 to reduce variance, ii) uses the full information of all of the samples which I have generated, and iii) can be evaluated quickly. The single-term estimator satisfies only iii), as it has cost $\mathcal{O} \left( N \right)$. A mid-way point would be to use $$ \frac{1}{N+1} \sum_{i = 0}^N \hat{c} (x_i, x_{i+1} | x_{0:N \setminus \{i, i + 1 \} } )$$ (interpreting indices $\mod (N+1)$), which also costs $\mathcal{O} \left( N \right)$ to form. $\endgroup$
    – πr8
    Nov 26, 2020 at 19:13

1 Answer 1


I think that the following should work: re-write

\begin{align} \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ) &= \frac{1}{2} \left\{ \frac{w (x_i)}{W - w (x_j)} - \frac{w (x_j)}{W - w (x_i)} \right\} \cdot \left\{ s(x_i) - s(x_j) \right\} \\ &= \frac{1}{2} \left\{ \frac{1}{W - w (x_j)} \cdot \left[ w(x_i) \cdot s(x_i) \right] + \frac{1}{W - w (x_i)} \cdot \left[ w(x_j) \cdot s(x_j) \right] \right\} \\ &- \frac{1}{2} \left\{ \frac{s(x_j)}{W - w (x_j)} \cdot w(x_i) + \frac{s(x_i)}{W - w (x_i)} \cdot w(x_j) \right\} \end{align}

which allows for the reformulation

\begin{align} \hat{\gamma} &\triangleq \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ) \\ &= \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \left\{ \frac{1}{W - w (x_j)} \cdot \left[ w(x_i) \cdot s(x_i) \right] \right\} \\ &- \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \left\{ \frac{s(x_j)}{W - w (x_j)} \cdot w(x_i) \right\} \\ &= \frac{N+1}{N} \cdot \left\{ \overline{\frac{1}{W - w}} \cdot \overline{ws} - \overline{\frac{s}{W-w}} \cdot \overline{w} \right\} \end{align}


\begin{align} \overline{\frac{1}{W - w}} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N \frac{1}{W - w (x_i) } \\ \overline{ws} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N w(x_i) s(x_i) \\ \overline{\frac{s}{W-w}} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N \frac{s(x_i)}{W - w(x_i)} \\ \overline{w} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N w(x_i) \end{align}

which suggests that the whole estimator can be computed in time $\approx 4N = \mathcal{O} \left( N \right)$, noting that $\overline{w} = \frac{W}{N + 1}$ does not need to be re-computed.

Some additional manipulations reveal that this estimator can be written as

\begin{align} \hat{\gamma} &= \left\{ \frac{1}{N+1} \sum_{i = 0}^N \frac{1}{W - w(x_i)} \right\} \cdot \left\{ \frac{1}{N} \sum_{i = 0}^N \left( s(x_i) - \overline{s} \right) \cdot \left( w(x_i) - \overline{w} \right) \right\} \\ &- \left\{ \frac{1}{N} \sum_{i = 0}^N \left( \frac{1}{W - w(x_i)} - \overline{\frac{1}{W - w}} \right) \cdot \left( s(x_i) - \overline{s} \right) \right\} \cdot \left\{ \frac{1}{N+1} \sum_{i = 0}^N w(x_i) \right\} \end{align}

which may suggest implementations of this estimator which achieve improved numerical stability.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.