# Metropolis-Hastings: target distribution with two modes; deterministic transformation

I'm trying to construct a Metropolis-Hastings algorithm to sample a target distribution $$p(x)$$ with two different and isolated modes. The example I'm working with is

$$$$p(x) = \frac{\left(e^{- \frac{1}{2} \left(x + 4.5\right)^{2}} + e^{- \frac{1}{2} \left(x - 14.5\right)^{2}}\right)}{2\sqrt{2\pi}},$$$$

which are two Gaussians centered at $$\mu = -4.5$$ and $$\nu = 14.5$$.

Note: in the image above the variable is called $$p$$ instead of $$x$$, and the probability distribution is $$f(p)$$ instead of $$p(x)$$.

If I wanted to sample this using a proposal distribution $$q(x'|x)$$ that is just the previous accepted state plus a uniform random number between -1 and 1,

$$$$x' = x + \epsilon, \qquad \text{ with } \epsilon \in [-1,1],$$$$

with $$\epsilon$$ a random number, then the algorithm would get stuck in one of the Gaussian distributions.

What I want to do is to modify this by doing a deterministic transformation at some steps that brings the system from one Gaussian to another. My problem is that it's not working (I'm not recovering the statistical momenta of the target distribution), and I don't know if I got something wrong conceptually or it is just an error of the code.

The transformation I want to use is

\begin{align} x' &= \nu + (-x-\nu)e^{-t} \qquad \text{if } x<0 \\ x' &= \mu + (-x-\mu)e^{-t} \qquad \text{if } x>0, \end{align}

for some $$t \geq 0$$. Note that for $$t=0$$ this is just a reflection with respect to $$x=0$$: $$x' = -x$$; and for $$t \rightarrow \infty$$, this brings the system to the peak of the other Gaussian: if $$x<0$$, then $$x' = \nu$$ (and if $$x>0$$, then $$x'=\mu$$). The $$t$$ is selected so as to maximize the acceptance probability,

\begin{align} A(x'|x) = \min\left\{1, \frac{p(x')}{p(x)} \frac{q(x|x')}{q(x'|x)} |J(x\rightarrow x')| \right\}, \end{align}

where $$|J(x\rightarrow x'|$$ is the Jacobian of the transformation from $$x$$ to $$x'$$. But for this to sample the target distribution $$p(x)$$ correctly, one has to add the inverse transformation as well,

\begin{align} x' &= -\nu + (-x+\nu)e^{t} \qquad \text{if } x>0 \\ x' &= -\mu + (-x+\mu)e^{t} \qquad \text{if } x<0, \end{align}

and this way the acceptance probability reduces to

\begin{align} A(x'|x) = \min\left\{1, \frac{p(x')}{p(x)} |J(x\rightarrow x')| \right\}. \end{align}

### The algorithm

So, the algorithm would be:

1. Set an initial state $$x_0$$ as the current state, $$x = x_0$$.

2a. Propose a new state $$x' = x + \epsilon$$ and accept it with probability

\begin{align} A(x'|x) = \min\left\{1, \frac{p(x')}{p(x)} \right\}. \end{align} 2b. Once every $$n$$ steps, instead of step 2a, choose randomly between the "forward" transformation

\begin{align} x' &= \nu + (-x-\nu)e^{-t} \qquad \text{if } x<0 \\ x' &= \mu + (-x-\mu)e^{-t} \qquad \text{if } x>0, \end{align}

or the inverse transformation

\begin{align} x' &= -\nu + (-x+\nu)e^{t} \qquad \text{if } x>0 \\ x' &= -\mu + (-x+\mu)e^{t} \qquad \text{if } x<0. \end{align}

and accept it with probability

