Typically, in Bayesian bootstrap, you have samples {$x_1,…,x_n$} of a random variable $X$.
Choose $\{p_1,…,p_n\}$ from a Dirichlet distribution, by sorting $\{0,1,u_1,…,u_{n-1}\}$ where $u_i$ are independent samples from $U \sim Uniform(0,1)$, and taking differences of sorted neighbors.
Use $\{p_1,…,p_n\}$ as weights for $\{x_1,…,x_n\}$ to calculate a sample from the statistic of interest $\theta$; i.e.
$$\hat\theta_j = \sum_{i=1}^n f(x_i) \cdot p_i$$
for bootstrap iteration $j$.
I already have weights $\{w_1,…,w_n\}$ from importance sampling, from which I compute
$$\hat\theta = \frac{\sum_{i=1}^n f(x_i) \cdot w_i}{\sum_{i=1}^n w_i} $$
With a new set of weights computed as $w_i' = w_i \cdot p_i$, computing
$$\hat\theta_j = \frac{\sum_{i=1}^n f(x_i) \cdot w_i'}{\sum_{i=1}^n w_i'} $$
doesn't work: the distribution of bootstrapped statistics is too narrow.
Resampling from the weighted samples to get unweighted $\{x_1',…,x_n'\}$ and applying Bayesian bootstrap to these also doesn't work: the distribution of bootstrapped statistics is too wide. (This is also true if I resample anew for each Bayesian bootstrap iteration.)
I think I need to reweight differently. How is this done?
EDIT: Actually, using weights $w_i' = w_i \cdot p_i$ was the right solution. I was comparing my bootstrap results against the results of using a bad, too-high Monte Carlo variance estimate. I'm not ready to call this foobar my fault, though, because it's stupidly hard to find the correct formula for Monte Carlo variance of expected values computed using "ratio" or "self-normalized" importance samples. For the record, I finally found it here, on page 63:
http://www.rni.helsinki.fi/~pek/test/mcmbc10.pdf
So thanks very much to Petri Koistinen. This is the formula for an unbiased Monte Carlo variance estimate, in case the link goes stale:
$$Var[\widehat\theta] \approx \frac{n}{n-1} \frac{\sum_{i=1}^n \left(f(x_i)-\widehat\theta\right)^2 w_i^2}{\left(\sum_{i=1}^n w_i\right)^2}$$
I think this is right because 1) Monte Carlo variance of expected values computed using unweighted samples is a special case; 2) changing the importance distribution induces the expected changes to it; 3) all the bootstrap methods I've implemented agree with the values I compute from it (finally); 4) it's invariant to weight scaling; and 5) it seems to have a decent proof behind it, in Koistinen's notes and another place I found it. (Though in the other place it was followed up by a bogus formula for computing a confidence interval, with a spurious multiplication by $\sqrt{n}$.) As a bonus, the leading term is the standard bias correction term for variance estimates.
You would think that, given that importance sampling is so often presented and used as a variance reduction technique, estimating the variance would be common.
Sorry for the noise, but at least it might help someone else in the future.
Best Answer
Think of the importance sampling as defining a new $f^*(x_i) = f(x_i)w_i/\sum w_i$. Then, when you do your Bayesian bootstrap, as in your formula for $\hat{\theta}_j$, do it with respect to $f^*$ instead of the original $f$.
The concept behind this is that the importance sampling "corrects" for a misspecified probability of observing the $x_i$ (misspecified because it isn't $1/n$.) This correction should be the same regardless of the weights generated for a particular bootstrap sample, because it's due to the underlying probability model. So you shouldn't combine the correction weights with the bootstrap weights to form new weights (your first approach). On the other hand, your second approach adds randomness by sampling according to the importance sampling weights rather than doing the importance sampling calculation, thus adding randomness to your results and causing the resulting distribution to be too wide. Your third approach adds randomness for the same reason, but to a different degree, and therefore has the same effect.