Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions #6845

amyoshino · 2023-08-04T23:07:16Z

What is this PR about?
Adds ICDF functions to Beta, Gamma, Chisquared and StudentT distributions.

This is a work in progress. The quantile functions depend on two main functions:

Inverse of the regularized incomplete beta function: betaincinv
Inverse to the regularized lower incomplete gamma function: gammaincinv

Those are implemented in SciPy but not in Pytensor yet. When importing the two functions from scipy.special and plugging them in the desired formula, I get a TypeError like this below:

"TypeError: ufunc 'betaincinv' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''"

It looks like the most straightforward way to make these icdf functions work is to implement them in Pytensor, in the same way the pt.betainc has been implemented

Any comments about if this is the right path? If so I am happy to open an issue in Pytensor as a feature request, and possibly helping on its development if nobody else wants to take it. But if there is any ideas on how to make it simpler than that, let me know.

References:

SciPy betaincinv: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.betaincinv.html
SciPy gammaincinv: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gammaincinv.html
ChiSquared: https://www.wolframalpha.com/input?i=chi+squared+distribution+quantile+function
Beta: https://www.wolframalpha.com/input?i=beta+distribution+quantile+function
Gamma: https://www.wolframalpha.com/input?i=gamma+distribution+quantile+function
StudentT: https://www.wolframalpha.com/input?i=student+distribution+quantile+function

Checklist

Explain important implementation details 👆
Make sure that the pre-commit linting/style checks pass.
Link relevant issues (preferably in nice commit messages)
Are the changes covered by tests and docstrings?
Fill out the short summary sections 👇

Major / Breaking Changes

...

New features

...

Bugfixes

...

Documentation

...

Maintenance

...

📚 Documentation preview 📚: https://pymc--6845.org.readthedocs.build/en/6845/

codecov · 2023-08-04T23:20:06Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 90.97%. Comparing base (a033261) to head (262d868).

❗ Current head 262d868 differs from pull request most recent head 00c1f34. Consider uploading reports for the commit 00c1f34 to get more accurate results

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #6845      +/-   ##
==========================================
- Coverage   92.29%   90.97%   -1.33%     
==========================================
  Files         100      100              
  Lines       16875    16887      +12     
==========================================
- Hits        15575    15363     -212     
- Misses       1300     1524     +224

Files	Coverage Δ
pymc/distributions/continuous.py	`97.81% <100.00%> (+0.02%)`	⬆️

... and 14 files with indirect coverage changes

amyoshino · 2024-01-07T19:01:37Z

@ricardoV94 now that pymc-devs/pytensor#502 is merged, we are able to complete this PR.

The tests are failing because betaincinv and gammaincinv cannot be imported from pytensor.tensor.math, I guess it is matter of time to pytensor to be updated. I will restart the checks when it happens. Locally all tests are passing.

Let me know what you think.

ricardoV94 · 2024-01-07T20:34:14Z

You can see which PyTensor version is being installed in the CI.

pymc/distributions/continuous.py

ricardoV94

Couple more comments. Looks great otherwise

pymc/distributions/continuous.py

tests/distributions/test_continuous.py

ricardoV94 · 2024-01-07T20:49:53Z

Yeah it's strange it's still installing 2.18.4, even though 2.18.5 was released already some days ago. We may need to be just a bit more patient?

amyoshino · 2024-01-07T22:01:23Z

Yeah it's strange it's still installing 2.18.4, even though 2.18.5 was released already some days ago. We may need to be just a bit more patient?

Right! That's what I thought, I saw the pytenor release 4 days ago and thought it was a good experience to try it now and see if the update immediately applies. No rush on it, it will work soon.

review-notebook-app · 2024-01-08T12:26:37Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

ricardoV94 · 2024-01-09T06:17:57Z

@maresb any idea why CI is not yet installing the last PyTensor? Was there an issue on the conda side?

maresb · 2024-01-09T14:11:39Z

There was a glitch in the deployment of 2.18.5 to anaconda.org.

I reran it, so hopefully it fixes things.

maresb · 2024-01-09T21:26:56Z

It uploaded successfully on anaconda.org, so I restarted the tests.

amyoshino · 2024-01-09T21:47:45Z

Interesting way to restart the tests! I searched for and the recommendation I found was to push an empty commit, but I think this is nicer.

ricardoV94 · 2024-01-09T22:07:49Z

You can also just ask a dev with write permissions. We can restart them from GitHub interface

