Tweedie loglikelihood returning inf for some values and impacting constrained optimization. Related to use of scipy's wrigh_bessel function in computing likelihood.

[diegodebrito]

Tweedie loglikelihood returning inf for some values and impacting constrained optimization. Related to use of scipy's wrigh_bessel function in computing likelihood.

[lorentzenchr]

The first (non-zero) value is x = 1.452e+5 and alpha=-1. Now, wright_bessel(1, 0, 1e5) = 2.356e+275 and np.finfo(np.float64).max = 1.798e+308. So the simple answer is that wright_bessel(1, 0, 1.4e5) is simply larger than the max floating point number.

I think the tweedie package treats it differently and works with the logarithm of Wright's Bessel function, thus avoiding the overflow.

[diegodebrito]

Thanks for the answer and for opening the Scipy issue, @lorentzenchr

In your opinion, would it be better to wait/push for the Scipy update and then modify the wright_bessel here? Or would it be a good idea to remove wright_bessel and try to add the more detailed calculations here? I noticed that glum and tweedie packages have the latter. I'm not familiar with Scipy's release cycle, so not sure what the timeline for that feature would be.

Also interested on your input here, @josef-pkt

[lorentzenchr]

That's a question for @josef-pkt to answer and it also depends on you expected timeline. I will only pursue the "proper" solution with scipy. Let's see how they respond.

[josef-pkt]

(I'm currently stuck in a very different neighborhood of statistics and have not looked at the details for this.)

If there is a workaround that can be added in statsmodels, then it would be good to have it because it would also work with older scipy versions, e.g. we copied some code from scipy to statsmodels.compat to make parts available faster.

However, if there is no good replacement for (compiled) scipy.special functions, then we just have to wait for the fix/enhancement in scipy and the release cycle.

[diegodebrito]

Great point, @josef-pkt. An alternative implementation wouldn't limit users to the eventual version of SciPy with the fix. I'm interest in working on this. We would love to have everything available in Statsmodels and not have to rely on the Tweedie package.

A question for you both regarding licenses @josef-pkt @lorentzenchr Can I read code available in other packages to help with this? Or should I rely on papers only? I think Tweedie and GLUM licenses are pretty permissive, but I don't want to start looking at the source code if I'm not allowed to.

[josef-pkt]

GLUM https://github.com/Quantco/glum is also BSD-3 licensed, so there is no restriction looking at their code, or even copying parts.
@thequackdaddy was the original contributor of tweedie in statsmodels. His package might not have much additional to our code, but I haven't looked at it.

[josef-pkt]

related:
tweedie family is missing get_distribution. We currently do not have a distribution class for it, AFAICS

[lorentzenchr]

With scipy/scipy#20646 merged, as of (future) scipy 1.14, statsmodels can replace log(wright_bessel) with log_wright_bessel. This should solve this issue.