I suspect this is a fairly unusual and exploratory question, so please bear with me.
I am wondering if one could apply the idea of importance sampling to Gibbs sampling. Here's what I mean: in Gibbs sampling, we change the value of one variable (or block of variables) at a time, sampling from the conditional probability given the remaining variables.
However, it may not be possible or easy to sample from the exact conditional probability. So instead we sample from a proposal distribution $q$ and use, for example, Metropolis-Hastings (MH).
So far, so good. But here is a divergent path: what happens if we, instead of using MH, use the same idea used in importance sampling, namely we sample from $q$ and keep an importance weight $p/q$ of the current sample?
In more detail: assume we have variables $x_1,\dots,x_n$ and a factored distribution $\phi_1,\dots,\phi_m$ so that $p \propto \prod_{i=1}^m \phi_i$. We keep the proposal probability $q_i$ used to sample the current value of each variable $x_i$. At each step we change a subset of the variables and update $p(x)/q(x)$ (only the factors of $p$ and $q$ that are affected). We take the samples and their importance weight to compute whatever statistic we are interested in.
Would this algorithm be correct? If not, any clear reasons why not? Intuitively it makes sense to me as it seems to be doing the same thing importance sampling does, but with dependent samples instead.
I did implement this for a Gaussian random walk model and observed that the weights become smaller and smaller (but not monotonically), so the initial samples end up having too much importance and dominate the statistic. I'm pretty certain the implementation is not buggy, because at each step I compare the updated weight to an explicit brute-force computation of it. Note that the weights do not go down indefinitely to zero, because they are $p/q$ where both $p$ and $q$ are products of a finite number of densities, and each sample is obtained from a Normal distribution that only rarely will be zero.
So I am trying to understand why the weights go down like that, and whether this is a consequence of this method being actually not correct.
Here's a more precise definition of the algorithm, applied to a Gaussian random walk on variables $X_1,\dots,X_n$. The code follows below.
The model is simply $X_i \sim \mathcal N(X_{i-1}, \sigma^2), i = 1,\dots,n$, with $X_0$ fixed to $0$.
The weight of the current sample is $\frac{\prod_i p(x_i)}{\prod_i q(x_i)}$, where $p$ are the Gaussian densities and $q$ are the distributions from which the current values have been sampled. Initially, we simply sample the values in a forward manner, so $q = p$ and the initial weight is $1$.
Then at each step I choose $j \in \{1,\dots,n\}$ to change. I sample a new value $x'_j$ for $X_j$ from $\mathcal N(X_{j-1},\sigma^2)$, so this density becomes the new used proposal distribution for $X_j$.
To update the weight, I divide it by the densities $p(x_j | x_{j-1})$ and $p(x_{j+1} | x_j)$ of old value $x_j$ according to $x_{j-1}$ and $x_{j+1}$, and multiply by the densities $p(x'_j | x_{j-1})$ and $p(x_{j+1} | x'_j)$ of new value $x'_j$ according to $x_{j-1}$ and $x_{j+1}$. This updates the numerator $p$ of the weight.
To update the denominator $q$, I multiply the weight by the old proposal $q(x_j)$ (thus removing it from the denominator) and divide it by $q(x'_j)$.
(Because I sample $x'_j$ from the normal centered on $x_{j-1}$, $q(x'_j)$ is always equal to $p(x'_j | x_{j-1})$ so they cancel out and the implementation does not actually use them).
Like I mentioned before, in the code I compare this incremental weight computation to the actual explicit computation just to be sure.
Here's the code for reference.
println("Original sample: " + currentSample);
int flippedVariablesIndex = 1 + getRandom().nextInt(getVariables().size() - 1);
println("Flipping: " + flippedVariablesIndex);
double oldValue = getValue(currentSample, flippedVariablesIndex);
NormalDistribution normalFromBack = getNormalDistribution(getValue(currentSample, flippedVariablesIndex - 1));
double previousP = normalFromBack.density(oldValue);
double newValue = normalFromBack.sample();
currentSample.set(getVariable(flippedVariablesIndex), newValue);
double previousQ = fromVariableToQ.get(getVariable(flippedVariablesIndex));
fromVariableToQ.put(getVariable(flippedVariablesIndex), normalFromBack.density(newValue));
if (flippedVariablesIndex < length - 1) {
NormalDistribution normal = getNormalDistribution(getValue(currentSample, flippedVariablesIndex + 1));
double oldForwardPotential = normal.density(oldValue);
double newForwardPotential = normal.density(newValue);
// println("Removing old forward potential " + oldForwardPotential);
currentSample.removePotential(new DoublePotential(oldForwardPotential));
// println("Multiplying new forward potential " + newForwardPotential);
currentSample.updatePotential(new DoublePotential(newForwardPotential));
}
// println("Removing old backward potential " + previousP);
currentSample.removePotential(new DoublePotential(previousP));
// println("Multiplying (removing from divisor) old q " + previousQ);
currentSample.updatePotential(new DoublePotential(previousQ));
println("Final sample: " + currentSample);
println();
// check by comparison to brute force calculation of weight:
double productOfPs = 1.0;
for (int i = 1; i != length; i++) {
productOfPs *= getNormalDistribution(getValue(currentSample, i - 1)).density(getValue(currentSample, i));
}
double productOfQs = Util.fold(fromVariableToQ.values(), (p1, p2) -> p1*p2, 1.0);
double weight = productOfPs/productOfQs;
if (Math.abs(weight - currentSample.getPotential().doubleValue()) > 0.0000001) {
println("Error in weight calculation");
System.exit(0);
}