4
$\begingroup$

I am trying to compute the log-likelihood in a Poisson regression in R. However, my computation produces negative infinity values for some observations. This is my code:

freq*exposure * log(lambda)- lambda - log(factorial(freq*exposure))

the $-\infty$ occurs in the factorial(freq*exposure) portion when the frequency and exposure are too large.

Does anyone know how to get around this? When I discretize my data into smaller bins the issue disappears, however, I get a worse fit.

UPDATE

The error occurs when I use these observations:

exposure = 2.994500e+03
frequency = 0.13224244
exposure*frequency = 396

UPDATE2

Based on @StupidWolf's suggestion I changed the log(factorial(freq*exposure)) in the log-likelihood to:

freq*exposure * log(lambda)- lambda - lfactorial(freq*exposure)

this solved the issue.

$\endgroup$
  • 1
    $\begingroup$ can you provide the values that give this error? one option is to use lfactorial() $\endgroup$ – StupidWolf Jul 3 at 12:25
  • $\begingroup$ @StupidWolf I have included the values. That sounds like a good idea. I should have thought of that $\endgroup$ – MarG Jul 3 at 12:46
  • 1
    $\begingroup$ Is there are a reason for not using dpois(freq*exposure, lambda, log=TRUE)? $\endgroup$ – Gordon Smyth Jul 3 at 12:49
  • $\begingroup$ @MarG see GordonSmyth's comments. I also wanted to ask why you don't want to use logLik() since it's a poisson regression $\endgroup$ – StupidWolf Jul 3 at 13:00
  • 1
    $\begingroup$ @StupidWolf see my previous answer. $\endgroup$ – MarG Jul 3 at 13:03
6
$\begingroup$

Welcome to the world of numerical overflow!

Using wolfram alpha I calculated $396! > 10^{859}$ i.e. a $1$ with $859$ zeros trailing behind it. The biggest number your computer can handle will be around $10^{308}$ (probably). Source. What's happened in your code is the computer has tried to compute factorial(exposure*frequency) and rounded this up to $+ \infty$, taking the negative log of this will return $-\infty$.

Using a command like lfactorial will circumvent this issue. A factorial is a product and a log of products is just a sum. lfactorial will compute the log factorial as

$$\texttt{lfactorial}(n) = \sum_{i=1}^n \log(i)$$ which is much more ''numerically stable'' than computing $$\texttt{log(factorial(n))} = \log \left\{ 1\times2\times\cdots\times n\right\}$$

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.