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

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.

• I'd be interested to know how this problem came about! Apr 3 at 20:02
• @jbowman It is important for me that the spectrum of the generated proposals have a certain profile. In order to ensure that, I need to minimize some kind of "energy" induced by the previously generated points and the new proposal. (But the locations of the previously generated samples are not altered during this minimization step.) Apr 3 at 20:10
• So, if I've followed you, which I may well not have, you are updating the proposal as you learn more about the distribution - which would, in the long run, lead to higher, perhaps much higher, acceptance probabilities? Apr 3 at 21:02
• @jbowman "you are updating the proposal as you learn more about the distribution" ... The target distribtion? No. It is actually quite simple: For a set of $k$ vectors $x_1,\ldots,x_k$, we have an energy $E(x_1,\ldots,x_k)$. The proposal scheme works as follows: Given $x_0$, generate new proposals according to $e^{-E(x_0,\;\cdot\;)/T}$, where $T>0$. Once a proposal is accepted, call it $x_1$, we formally restart MH, but this time we generate new proposals according to $e^{-E(x_0,\:x_1,\;\cdot\;)/T}$. And so on. Apr 3 at 21:28
• @jbowman The question is: Does this work (i.e. can I build asymptotically unbiased estimators from the generated samples)? Maybe projecteuclid.org/journals/bernoulli/volume-7/issue-2/… is the answer to that question, but I'm not 100% sure. Apr 3 at 21:31

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.

• Thank you for your answer. Did you take a look into the comments? Does the situation change with the specific form of the proposal? May 6 at 11:14
• Also: Doesn't this paper describe such an adaptive Metropolis-Hastings variant? It seems like that we might be able to show that the process obtained still has Ergodic properties; even when it is not Markov anymore. May 6 at 11:26
• @0xbadf00d: indeed adaptive MCMC is a valid approach. However, in the paper you've cited, they show that the proposal distribution converges which is key. In my pathological case, it does not. If I'm understanding your proposal, it doesn't seem like the proposal distribution will converge, in fact it will different conditional on the current state? This will require some non-trivial conditions to make work, although not impossible (e.g. Hamiltonian Monte Carlo). May 6 at 11:35
• Can you clarify what exactly you mean by claiming that the "proposal distribution does converge"? We have a sequence of proposal kernels $Q_k$; do you mean that there must be a kernel $Q$ so that $Q_k(x,\;\cdot\;)$ converges to $Q(x,\;\cdot\;)$ for all $x$? If so, in which sense must the convergence take place (weakly, total variation, ...)? May 6 at 12:02
• In the paper you cited, they mention how the covariance matrix eventually converges essentially as an approximation of the target function. If I understand your algorithm, the proposal function is dependent on the current state and we should not generally expect it to converge unless there are other stipulations I am missing? May 6 at 12:12

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