Gibbs algorithm using negative binomial produces NAs

Question

I have the following full conditionals distributions: $$ X_2|X_1=x_1\sim Bin(x_1,p)\\ X_1|X_2=x_2\sim NegBin(x_2,p) $$ So I'm using the following code to generate a sample from each one:

Gibbs2 <- function(x1,x2,p,niter){
lambda = matrix(0, nrow=niter)
betax= matrix(0, nrow=niter)
lambda[1]= x1
betax[1]= x2
for (i in 2:niter) {
  lambda[i]= rbinom(1, betax[i-1],p)
  betax[i]= rnbinom(1, lambda[i],p)
}
return(theta = list(lambda=lambda,betax=betax))

Using the function returns the following:

    > cadenas2<-Gibbs2(20,40,0.5,100)
There were 50 or more warnings (use warnings() to see the first 50)

But when I change the function to both distributions being binomial, so using rbinom for both cases, then it works. Does anyone know what's wrong with my code? Why does rnbinom produce such error?

Answer 1

In short, the problem behind the question is statistically meaningless if probabilistically interesting.

When using the definition of the Negative Binomial random variate as the number of failures, the two conditional distributions are actually incompatible with one another, that is, there is no joint distribution on $(X_1,X_2)$ producing these conditionals. If instead one uses the definition of the Negative Binomial random variate as the number of attempts, then they are compatible: \begin{align*} f(x_1|x_2) &= {x_1-1 \choose x_2-1} p^{x_2} (1-p)^{x_1-x_2}\Bbb 1_{x_1\ge x_2}\\ f(x_2|x_1) &= {x_1 \choose x_2} p^{x_2} (1-p)^{x_1-x_2}\Bbb 1_{x_1\ge x_2}\\ \end{align*} implies that $$ \dfrac{f(x_1|x_2)}{f(x_2|x_1)} = \dfrac{{x_1-1 \choose x_2-1} p^{x_2} (1-p)^{x_1-x_2}}{x_1 \choose x_2} p^{x_2} (1-p)^{x_1-x_2}=\dfrac{x_2}{x_1} $$ factorises as a product of a function of $x_1$ and of a function of $x_2$. In this case, the only issue to define properly $f(x_1|x_2)$ when $x_2=0$ since this is not the case for rnbinom(1,0,p). I suggest using a Dirac mass at x_1=1.

if(lambda[i]>0) betax[i]=lambda[i]+rnbinom(1, lambda[i],p)

In which case the corrected R code produces a sequence

However, given the above ratio I strongly suspect the joint measure on $(X_1,X_2)$ is not proper i.e. sums up to $\infty$. The Markov chain is thus either null recurrent or transient.

Note that if both conditionals are Binomial distributions, $\mathcal B(x_1,p)$ and $\mathcal B(x_2,p)$, respectively, then

1. the conditionals are also incompatible with one another
2. the resulting Markov chain always converges monotonically (i.e., decreases) to zero, $(0,0)$.

Answer 2

The lambda and betax can become zero during the sampling, in which case you have a problem because the size in both rbinom() and rnbinom() must be strictly positive.