Looking for a recursive formula for asymptotic variance of importance sampling estimator (self-normalized) Looking for a recursive formula to approximate variance of importance sampling estimator $Var_q\big[\delta_{IS}\big]\approx\sum_{i=1}^n\tilde w(X_i)^2\big[h(X_i)-\delta_{IS}\big]^2$. This is an approximation of $$Var\big[\delta_{IS}\big]=E_q\Bigg\{w(X)^2\Big[h(X)-E_q\big[h(X)\big]\Big]^2\Bigg\}$$
I'm using this to make running variance plot.
 A: Define
\begin{align}
w_i &= w \left( X_i \right) \\
h_i &= h \left( X_i \right) \\
\mu^w_N &= \frac{1}{N} \sum_{i = 1}^N w_i \\
\mu^{wh}_N &= \frac{1}{N} \sum_{i = 1}^N w_i h_i \\
\tilde{\mu}^h_N &= \frac{\mu^{wh}_N}{\mu^w_N}
\end{align}
so that $\tilde{\mu}^h_N$ is the standard SNIS estimator.
The standard estimator of the asymptotic variance for SNIS is then given by
\begin{align}
V_N &= \frac{\frac{1}{N}\sum_{i = 1}^N \left\{  w_i^2 \cdot \left( h_i - \tilde{\mu}^h_N \right)^2 \right\} }{\left( \frac{1}{N} \sum_{i = 1}^N w_i \right)^2}
\end{align}
Expand this to see that
\begin{align}
V_N &= \frac{\left( \frac{1}{N} \sum_{i = 1}^N w_i^2 h_i^2 \right) - 2 \cdot \tilde{\mu}^h_N \cdot \left( \frac{1}{N} \sum_{i = 1}^N w_i^2 h_i \right) + \left( \tilde{\mu}^h_N \right)^2 \cdot \left( \frac{1}{N} \sum_{i = 1}^N w_i^2 \right)}{\left( \frac{1}{N} \sum_{i = 1}^N w_i^2 \right)} \\
&= \frac{\mu^{wwhh}_N - 2 \tilde{\mu}^h_N \cdot \mu^{wwh}_N + \left( \tilde{\mu}^h_N \right)^2 \cdot \mu^{ww}_N}{\left( \mu^{w}_N \right)^2} \\
&= \frac{\left( \mu^w_N \right)^2 \cdot \mu^{wwhh}_N - 2 \mu^w_N \cdot \mu^{wh}_N \cdot \mu^{wwh}_N + \left( \mu^{wh}_N \right)^2 \cdot \mu^{ww}_N}{ \left( \mu^w \right)^4}
\end{align}
where
\begin{align}
\mu^{ww}_N &= \frac{1}{N} \sum_{i = 1}^N w_i^2 \\
\mu^{wwh}_N &= \frac{1}{N} \sum_{i = 1}^N w_i^2 h_i \\
\mu^{wwhh}_N &= \frac{1}{N} \sum_{i = 1}^N w_i^2 h_i ^2.
\end{align}
As such, the estimator $V_N$ can be computed as a function of the statistics $\mu^w_N, \mu^{ww}_N, \mu^{wh}_N, \mu^{wwh}_N, \mu^{wwhh}_N$. Each of these statistics is a sample average, which can easily be updated recursively.
A modification of the above procedure which may provide improved numerical stability would be to note that
\begin{align}
\mu^{ww}_N &= \left( \mu^{w}_N \right)^2  + \frac{1}{N} \sum_{i = 1}^N \left( w_i - \mu^{w}_N  \right)^2 \\
\mu^{wwh}_N  &= \mu^{wwh}_N  \cdot \mu^{wh}_N  + \frac{1}{N} \sum_{i = 1}^N \left( w_i - \mu^{ww}_N  \right) \left(w_i h_i - \mu^{wh}_N \right) \\
\mu^{wwhh}_N &= \left( \mu^{wh}_N \right)^2 + \frac{1}{N} \sum_{i = 1}^N \left( w_i h_i - \mu^{wh}_N \right)^2,
\end{align}
define
\begin{align}
\beta^w_N &= \frac{1}{N - 1} \sum_{i = 1}^N \left( w_i - \mu^{w}_N \right)^2 \\
\rho^{w, wh}_N &= \frac{1}{N - 1} \sum_{i = 1}^N \left( w_i - \mu^{w}_N  \right) \left(w_i h_i - \mu^{wh}_N  \right) \\
\beta^{wh}_N &= \frac{1}{N - 1} \sum_{i = 1}^N \left( w_i h_i - \mu^{wh}_N  \right)^2,
\end{align}
and use standard, stable, recursive update formulas for variances and covariances to update $\left(\beta^w_N, \rho^{w, wh}_N, \beta^{wh}_N \right)$ in addition to $\left( \mu^{w}_N, \mu^{wh}_N \right)$. These estimators are all unbiased for their corresponding ideal quantities and are otherwise well-behaved, and so may be useful for other diagnostic purposes.
