Negative Binomial Reparameterization

In a previous post, I showed that direct estimation of the p and r parameters in a negative binomial distribution could involve bad autocorrelation in the MCMC chains, and I suggested that there must be some standard reparameterization to solve the problem, and asked for a pointer. Dr. John Davey of the University of Edinburgh was good enough to point the way. The solution is indeed straight forward: In the direct estimation, the p and r parameters are given priors, while the mean m (and variance v) is derived from p and r. In the reparameterization, the m and r parameters are given priors, while the p (and variance v) is derived from m and r. Autocorrelation in the chains is greatly reduced. Here is an example.

The model specifications:

Parameterized by p and r:

model {
    for( i in 1 : N ) {
      y[i] ~ dnegbin( p , r )
    }
    p ~ dbeta(1.001,1.001)
    r ~ dgamma(0.01,0.01)
    m <- r*(1-p)/p
    v <- r*(1-p)/(p*p)
}

Parameterized by m and r:

model {
    for( i in 1 : N ) {
      y[i] ~ dnegbin( p , r )
    }
    p <- r/(r+m)
    r ~ dgamma(0.01,0.01)
    m ~ dgamma(0.01,0.01)
    v <- r*(1-p)/(p*p)
}

Results:

Parameterized by p and r:

Parameterized by m and r:

The posteriors are a bit different for the two parameterizations, because I used "generic" priors without any attempt to transform them to be equivalent.

Thanks to John Davey for this pointer!