Running Metropolis-Hastings algorithm with changing proposal kernel; each time the kernel is changing starting the algorithm afresh. Does it work?

Question

I have a Markov kernel $Q$ from which I would like to generate proposals for the Metropolis-Hastings algorithm. The problem is: When the proposal is accepted, the "internal state" of $Q$ changes. This means that if proposals $y_1,\ldots,y_n$ are accepted, the internal state of $Q$ depends on $y_1,\ldots,y_n$. I know, this means that we cannot use $Q$ as a proposal kernel for the Metropolis-Hastings algorithm.

However, my simple solution to that problem is the following: Before the first sample is accepted, the state of $Q$ does not change. Now, I simply run the Metropolis-Hastings algorithm with proposal kernel $Q$ until the first proposal is accepted. Then I stop. Then I start the Metropolis-Hastings algorithm again, but with the different proposal kernel given by the modified kernel $Q$.

Is this process still guaranteed to work? Are the accepted samples distributed according to the target density after a sufficient long period of time?

EDIT:

I think we can describe the algorithm I've got in mind as follows:

• Let $E$ be the state space and $Q_k$ be a Markov kernel with source $E^k$ and target $E$
• Start with any $x_0\in E$
• Run Metropolis-Hastings with initial state $x_0$ and proposal kernel $Q_1$ for a single iteration
• Let $y_1\in E$ denote the proposed sample and $x_2\in E$ the state after the iteration (So, $x_1=y_1$ if the proposal was accepted and $x_1=x_0$ otherwise)
• Now run Metropolis-Hastings with initial state $x_1$ and proposal kernel $Q_2(x_0,y_1,\;\cdot\;)$ (remark: I'm unsure whether it wouldn't be better to replace this with $Q_2(x_0,x_1,\;\cdot\;)$)
• and so on ...

It would be interesting to know whether - under certain assumptions on $Q_1,Q_2,\ldots$ - the samples $x_b,x_{b+1},\ldots$ are still distributed according to the target density for sufficiently large $b$.

EDIT 2

You can assume that $Q_k(x_1,\ldots,x_k;\;\cdot\;)$ has density $$q_k(x_1,\ldots,x_k;\;\cdot\;):=e^{-\beta f_k(x_1,\ldots,x_k;\;\cdot\;)}$$ with respect to the Lebesgue measure on $[0,1)^d$, where $f_k(x_1,\ldots,x_k;\;\cdot\;)$ is nonnegative. Also: The $Q_i$ are constructed in a way so that at the last iteration $k_{\text{max}}$, we have $f_{k_{\text{max}}}(x_1,\ldots,x_{k_{\text{max}}})=0$. So, the sequence $Q_1,\ldots,Q_{k_{\text{max}}}$ somehow converges; maybe this is enough to show that everything works.

Answer 1

If I understand the proposed algorithm correctly, we can prove that this doesn't generally sample from the target distribution by way of counter example. And while this algorithm will be a bit pathological, I believe it does do a good job of illustrating what can easily go wrong with this method.

Consider the case that the target distribution is a N(0, 1) distribution. For $Q_1$, we select a N(0, 0.5) distribution. For $Q_2$, we select a N(0, 0.001) if $x_1$ was positive and N(0, 0.5) if $x_2$.

In the target distribution, 50% of the values should be positive. However, in our algorithm, if the current state becomes positive, it gets stuck in the positive region for a very long time due to the small step size. If the current state is negative, it moves quickly due to the large step size. Thus, the samples from our algorithm will be disporpotionally positive, even after sufficient burnin. Therefore, we never approach the target distribution.

Answer 2

Given the large amount of comments I thought that it might be better when I place some of them in an answer.

What I understand from the comments following the question is that the idea is about performing some form of MCMC-sampling, but with a Kernel that adapts after each new sample. And, the motivation is to ensure that the sample will satisfy a certain condition.

Simple example

A simple example would be the sampling of a distribution symmetric around zero (e.g. a standard normal distribution) but with the constraint that the sample mean needs to be zero (typically, if the distribution has zero mean, a sample from it does not need to be zero).

If we have a kernel that proposes every time that the sample mean is unequal to zero a sample whose value makes the sample mean equal to zero. Then we ensure that this property is fulfilled (at least every odd step).

The set of samples that can be sampled will not be iid variables and will not be the same as sampling the target distribution without the constraint. However, if the target distribution is symmetric around zero, then this algorithm will generate a sample whose empirical distribution approaches the true distribution.

Below is a code example that has the proposed sample based on the complete history $$x^\star|x_0,x_1,\dots,x_t \sim \begin{cases} -\sum_{i=0}^t x_i & \qquad \text{if $\sum_{i=0}^t x_i \neq 0$} \\ N(x_t,0.04) & \qquad \text{if $\sum_{i=0}^t x_i = 0$} \\ \end{cases}$$

Below is an example of the histogram of a sample of size 50000 when the target function is a standard normal.

This sample is in not a typical sample from a normal distribution. The sample will be having zero mean with probability 1, whereas a sample from a normal distribution will be having zero mean with probability 0 (and also the samples will be relatively symmetric).

However, the empirical distribution will approach the the distribution function of the target distribution. So in that sense this sampling 'works'.

For different more complex cases it will depend. For example, when we use the method above with a non-symmetric distribution, then it stops 'working'.

set.seed(1)

newsample = function(old_sample, LikelihoodFunction) {
    L = length(old_sample)
    m = sum(old_sample)

    ### if the current sample has not zero mean then the suggestion will always be a sample that makes the mean zero
    if (m!=0) {
      suggest = -m
    } else {
      suggest = rnorm(1,old_sample[L],0.2)
    }

    ### compare the likelihood and base the next sample on it
    u = runif(1)
    l1 = LikelihoodFunction(suggest)
    l2 = LikelihoodFunction(old_sample[L])
    if (u<l1/l2) {
      return(suggest)
    } else {
      return(old_sample[L])
    }    
}

### start the mcmc in the point 0
sample = c(0)

### generate 50000 samples
for (i in 1:50000) {
   sample = c(sample, newsample(sample, LikelihoodFunction = function(x) {dnorm(x)}))
}

### plot histogram with curve for target density as comparison 
hist(sample, seq(-5,5,0.1), xlim =c(-3,3), freq = 0)
xs = seq(-4,4,0.01)
lines(xs, dnorm(xs))