Metropolis-Hastings simulation of independent geometric random variables

Question

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?

Answer 1

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