nbinom stuff in TMB

[Ben Bolker]

Hi, folks,

the (R)TMB versions of dnbinom() [dnbinom2 or dnbinom_robust] are not
as robust when approaching the Poisson limit as the base-R function.

With a bit of work I could write my own adjoint code incorporating
the special cases shown below, if I needed it ...

Or, should the built-in TMB code be enhanced?

Is this is an issue that affects anyone else?
Should this be a TMB issue instead?


## demonstration code

x <- 1; mu <- 1
var_minus_mu <- 10^seq(-20, 0, length = 101)
## var = mu*(1+mu/theta) -> var_minus_mu = mu^2/theta -> theta =
mu^2/var_minus_mu
thetavec <- mu^2/var_minus_mu
varvec <- mu + var_minus_mu ## will underflow ...

fixtmb <- function(x) Im(unclass(x))
v1 <- dnbinom(x = x, mu = mu, size = thetavec, log = TRUE)
v2 <- RTMB::dnbinom2(x = x, mu = mu, var = varvec, log = TRUE) |> fixtmb()
v3 <- RTMB::dnbinom_robust(x = x, log_mu = log(mu),
log_var_minus_mu = log(var_minus_mu),
log = TRUE) |> fixtmb()

par(las = 1, bty = "l")
matplot(thetavec, cbind(v1, v2, v3), type = "l", log = "x",
ylim = c(-2, 0), lwd = 2, col = c(1,2,4),
ylab = "log LL")
abline(h = dpois(1, 1, log = TRUE), lty = 2)



R's version of dnbinom (dnbinom_mu, here:

https://github.com/r-devel/r-svn/blob/b130c55d9f3f018a339432a1356fd5ef79cb357e/src/nmath/dnbinom.c#L70

uses a different expression when size >> x:

if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/
return R_D_exp(size * (size < mu ? log(size/(size+mu)) : log1p(-
mu/(size+mu))));
if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
/* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
double p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size)));
return R_D_exp(x * p - mu - lgamma1p(x) +
log1p(x*(x-1)/(2*size)));

[Ben Bolker]

Is there any interest in this? Should I post an issue/start working on it?

[Kasper Kristensen]

Ben, yes definitely! This has to fixed.

The apparently not so robust dnbinom is here:

https://github.com/kaskr/adcomp/blob/943a9ecf02838381c1d7094441b2880490cae53f/TMB/inst/include/tiny_ad/robust/distributions.hpp#L47

The inaccuracy is from calculating

lgamma(x + n) - lgamma(n)

for large n.

It should be reasonably easy to fix given that this code is used with forward mode AD, so branching is not an issue.

[Ben Bolker]

lgamma(x + n) - lgamma(n) - lgamma(x + 1.) is lchoose(x+n-1, n-1)  (we could also probably get there from lbeta).  Don't know which of these pieces might be available in whatever C libraries we have access to (do we have Rmath?)   Do we need digamma for the gradients?