first draft of Binomial.icdf issue #6612 - PyData 2024 Hackathon

[niknow]

Description

Related Issue

• Closes #
• Related to Add icdf functions for distributions #6612

Checklist

• Checked that the pre-commit linting/style checks pass
• Included tests that prove the fix is effective or that the new feature works
• Added necessary documentation (docstrings and/or example notebooks)
• If you are a pro: each commit corresponds to a relevant logical change

Type of change

• New feature / enhancement
• Bug fix
• Documentation
• Maintenance
• Other (please specify):

---

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

[welcome]

]
💖 Thanks for opening this pull request! 💖 The PyMC community really appreciates your time and effort to contribute to the project. Please make sure you have read our Contributing Guidelines and filled in our pull request template to the best of your ability.

[ricardoV94]

This is clever but widely inefficient for large n?

[ricardoV94]

Also would probably need to pass axis=-1 to argmax

[niknow]

You are totally right, it is slow for large n as it is in O(n). I did actually flag this point at the hackathon and would agree that for that reason this is an ad-hoc implementation, but was encouraged to submit it anyway. ;) If you feel it is too slow for production I'm totally fine if you reject this PR.

I would note though that for very large n, one is usually better off using a normal approximation instead anyway.

And on the upside: This trick can be used to quickly get an implementation for the icdf of any discrete distribution, so one idea we had at the event was to use that as a first draft and make it quicker later if needed. One might be able to get the complexity down to O(log(n)) one replaces this with some sort of bisection.

Having a fast icdf just for Binomial might require quite a bit more effort, e.g. some numerical root finder of the CDF (which for Binomial is the regularized incomplete beta function). So I think the options are:

1. take this idea as a draft and reject for production use (maybe keep in mind using for unit tests)
2. take this idea and accept that it is slow for now
3. try to tweak to O(log(n))
4. try to lift and shift a more tailored approach based on beta function (will probably require a bit higher effort)

Of course open to ideas if you see another option? ;-)

[ricardoV94]

Thanks for the explanation!

I think we should try 3/4, and we can use this approach to obtain and independent reference for testing?

Perhaps the binary search wouldn't be too hard to implement with a ScalarLoop (so it is trivial to vectorize as an Elemwise operation)?

https://github.com/pymc-devs/pytensor/blob/main/pytensor/scalar/loop.py

The requirement is that the CDF expression be composed only of Scalar (or Elemwise in the batched version) operations.

The code may be a bit daunting but here are some cases where we use it: https://github.com/pymc-devs/pytensor/blob/efa845a3484915e4e15a928918fa97d081886d50/pytensor/scalar/math.py#L870

Or tests that may be more readable: https://github.com/pymc-devs/pytensor/blob/main/tests/scalar/test_loop.py

[ricardoV94]

A quick solution may be to wrap scipy binom_ppf: https://github.com/scipy/scipy/blob/0c4cfab8b27d85785d28ec0e4b8d1ce536276c0e/scipy/stats/_new_distributions.py#L500-L501

Which is calling boost under the hood. It uses an iterative solver based on a good first guess.