# Negative infinity produced when computing log-likelihood in Poisson Regression R

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.

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

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\}$$