Add icdf for Wald distribution

[TalenMud]

Description

Adds icdf (inverse CDF / quantile function) for the Wald (Inverse Gaussian) distribution, as part of the broader effort to implement ICDF methods across all distributions.

Currently a draft skeleton with parameter checks is in place. ICDF formula and tests to follow.

Related Issue

• 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):

[read-the-docs-community]

Documentation build overview

📚 pymc | 🛠️ Build #33913954 | 📁 Comparing c86a83a against latest (c870b5c)

  🔍 Preview build  

16 files changed · ± 16 modified

± Modified

• glossary.html
• api/math.html
• contributing/jupyter_style.html
• api/generated/pymc.backends.zarr.ZarrTrace.html
• api/generated/pymc.sample_posterior_predictive.html
• api/generated/pymc.sample_prior_predictive.html
• api/generated/pymc.sampling.jax.sample_blackjax_nuts.html
• api/generated/pymc.sampling.jax.sample_numpyro_nuts.html
• api/generated/pymc.to_inference_data.html
• learn/core_notebooks/dimensionality.html
• and 6 more...

[TalenMud]

Hi, I've been researching the Wald ICDF implementation and found that unlike many other distributions it has no closed-form quantile function. The most reliable approach appears to be a Newton's method iteration started from the mode of the distribution, as described in Giner & Smyth (2016) "statmod: Probability Calculations for the Inverse Gaussian Distribution" https://arxiv.org/abs/1603.06687. This guarantees monotonic convergence for all parameter values.

Before proceeding I wanted to check, is the preferred approach to implement this as a proper Newton iteration using pytensor.scan, or would a scipy wrapper via SplineWrapper be more appropriate for PyMC's codebase? I'm keen to go with the full Newton implementation as it gives full machine accuracy and works for all parameter values, but happy to go the scipy wrapper route if that's what the maintainers prefer.

[ricardoV94]

@TalenMud sorry for the delay and thanks for the initiative. Yes we've ran out of closed form solutions. You can check how #8339 does a binary search for discrete variables. For continuous you'll likely need something like you said.

If you want to look at an implementation boost has one here: https://github.com/boostorg/math/blob/8ee12a5355935cbaac5d5338372d0d0e3311b473/include/boost/math/distributions/inverse_gaussian.hpp#L343-L402

They are usually not how you'd write in PyMC/PyTensor but can be helpful as a reference. Feel free to take a stab. Note PyTensor has a minimize/root Ops that we can use if that's simpler than rolling a custom Scan

[TalenMud]

Thanks for the pointers! I'll take a look at the Boost implementation and #8339 for the pattern. I'll also check out the PyTensor root Ops as the approach, it sounds better than rolling a custom Scan.

[TalenMud]

@ricardoV94 ready for review! Newton's method implementation starting from the mode per Giner & Smyth (2016). Validates against scipy.stats.invgauss.ppf to 7 significant figures. Happy to add tests once the approach is confirmed. Also I decided to follow #8339 and use a pytensor scan.

[ricardoV94]

@ricardoV94 ready for review! Newton's method implementation starting from the mode per Giner & Smyth (2016). Validates against scipy.stats.invgauss.ppf to 7 significant figures. Happy to add tests once the approach is confirmed. Also I decided to follow #8339 and use a pytensor scan.

I'm a bit busy but it's on my stack!

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.77%. Comparing base (c870b5c) to head (c382f68).
⚠️ Report is 37 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #8357      +/-   ##
==========================================
+ Coverage   91.73%   91.77%   +0.03%     
==========================================
  Files         128      128              
  Lines       20697    21124     +427     
==========================================
+ Hits        18987    19386     +399     
- Misses       1710     1738      +28     

Files with missing lines	Coverage Δ
pymc/distributions/continuous.py	97.07% <100.00%> (-1.20%)	⬇️

... and 15 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ricardoV94]

In the meantime we'll need a test either way

[TalenMud]

Working on the Tests now!

[TalenMud]

@ricardoV94 tests added! A few notes on the implementation decisions:

1. Test scope: the test covers lam=1 (default) only. During testing I discovered that Wald.logcdf gives incorrect results for lam != 1 compared to scipy — this is a pre-existing issue unrelated to this PR. Happy to open a separate issue for that if useful.

2. Precision: relaxed to decimal=4 rather than the default 6. The relative error is excellent (~3.4e-08) but for large parameter values (e.g. mu=100) the absolute difference exceeds 6 decimal places. This is expected behaviour for an iterative numerical method.

3. Taylor approximation: added pt.switch for large kappa values to avoid catastrophic cancellation in the mode calculation mu * (sqrt(1 + kappa^2) - kappa) when kappa > 1e3.

4. Scan deprecation: fixed return_updates=False to conform to the new PyTensor API.

[ricardoV94]

Wald.logcdf gives incorrect results for lam != 1 compared to scipy — this is a pre-existing issue unrelated to this PR. Happy to open a separate issue for that if useful.

Would be great to fix that

[TalenMud]

Would be great to fix that

Would you like me to open another issue for it?
Since I think it'd be out of scope for this PR, what do you think?

[ricardoV94]

Would be great to fix that

Would you like me to open another issue for it? Since I think it'd be out of scope for this PR, what do you think?

Sure