ricardoV94 · 2024-01-29T17:49:53Z

Is the nan coming from pm.logcdf(pm.Gamma.dist(alpha=2, beta=5), np.inf)?

ricardoV94 · 2024-01-29T18:04:00Z

Sounds like PyTensor gammainc C code does not have a special case for np.inf:

https://github.com/pymc-devs/pytensor/blob/19e1a982e5a2e9ba2bd9d708133abf102c6fb835/pytensor/scalar/c_code/gamma.c#L217-L223

In contrast this scipy has a special case?
https://github.com/scipy/scipy/blob/4ecf5638c36c3a87b473056f5fffecfeddffd368/scipy/special/cephes/igam.c#L149-L151

amyoshino · 2024-02-04T21:49:39Z

pm.logcdf(pm.Gamma.dist(alpha=2, beta=5), np.inf)

Indeed, this code returns nan. I tried to fix it writing the special case for np.inf in the gamma.c function but it didn't work out yet. Maybe there are other places I need to change in order to make it work.

`DEVICE double GammaP (double n, double x)
{ /* --- regularized Gamma function P /
if ((n <= 0) || (x < 0)) return NPY_NAN; / check the function arguments /
if (x <= 0) return 0; / treat x = 0 as a special case */
if (x < n+1) return _series(n, x) *exp(n *log(x) -x -logGamma(n));
if (isinf(n)) {
if (isinf(x)) return NPY_NAN;
return 0;
}
if (isinf(x)) return 1;
return 1 -_cfrac(n, x) *exp(n log(x) -x -logGamma(n));
} / GammaP() */

/--------------------------------------------------------------------/

DEVICE double GammaQ (double n, double x)
{ /* --- regularized Gamma function Q /
if ((n <= 0) || (x < 0)) return NPY_NAN; / check the function arguments /
if (x <= 0) return 1; / treat x = 0 as a special case */
if (x < n+1) return 1 -_series(n, x) *exp(n *log(x) -x -logGamma(n));
if (isinf(n)) {
if (isinf(x)) return NPY_NAN;
return 1;
}
if (isinf(x)) return 0;
return _cfrac(n, x) *exp(n log(x) -x -logGamma(n));
} / GammaQ() */`

Even with this, the pm.logcdf(pm.Gamma.dist(alpha=2, beta=5), np.inf) still returns nan.

ricardoV94 · 2024-02-04T21:57:56Z

You need to delete pytensor cache or bump the Op c_cache_version number. Try pytensor-cache clear or pytensor cache-clear I never remember. If clear doesn't work use purge

amyoshino · 2024-02-05T22:44:04Z

Thanks for the pointers! You were on point in all suggestions, the logcdf problem and the solution to it.

After changing the gamma.c function and pytensor-cache clear and pytensor-cache purge (I needed that too), the Truncated icdf worked as expected and all the text_mixture.py tests pass.

I will open an issue and a PR in pytensor to change the gamma.c code.

ricardoV94 · 2024-02-05T23:07:40Z

That PR should probably do a similar fix for the complementary gammainc function (gammaincc?)

Edit: I see you were already onto it

ricardoV94 · 2024-03-09T19:57:28Z

@amyoshino PyTensor new release should be in conda soon. We will need to update the pin on PyTensor version here, as it was a "major" release (or in a separate PR if there are other things that need fix)

https://github.com/pymc-devs/pytensor/releases/tag/rel-2.19.0

amyoshino · 2024-03-16T01:18:06Z

@amyoshino PyTensor new release should be in conda soon. We will need to update the pin on PyTensor version here, as it was a "major" release (or in a separate PR if there are other things that need fix)

https://github.com/pymc-devs/pytensor/releases/tag/rel-2.19.0

That's awesome!! Looking forward to it to be updated!
Let me know if I need to do anything with respect to this PR as part of the process.

ricardoV94 · 2024-03-16T07:23:23Z

You need to bump the PyTensor dependency constraints on the PyMC side

amyoshino · 2024-03-16T14:11:54Z

Ah, I see! Alright, let me try that then! Thanks for the pointer.

amyoshino · 2024-03-19T22:30:39Z

Just pushed the bump dependency modifications in here. Seems like there is one test failing though, which wasn't failing when I ran the test suite locally.

ricardoV94 · 2024-03-20T18:36:19Z

@amyoshino I'm taking core of the Pytensor bump in #7203