\begin{align} A(x'|x) = \min\left\{1, \frac{p(x')}{p(x)} |J(x\rightarrow x')| \right\}. \end{align}

with $$|J(x\rightarrow x')| = e^{-t}$$ for the "forward" transformation and $$|J(x\rightarrow x')| = e^{t}$$ for the inverse transformation. I've checked that the value of $$t$$ that maximizes the acceptance probability is around $$t=2.2$$.

1. If accepted, set $$x'$$ as the new current accepted state, $$x = x'$$. Otherwise, set the previous state as the new current state, $$x=x$$.

2. Repeat steps 2-4.

Is there anything wrong conceptually?

### Edit: Remarks on the acceptance probability

As Xi'an pointed out in the comment, if we start from $$x=\nu=14.5$$ and apply the forward transformation (to go to the other Gaussian at $$\mu = -4.5$$)

$$$$x' = \mu - (x + \mu)e^{-t},$$$$

we can see that we won't get close to $$\mu$$ unless $$e^{-t}$$ is very small. However the goal is not really to land at $$x' = \mu$$, but just at the region of the other Gaussian (so that the next random steps using $$x'=x+\epsilon$$ can do the job sampling that Gaussian).

We can plot the probability that this move will be accepted as a function of $$t$$, $$P(t) = e^{-t}p(x')/p(x)$$, starting from $$x := x_0 = 14.5 (=\nu)$$:

The probability is kind of lowish from here, around $$6\%$$ for $$t=2.2$$ (this is where I took the value for $$t$$ for the algorithm from, by the way). And we would land at $$x' = -5.6$$ for this value of $$t$$, which would do the job. For initial states $$x$$ not at the peak $$\nu$$ of the Gaussian, the probability of going to the other Gaussian is higher. Starting at $$x = 13$$:

I have checked that the overall acceptance of the forward transformation on a run of the algorithm is around $$13\%$$. The one of the inverse is around $$20\%$$.

### Results

With this, the algorithm is able to sample the two different Gaussians, but there is statistical discrepancy in the momenta of the distribution (around 10 $$\sigma$$s), and I'm fairly sure that it's not due to autocorrelation effects.

### The code

In case it helps, this is the code:

#function that returns p(x) for a given x, mu and nu
def p(x,mu,nu):
return (0.199471140200716*np.exp(-0.5*(x - mu)**2) + 0.199471140200716*np.exp(-0.5*(x - nu)**2))

mu = -4.5 #Left Gaussian
nu = 14.5 #Right Gaussian
t = 2.2 # time of the transformation
N = 1000000 #number of iterations
n = 10 #number of random steps per deterministic transformation step

x = [] #History of accepted states
x.append(-5.0) #Initial state

for i in range(0, N):

# n-1 out of n steps, perform the random transformation x' = x + np.random.uniform(-1,1)
if(i%n <= n-2):

x_prime = x[i] + np.random.uniform(-1,1)

alpha = np.random.uniform(0,1) #random number for the acceptance probability

# Acceptance step
if( alpha <  p(x_prime,mu,nu)/p(x[i],mu,nu) ):
x.append(x_prime) #if accepted, append x_prime
else:
x.append(x[i]) #if rejected, append x[i] again

# Once per n steps, perform the deterministic forward or inverse transformation
else:

u = np.random.uniform(0,1) #random number to decide forward or inverse transformation

#Perform forward transformation
if(u < 0.5):

if(x[i] < 0 ):
x_prime = nu + (-x[i]-nu)*np.exp(-t)
else:
x_prime = mu + (-x[i]-mu)*np.exp(-t)

alpha = np.random.uniform(0,1) #random number for the acceptance probability

# Acceptance step. I added a cutoff x_prime in [-15,25] to avoid underflow and overflow. Add the jacobian e^(-t)
if( x_prime > -15 and x_prime < 25 and np.log(alpha) < np.log(p(x_prime,mu,nu)/p(x[i],mu,nu))+(-t) ):
x.append(x_prime)
else:
x.append(x[i])

#Perform inverse transformation
else:

if(x[i] > 0 ):
x_prime = -nu + (-x[i]+nu)*np.exp(t)
else:
x_prime = -mu + (-x[i]+mu)*np.exp(t)

alpha = np.random.uniform(0,1) #random number for the acceptance probability

if( x_prime > -15 and x_prime < 25 and np.log(alpha) < np.log(p(x_prime,mu,nu)/p(x[i],mu,nu))+(t) ):
x.append(x_prime)
else:
x.append(x[i])

# Compute statistical momenta and errors

# <x>
x_mean = np.mean(x[100000::25])
x_mean_error = np.std(x[100000::25])/np.sqrt(len(x[100000::25]))

# <x^2>
x2_mean = np.mean(np.asarray(x[100000::25])**2)
x2_mean_error = np.std(np.asarray(x[100000::25])**2)/np.sqrt(len(x[100000::25]))

# <x^3>
x3_mean = np.mean(np.asarray(x[100000::25])**3)
x3_mean_error = np.std(np.asarray(x[100000::25])**3)/np.sqrt(len(x[100000::25]))

# Compute discrepancy with analytical results. mom1, mom2 and mom3 are the analytical results computed with SymPy.
mom1Hist.append(abs(mom1.evalf() - x_mean)/x_mean_error )
mom2Hist.append(abs(mom2.evalf() - x2_mean)/x2_mean_error )
mom3Hist.append(abs(mom3.evalf() - x3_mean)/x3_mean_error )



### Edit: Discrepancies

The statistical momenta I'm checking are $$,,$$ and $$<|x|>$$. The analytical results are:

And the averages I get with the algorithm are:

Which were obtained on a run with $$N = 20\times10^6$$ iterations and picking one state per 100, starting at the state 100000 (to avoid any effect due to the burn-in period), to avoid autocorrelation errors. So I did the average with 199000 uncorrelated states. I did some checks and the autocorrelation time is around 35, so picking 1 state per 100 should be safe to forget about autocorrelation (that's what I hoped at least). In Python code:

# Compute statistical momenta and errors

# <x>
x_mean = np.mean(x[100000::100])
x_mean_error = np.std(x[100000::100])/np.sqrt(len(x[100000::100]))

# <x^2>
x2_mean = np.mean(np.asarray(x[100000::100])**2)
x2_mean_error = np.std(np.asarray(x[100000::100])**2)/np.sqrt(len(x[100000::100]))

# <x^3>
x3_mean = np.mean(np.asarray(x[100000::100])**3)
x3_mean_error = np.std(np.asarray(x[100000::100])**3)/np.sqrt(len(x[100000::100]))

# <|x|>
x1abs_mean = np.mean(abs(np.asarray(x[100000::100]))**3)
x1abs_mean_error = np.std(abs(np.asarray(x[100000::100]))**3)/np.sqrt(len(x[100000::100]))


The discrepancies I get in $$\sigma$$s are, respectively,

Which I obtained by computing

$$$$\text{Discrepancy in }\sigma = \frac{\text{analytical}-\text{average}}{\text{error of average}}$$$$

And this discrepancy I get is the thing that worries me. I also checked with a package that accounts for the autocorrelation time error (which implements the automatic windowing procedure proposed by Ulli Wolff on his paper and accounts for the corresponding error), but I still get these discrepancies, so this makes me think that taking 1 state per 100 is safe.

Is it normal to obtain such discrepancies with uncorrelated samples? Or do you think I am wrong with my error analysis?

### Second Edit: Bad sampling with the deterministic transformation (probably due to sample-space overlap between forward and inverse transformation)

I have realized that when one plots just 1 state every $$n$$ states (and thus selecting the states that result from the proposal of the deterministic transformation), the discrepancy between the target distribution and the sampled one becomes apparent,

But when one plots every state, as I did in the section Results above, this discrepancy gets "corrected" because the step $$x' = x + \epsilon$$ samples the Gaussian correctly and overshadows the bad sampling, making it difficult to spot visually. Since in my Python code above I put $$n=10$$, just the 10% of the steps seem not to be sampling correctly.

This effect does not seem to be due to autocorrelation, since plotting 1 every $$6n$$ states (in my case, 1 every 60 states) the thing does not change much:

I read in this paper (Section 2.1.1, Informal discussion, without proof) that the regions of sample space covered by the forward and backward transformations have to be disjoint in order for detailed balance to hold. That is, if we had $$x' = x + \epsilon$$ with $$\epsilon > 0$$ as forward and $$x' = x - \epsilon$$ as inverse, then the forward would cover $$\mathcal X \in [x, \infty)$$, while the inverse would cover $$\mathcal X \in (-\infty, x]$$.

In the transformation of my code above, this condition is not met. So I tried to do a different transformation to see if this solved the problem:

Forward transformation: \begin{align} x' &= -3 x \qquad \text{if } x < 0 \\ x' &= - x/3 \qquad \text{if } x > 0 \end{align}

Inverse transformation: \begin{align} x' &= -x/3 \qquad \text{if } x > 0 \\ x' &= - 3x \qquad \text{if } x < 0 \end{align}

And the disagreement above vanishes:

### Second Edit: Disagreement also in Xi'an 's answer model

According to this non-overlap requirement, Xi'an's example would not be sampling correctly either. I checked with his code with a slight modification in the non-deterministic transformation:

\begin{align} x' = x + \epsilon \end{align}

with $$\epsilon \in [-0.2,0.2]$$ instead of $$[-1,1]$$, so that it is more difficult to jump to the other Gaussian if it's not with the deterministic transformation. With this, I found, taking 1 every $$n$$ states:

The effect is slight and completely overshadowed if one just plots every state, as in Xi'an's answer. However, if one increases the overlap region and changes the transformations to

$$x\longmapsto 2\times(-1)^{\mathbb I(x>0)}-x/2$$ $$x\longmapsto 2\times(-2)^{\mathbb I(x<0)}-2x$$

so to increase the overlap region, the disagreement is more apparent:

If one drops the first term of the transformations,

$$x\longmapsto -x/2$$ $$x\longmapsto -2x$$

hence meeting the no-overlap condition, the disagreement vanishes:

Code to reproduce, with comments where I changed something:

gnorm<-function(x)(dnorm(x,-3.5)+dnorm(x,4.5))/2.0 #just added the normalizing /2.0 factor
mh=rep(0,T<-1e7)
for(t in 2:T){
u=runif(5)
if(t %% 10 < 9){ #so that there are 9 steps of random and 1 deterministic
mhp=mh[t-1]+(2*u[2]-1)*0.2 # x' = x + [-0.2, 0.2]
mh[t]=ifelse(u[3]<gnorm(mhp)/gnorm(mh[t-1]),mhp,mh[t-1])
}else{
dx=1+(u[4]<.5)
mhp=ifelse(dx==1,
2*ifelse(mh[t-1]<0,1,-1)-mh[t-1]/2, # added a factor 2 so that disagreement is more apparent
2*2*ifelse(mh[t-1]<0,-1,1)-2*mh[t-1]) #same, to have the correct inverse
mh[t]=ifelse(u[5]<dx*gnorm(mhp)/gnorm(mh[t-1])/(3-dx),mhp,mh[t-1])
}}


and I plotted 1 every 10 states,

hist(mh[((seq_along(mh)) %% 10) == 9],breaks=350,freq = FALSE, col="black")

x <- seq(-10,10,0.001)


Is this the issue? Or I'm still mising something?

• Thanks for the answer. Regarding the Jacobian, I checked with a Gaussian as target distribution, $$p(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}$$ and with transformations \begin{align} x' = x e^{-t} \\ x' = x e^{t} \end{align} I recover the statistical momenta of the distribution nicely. If I put the inverse of the Jacobian, I don't recover them. As for the reversibility... that is what I'm afraid of, yes. But the thing is that it works with just one Gaussian... – David Albandea Jul 1 '20 at 11:24
• Maybe by knowing at which of the peaks I am (if x>0 or x<0) in order to do the transformation to the other Gaussian introduces something to the problem that I don't understand... I don't know – David Albandea Jul 1 '20 at 11:24
• I did an edit to clarify the way I obtain the discrepancies. I also wrote some remarks on the last thing you said about the transform from one Gaussian to another. – David Albandea Jul 2 '20 at 11:34
• Alright, I'm just happy to hear that there is, in principle, nothing wrong with the algorithm. Thanks for the help! – David Albandea Jul 2 '20 at 15:49
• The Monte Carlo error may well be much larger than expected. – Xi'an Jul 2 '20 at 15:53

For comparison, here is an R code targeting a $$\frac{1}{2} \mathcal N(7/2,1)+\frac{1}{2}\mathcal N(9/2,1)$$ mixture, using the deterministic switch $$x\longmapsto (-1)^{\mathbb I(x>0)}-x/2$$ and its reciprocal $$x\longmapsto (-2)^{\mathbb I(x<0)}-2x$$

gnorm<-function(x)dnorm(x,-3.5)+dnorm(x,4.5)
mh=rep(0,T<-1e7)
for(t in 2:T){
u=runif(5)
if(u[1]<.5){
mhp=mh[t-1]+2*u[2]-1
mh[t]=ifelse(u[3]<gnorm(mhp)/gnorm(mh[t-1]),mhp,mh[t-1])
}else{
dx=1+(u[4]<.5)
mhp=ifelse(dx==1,
ifelse(mh[t-1]<0,1,-1)-mh[t-1]/2,
2*ifelse(mh[t-1]<0,-1,1)-2*mh[t-1])
mh[t]=ifelse(u[5]<dx*gnorm(mhp)/gnorm(mh[t-1])/(3-dx),mhp,mh[t-1])
}}


which produces a very good fit of the target

and close approximations of the first moments

> mean(mh) # should be 0.5
[1] 0.500631
> mean(mh^2) # should be 17.25
[1] 17.26589

• Can this be regarded as a toy example of reversible jump MCMC? – David Albandea Jul 11 '20 at 16:42
• @DavidAlbandea: RJMCMC is usually reserved for MCMC working over several parameter spaces and different models, e.g., when "the number of things you do not know is one of the things you do no know" to quote Peter Green. – Xi'an Jul 11 '20 at 17:52