# Metropolis-Hastings simulation of independent geometric random variables

Consider the following Metropolis-Hastings scheme to sample independent geometric random variables $X = (X_1, \dots, X_N)$, where each $X_j$ has pmf $\mathbb{P}(X_j = x) = p(1-p)^x$ for $x \geq 0$. At iteration $t$ we either propose (with probability $1/2$) an adding move, where for $j = 1,\dots, N$, $$x_{j,t}^\prime = \begin{cases} x_{j, t-1} + 1 & \text{with probability } 1/2,\\ x_{j, t-1}& \text{otherwise}. \end{cases}$$ or a removing move (with probability $1/2$), where for $j = 1,\dots, N$, $$x_{j,t}^\prime = \begin{cases} x_{j, t-1} - 1 & \text{with probability } 1/2 \text{ if x_{j, t-1} > 0},\\ x_{j, t-1}& \text{otherwise}. \end{cases}$$ The Metropolis-Hastings acceptance probability for an adding move is below, where $\delta(j)$ is a binary function that denotes whether we proposed an increase for RV $X_j$. $$\min \left\{1, \frac{1/2}{1/2} \prod_{j = 1}^N (1-p)^{\delta(j)} \frac{(1/2)^{1 - \delta(j)}(1/2)^{\delta(j)}}{(1/2)^{\delta(j)}(1/2)^{1 - \delta(j)}} \right\} = \min\left\{1, \prod_{j=1}^N (1-p)^{\delta(j)}\right\},$$ where the $(1-p)^{\delta(j)}$ term is the ratio of the target distribution, the first fraction is the ratio of the probabilities of proposing an adding move and a removing move, and the second fraction is the ratio of the proposal probabilities given a proposing or removing move. In the case where $\delta(j) = 0$ for all $j$ then the acceptance probability formula is different because $x_t = x_{t-1}$ and this could be the result of an adding or a removing move, but the acceptance probability still equals 1.

The probability of accepting a removing move is $$\min \left\{1, \frac{1/2}{1/2} \prod_{j = 1}^N \left[ (1-p)^{-\gamma(j)} \frac{(1/2)^{1 - \gamma(j)}(1/2)^{\gamma(j)}}{(1/2)^{\gamma(j)}(1/2)^{1 - \gamma(j)}}\right]^{\mathbb{I}\{x_{j,t-1} > 0\}} \right\} = \min\left\{1, \prod_{j=1}^N (1-p)^{-\gamma(j) \mathbb{I}\{x_{j,t-1} > 0\}}\right\},$$ where $\gamma(j)$ denotes whether we proposed a decrease for $X_j$.

If I run this with $N=1, p=1/2$ then I produce a sample from a geometric($1/2$) distribution, but with $N=2, p = 1/2$ I produce two independent samples from a geometric($5/8$) distribution. Any idea why?

• Could you explain why you introduce a possibility of keeping the same value in the proposal? This reduces efficiency for no clear reason. I also fail to understand why you define two Metropolis-Hastings moves rather than a single one that would propose moving by $\pm 1$ or remaining at the same value. – Xi'an Jan 6 '16 at 13:10
• I accept that this is not the most efficient method for performing this sampling, I posted the question because I'm really curious where the problem is: why does it work for j = 1 but not for higher j? – Alex Jan 6 '16 at 13:29
• I mean N not j in the previous comment. – Alex Jan 6 '16 at 14:16

If you consider each Metropolis-Hastings scheme separately with a separate acceptance probability, i.e., distinguish adding from removing as two separate Metropolis-Hastings schemes, the chain associated with each scheme is not irreducible and hence the validity of the overall algorithm is not established. Actually, each proposal is invalid in terms of the intended target since it either never decreases or never increases.

To validate the scheme properly, you need to use the mixture of both proposals in the ratio. In this case, the proposal is $$Q(\mathbf{x},\mathbf{x}')=3^{n_{00}}\big/ 2^{n_e}4^{n_d+n_{00}}$$ where $n_{00}$ is the number of pairs $(x_i,x_i')$ that are equal to $(0,0)$, $n_e$ the number of pairs $(x_i,x_i')$ that are identical but not to $(0,0)$, and $n_d$ the number of pairs $(x_i,x_i')$ that differ. Unless I am confused this leads to a Metropolis-Hastings acceptance probability equal to the ratio of the target densities, i.e. $$\min\{1,(1-p)^{n_+-n_-}\}$$where $n_+$ is the number of moves up and $n_-$ the number of moves down. This is because $$Q(\mathbf{x},\mathbf{x}')=Q(\mathbf{x}',\mathbf{x})$$

Here is a naïve R code running the above algorithm:

p=.7
N=10
T=1e4
mc=matrix(sample(0:10,N*T,rep=TRUE),ncol=N,nrow=T)
for (t in 2:T){
for (i in 1:N){
if (mc[t-1,i]==0){ mc[t,i]=sample(c(0,1),1,prob=c(3,1))}else{
mc[t,i]=mc[t-1,i]+sample(c(-1,0,1),1,prob=c(1,2,1))}}
prob=(1-p)^{sum(mc[t,]-mc[t-1,])}
if (runif(1)>prob) mc[t,]=mc[t-1,]}


and producing the correct expectation:

> apply(mc+1,2,mean)
[1] 1.4405 1.5032 1.4090 1.4497 1.4254 1.5170 1.4529 1.5191 1.4581 1.4544

• Thanks, I accept that the scheme you propose is correct, but I don't see why the chain I discussed is irreducible. If the sampler is in state $(x_1, x_2)$ then state $(x_1 + k_1, x_2 + k_2)$ can be reached by $|k_1|$ increases/decreases to just $x_1$ (increase/decrease depending on whether $k_1$ is positive or negative) and $|k_2|$ increases/decreases to just $x_2$. – Alex Jan 7 '16 at 20:40
• The mixture scheme is irreducible. Both Metropolis schemes are by themselves transient since the proposal Q either always goes up or goes down. With the given Q, Q(x,x+1)=1/2 implies that Q(x,x-1)=0, hence the Metropolis scheme is not stationary against the geometric target. – Xi'an Jan 7 '16 at 20:45
• To clarify, the sampler chooses whether to propose an adding move or a removing move at each iteration. So if $N = 2, x = (1, 1)$ then the possibilities for $x^\prime$ are $(0, 0), (0, 1), (1, 0), (1, 1), (2, 1), (1, 2), (2, 2)$, where the first three are the result of removing moves, $(1, 1)$ could be a result of either adding or removing, and the last three are the result of adding moves. With the exception of $(1, 1)$ each transition has probability $1/8$ in this case. – Alex Jan 8 '16 at 8:26
• @Alex: let me try one last time: if you consider the iteration as proposing adding, removing or unchanging, there is a single Metropolis-Hastings acceptance probability and the corresponding Metropolis-Hastings algorithm is perfectly valid. This is the one I coded. If instead you separate the adding from the removing moves, meaning all components in your vector are subjected to the same type of move, then the corresponding Metropolis-Hastings algorithm is not valid. – Xi'an Jan 9 '16 at 13:27