ricardoV94 · 2024-03-20T18:37:57Z

Just pushed the bump dependency modifications in here. Seems like there is one test failing though, which wasn't failing when I ran the test suite locally.

May have been an unlucky seed, rerunning it

amyoshino · 2024-03-21T03:38:54Z

@amyoshino I'm taking core of the Pytensor bump in #7203

In much better hands now 😄

Just pushed the bump dependency modifications in here. Seems like there is one test failing though, which wasn't failing when I ran the test suite locally.

May have been an unlucky seed, rerunning it

It was indeed a bad luck, the only failing part was the mypy, I guess because of the very outdated repo. I rebased and left only the previous commits here, so it is failing now because PyTensor is still 2.18. I think the way to go then is to wait for your PR to be merged and try again, right?

ricardoV94 · 2024-03-21T17:57:10Z

I think the way to go then is to wait for your PR to be merged and try again, right?

Right. I was hoping it was an easier bump, but it also unpins a Python major version :D

amyoshino · 2024-03-22T18:04:02Z

Cool! It looks like the bump is merged! Looking forward to restart the checks when you think is appropriate.

ricardoV94 · 2024-03-22T19:07:45Z

Cool! It looks like the bump is merged! Looking forward to restart the checks when you think is appropriate.

You can rebase anytime to incorporate the bump

amyoshino · 2024-03-22T21:58:47Z

Cool! It looks like the bump is merged! Looking forward to restart the checks when you think is appropriate.

You can rebase anytime to incorporate the bump

Ah, true! I got so excited that the bump was merged and this potentially can be merged after so long that I forgot that 😅 My apologies haha

EDIT: Good news! I think it is ready to go now

…ions removing unnecessary parameterization test on studentT test removing unnecessary test on chisquared test

…check_icdf_parameters

ricardoV94

Amazing work @amyoshino !!

twiecki · 2024-03-25T13:19:58Z

Congrats, this was a tricky one!

amyoshino · 2024-03-25T22:04:09Z

That was fun! 😄
Thank you for the support through this journey!

amyoshino mentioned this pull request Aug 15, 2023

[Feature Request]: Add betaincinv and gammaincinv pymc-devs/pytensor#413

Closed

amyoshino mentioned this pull request Nov 16, 2023

Implement betaincinv and gammainc[c]inv functions pymc-devs/pytensor#502

Merged

6 tasks

amyoshino changed the title ~~WIP: Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions~~ Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions Jan 7, 2024

amyoshino marked this pull request as ready for review January 7, 2024 19:02

ricardoV94 requested changes Jan 7, 2024

View reviewed changes

pymc/distributions/continuous.py Outdated Show resolved Hide resolved

ricardoV94 added enhancements logprob labels Jan 7, 2024

ricardoV94 reviewed Jan 7, 2024

View reviewed changes

ricardoV94 added the bug label Jan 7, 2024

amyoshino force-pushed the icdf_gamma_beta branch from e971c9d to 11942b9 Compare January 8, 2024 12:26

amyoshino marked this pull request as draft January 8, 2024 12:27

amyoshino force-pushed the icdf_gamma_beta branch 5 times, most recently from c1d3d35 to babc479 Compare January 9, 2024 02:00

amyoshino marked this pull request as ready for review January 9, 2024 02:42

maresb closed this Jan 9, 2024

maresb reopened this Jan 9, 2024

This was referenced Feb 6, 2024

BUG: gammainc C code does not have a special case for np.inf pymc-devs/pytensor#633

Closed

Add inf special cases in C implementation of gamma related functions pymc-devs/pytensor#634

Merged

amyoshino force-pushed the icdf_gamma_beta branch from 9f4352e to 262d868 Compare March 21, 2024 03:15

amyoshino added 2 commits March 22, 2024 19:03

Add icdf functions for Beta, Gamma, Chisquared and Students distribut…

bdd9a65

…ions removing unnecessary parameterization test on studentT test removing unnecessary test on chisquared test

fixing check_parameters occurences in icdf functions, replacing with …

00c1f34

…check_icdf_parameters

amyoshino force-pushed the icdf_gamma_beta branch from 262d868 to 00c1f34 Compare March 22, 2024 22:04

ricardoV94 approved these changes Mar 25, 2024

View reviewed changes

ricardoV94 merged commit abe7bc9 into pymc-devs:main Mar 25, 2024
21 checks passed

ricardoV94 mentioned this pull request Mar 25, 2024

