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.
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.
