Preventing numerical 0s in hurdle model

Modeling
1

Here’s an example version of the code I’m running for a hurdle Negative-binomial model:

model {
  for (n in 1:N){
      if(Y[n]==0){
        target += -log1p_exp(-logittheta) ;
      }else{
        print("n = ",n) ;
        print("Prob > 0 is ",1-(1-mucond[n])^psi) ;
        print("psi is ",psi) ;
        print("mucond is ",mucond[n]) ;
        print("lccdf ",neg_binomial_lccdf(0|psi,(1-mucond[n])/mucond[n])) ;
        
        target += -log1p_exp(logittheta) 
        + neg_binomial_lpmf(Y[n]|psi,(1-mucond[n])/mucond[n]) 
        - neg_binomial_lccdf(0|psi,(1-mucond[n])/mucond[n]) ;
      }
    }
}

The issue I am running into is the the neg_binomial_lccdf is returning negative infinity because mucond is very very close to zero, so I am getting some rounding errors. I know this because I have checked each quantity used in computing neg_binomial_lccdf. Here’st the error:

Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

Any ideas on how to handle this?

Some more example output:

n = 99
Prob > 0 is 0
psi is 1
mucond is 3.46945e-18
lccdf -inf

1 Like
2

Hey there!

I guess the key is to find suitable priors (and perhaps initial values). If the priors are very diffuse on the log scale, the chance they will be super close to zero is pretty high. Also, the parameterization of your predictors matters!

Something that is very hacky is to add a small number to mu… however, this will obviously bias your results…

Btw, following the Stan manual, I coded a hurdle Neg-Binomial hurdle like this:

data {
  int<lower=0> N;
  int Y[N];
  vector[N] x;
}
parameters {
  real alpha;
  real beta;
  real logit_theta;
  real<lower=0> phi;
}
transformed parameters{
  vector[N] mu = exp(alpha + beta*x) + 1e-10;
}
model {
  for (n in 1:N){
      if(Y[n]==0){
        target += log_inv_logit(logit_theta);
      } else {
        target += log1m_inv_logit(logit_theta) + neg_binomial_2_lpmf(Y[n] | mu[n], phi) - neg_binomial_2_lccdf(0 | mu[n], phi);
      }
    }
  alpha ~ normal(-1.5, 2.5);
  beta ~ normal(0, 1.5);
  logit_theta ~ normal(-0.5, 0.5);
  phi ~ exponential(1);
}

…which kind of worked for the example data I used.

The goal would probably be to have a more stable neg_binomial_lccdf function in Stan… :/

@martinmodrak might have a better idea (he knows about numerical difficulties of the NegBinomial…)

Cheers,
Max

1 Like
3

Agree that prior that would disallow mu_cond to get too close to zero and hence (1-mucond[n])/mucond[n]) not extremely huge and hence mean larger than zero is probably reasonable (I guess you can a priori say that the “over-the-hurdle” part of the model is not negligible).

Not sure to what extent could this be blamed on numerics: maybe negative infinity is actually the correct answer here? (but probably not). When I struggled with the neg binomial I didn’t need (and hence didn’t touch) the lcdf and lccdf functions. It would help if you could check if some other trusted implementation (e.g. R or Python or whichever language you use) arrives at a finite result - if so, this is likely a bug and I would kindly ask you to report it at https://github.com/stan-dev/math/issues

EDIT: One more thing: If I am not mistaken, then neg_binomial_lccdf(0|a, b) == log(1 - exp(neg_binomial_lpmf(0| a, b)), right? But the formula at https://mc-stan.org/docs/2_23/functions-reference/negative-binomial-distribution.html actually simplifies quite neatly when n = 0 to (please double check my computations) so we get:

neg_binomial_lccdf(0|a, b) == log1m( (b / (b + 1)) ^ a) == log1m_exp(a * (log(b) - log1p(b))

At least one of the latter two formulae could be numerically much more stable.

EDIT 2: Noticed that you already calculate the probability > 0 (and it seems to agree with my computations, so I am glad I didn’t mess up somehwere). The problem probably is elsewhere: if mucond <- 3.46945e-18 then 1 - mucond is numerically exactly 1 and we’ve already lost so much precision.
We then get b = 288230122930147456 and b / (b + 1) is numerically exactly 1 and log(b) - log1p(b) is numerically exactly 0 and this propagets to -Inf as result.

But if we calculate the log of the probability that Y > 0 as log1m_exp(psi * log1m(mucond)) we get -40.202535593937661 which looks about right!

We also can get this result from b if we compute neg_binomial_lccdf(0|a, b) == log1m_exp(-a * log1p(1/b)). So this means that there is a numerical bug in neg_binomial_lccdf (issue here: https://github.com/stan-dev/math/issues/2031).

Best of luck with your model!

2 Likes
4

Thanks, Martin. Using:

log1m_exp(psi * log1m(mucond))

Resolved my issue.

1 Like