Add icdf functions for distributions #6612

Open

38 tasks

amyoshino deleted the icdf_gamma_beta branch March 28, 2024 17:03

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions #6845

Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions #6845

amyoshino commented Aug 4, 2023 •

edited

Loading

codecov bot commented Aug 4, 2023 •

edited

Loading

amyoshino commented Jan 7, 2024 •

edited

Loading

ricardoV94 commented Jan 7, 2024

ricardoV94 left a comment

ricardoV94 commented Jan 7, 2024

amyoshino commented Jan 7, 2024

review-notebook-app bot commented Jan 8, 2024

ricardoV94 commented Jan 9, 2024

maresb commented Jan 9, 2024 •

edited

Loading

maresb commented Jan 9, 2024

amyoshino commented Jan 9, 2024

ricardoV94 commented Jan 9, 2024 •

edited

Loading

ricardoV94 commented Jan 29, 2024 •

edited

Loading

ricardoV94 commented Jan 29, 2024

amyoshino commented Feb 4, 2024

ricardoV94 commented Feb 4, 2024

amyoshino commented Feb 5, 2024

ricardoV94 commented Feb 5, 2024 •

edited

Loading

ricardoV94 commented Mar 9, 2024

amyoshino commented Mar 16, 2024

ricardoV94 commented Mar 16, 2024

amyoshino commented Mar 16, 2024

amyoshino commented Mar 19, 2024

ricardoV94 commented Mar 20, 2024

ricardoV94 commented Mar 20, 2024

amyoshino commented Mar 21, 2024 •

edited

Loading

ricardoV94 commented Mar 21, 2024

amyoshino commented Mar 22, 2024

ricardoV94 commented Mar 22, 2024

amyoshino commented Mar 22, 2024 •

edited

Loading

ricardoV94 left a comment

twiecki commented Mar 25, 2024

amyoshino commented Mar 25, 2024

Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions #6845

Add icdf functions for Beta, Gamma, Chisquared and StudentT distributions #6845

Conversation

amyoshino commented Aug 4, 2023 • edited Loading

Major / Breaking Changes

New features

Bugfixes

Documentation

Maintenance

codecov bot commented Aug 4, 2023 • edited Loading

Codecov Report

amyoshino commented Jan 7, 2024 • edited Loading

ricardoV94 commented Jan 7, 2024

ricardoV94 left a comment

Choose a reason for hiding this comment

ricardoV94 commented Jan 7, 2024

amyoshino commented Jan 7, 2024

review-notebook-app bot commented Jan 8, 2024

ricardoV94 commented Jan 9, 2024

maresb commented Jan 9, 2024 • edited Loading

maresb commented Jan 9, 2024

amyoshino commented Jan 9, 2024

ricardoV94 commented Jan 9, 2024 • edited Loading

ricardoV94 commented Jan 29, 2024 • edited Loading

ricardoV94 commented Jan 29, 2024

amyoshino commented Feb 4, 2024

ricardoV94 commented Feb 4, 2024

amyoshino commented Feb 5, 2024

ricardoV94 commented Feb 5, 2024 • edited Loading

ricardoV94 commented Mar 9, 2024

amyoshino commented Mar 16, 2024

ricardoV94 commented Mar 16, 2024

amyoshino commented Mar 16, 2024

amyoshino commented Mar 19, 2024

ricardoV94 commented Mar 20, 2024

ricardoV94 commented Mar 20, 2024

amyoshino commented Mar 21, 2024 • edited Loading

ricardoV94 commented Mar 21, 2024

amyoshino commented Mar 22, 2024

ricardoV94 commented Mar 22, 2024

amyoshino commented Mar 22, 2024 • edited Loading

ricardoV94 left a comment

Choose a reason for hiding this comment

twiecki commented Mar 25, 2024

amyoshino commented Mar 25, 2024

amyoshino commented Aug 4, 2023 •

edited

Loading

codecov bot commented Aug 4, 2023 •

edited

Loading

amyoshino commented Jan 7, 2024 •

edited

Loading

maresb commented Jan 9, 2024 •

edited

Loading

ricardoV94 commented Jan 9, 2024 •

edited

Loading

ricardoV94 commented Jan 29, 2024 •

edited

Loading

ricardoV94 commented Feb 5, 2024 •

edited

Loading

amyoshino commented Mar 21, 2024 •

edited

Loading

amyoshino commented Mar 22, 2024 •

edited

Loading