# Gibbs algorithm using negative binomial produces NAs

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= x1
betax= 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?

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)$$.
• Could you explain why does that happen with both conditionals being binomial? Do you have any sources where I can read deeper into that? Jul 14, 2019 at 21:00
• @andrespm: it is straightforward: if $X_1^t\sim\cal B(X_2^{t-1},p)$ and $X_2^t\sim\cal B(X_1^{t},p)$,$$X_2^t\le X_1^{t}\le X_2^{t-1}$$and$$X_1^{t}\le X_2^{t-1}\le X_1^{t-1}$$ Jul 15, 2019 at 4:16

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.

• Why does it work when both are generated from a binomial distribution though? Jul 14, 2019 at 4:19
• Because the negative binomial distribution gives you the number of failures before the first success in a binomial experiment. And the binomial distribution gives you the number of successes in $N$ trials. Jul 14, 2019 at 6:59
• @DimitrisRizopoulos: as indicated in my answer, the simulation scheme does not "work" either with binomial conditionals as it converges to zero. It simply does not produce NAs, because rbinom(1,0,p) returns a zero. Jul 23, 2019 at 5:55