scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
13.07k stars 5.19k forks source link

ENH: gammainc to support not regularized lower incomplete gamma function. #18619

Open JoeJoe1313 opened 1 year ago

JoeJoe1313 commented 1 year ago

Is your feature request related to a problem? Please describe.

No response

Describe the solution you'd like.

It would be nice to add an argument regularized: bool and have the possibility to choose between regularized calculation and not regularized.

Describe alternatives you've considered.

No response

Additional context (e.g. screenshots, GIFs)

No response

entaylor commented 5 months ago

+1 for this; commenting for updates.

While i have been googling around for a solution to this, it seems it is also possible to obtain values for a < 0 through recurrence relations. (See eg here). If/when someone is coding up the un-regularised solution, it might also be worthwhile implementing the iterative solution via recurrence for such cases?

My understanding is that sympy has this capability, which might provide a convenient check/benchmark.

lucascolley commented 5 months ago

@mdhaber you wrote a little about gamma and regularized in https://github.com/data-apis/array-api/issues/725. Any thoughts here?

mdhaber commented 5 months ago

+1. I imagine this is on the radar of @steppi and @izaid for the new library.

izaid commented 5 months ago

It's on the radar. https://github.com/scipy/scipy/pull/20539 adds normalised Legendre polynomials via a keyword argument, we'd use some of that to add a keyword argument for regularised gamma functions

entaylor commented 4 months ago

Thanks @mdhaber and @izaid for your efforts and your quick response on this! I have continued to try to find potential workarounds, but have been stymied at every turn! Adding this functionality will unlock my problem.

Not to be a pest, but just to try to manage my own projects, is it possible to know a rough timeline for this?

mdhaber commented 4 months ago

I'll let the others reply about timeline, but perhaps I can help in the meantime. Can you describe a bit more of what you need? What workarounds are not working? For example, why is numerical integration not working (if it's something you need immediately)? If it's accuracy, why not use mpmath? If speed is the main concern, do you also need high accuracy? Never mind, I see that the integrand at x=0 is infinite for negative a, and numerical integration doesn't give a useful result. I'm going to guess that someone else will reply shortly.

entaylor commented 4 months ago

Thanks, @mdhaber . Yes, the first requirement is to get the unregularised incomplete gamma function, to allow us to consider a < 0. As i say, this is possible with e.g. sympy or mpmath. But the second requirement is to have something that we can use with pymc or numpyro, which seem (unless i'm mistaken?) to inherit the scipy gamma function definitions via jax.

I think our alternative is to try to code our own implementation of the unregularised incomplete gamma function using jax, but ... i'm a bit daunted by that prospect!