Probability Distribution -- Gamma Distribution

[Calluumm]

Adds a Gamma Distribution to Stats
This is a modernisation/continuance of #278 all credit to the original author.

Comparatively to #278;

• Tests have been standardised to current stdlib standards
• Documentation lightly amended
• Filepaths adjusted for current structure

ctest distribution_gamma passed fine
the function itself was previously tested against SciPy so I did again it's still accurate, I further tested against FSML's function and again it was identical.
specialfunctions did have to be linked into stats for the lower incomplete gamma function, if this isn't standard it could be duplicated into stats but I'm not sure what the practice is there for this lib

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 68.55%. Comparing base (de3e59f) to head (29c1580).
⚠️ Report is 21 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1099   +/-   ##
=======================================
  Coverage   68.55%   68.55%           
=======================================
  Files         396      396           
  Lines       12746    12746           
  Branches     1376     1376           
=======================================
  Hits         8738     8738           
  Misses       4008     4008           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

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

[Copilot]

Pull request overview

Adds a Gamma distribution implementation to the stats component, including build integration, tests, and user-facing specifications.

Changes:

• Introduces stdlib_stats_distribution_gamma with rvs_gamma, pdf_gamma, and cdf_gamma.
• Updates build system to compile/link the new stats module (including linking specialfunctions).
• Adds a stats test target and FORD spec documentation + index entry.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 15 comments.

Show a summary per file

File	Description
src/stats/stdlib_stats_distribution_gamma.fypp	Implements Gamma RNG/PDF/CDF (uses specialfunctions for incomplete gamma).
src/stats/CMakeLists.txt	Adds the module to stats sources and links stats against specialfunctions.
test/stats/test_distribution_gamma.fypp	Adds regression tests for RNG/PDF/CDF and a chi-squared RNG check.
test/stats/CMakeLists.txt	Registers the new gamma distribution test.
doc/specs/stdlib_stats_distribution_gamma.md	Adds module specification and examples.
doc/specs/index.md	Adds Gamma distribution spec to the documentation index.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[jalvesz]

Hi @Calluumm, thanks for taking care of this upgrade!
Regarding skipping some of the real kinds, instead of

#:set WITH_QP = False

You should prefer

#:for k, t, s in ...
#:if k in ['sp', 'dp']
...
#:endif
#:endfor

Please notice that the variable you were touching is a global preprocessing variable, which means that if you change it here, whichever other file is treated after will inherit that value. This is very dangerous and difficult to track when fypp preprocessing in parallel.

[Calluumm]

Thanks,
They were relics of the original I wasn't sure how to fix;
I've given that solution a go and it tested fine

[Copilot]

Pull request overview

Copilot reviewed 6 out of 6 changed files in this pull request and generated 7 comments.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Pull request overview

Copilot reviewed 6 out of 6 changed files in this pull request and generated 9 comments.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[jalvesz]

Is there a reason for so many of these functions to be impure ? here for instance, ieee_value and log_gamma can be used within elemental functions. I think it would be good to do a pass to double check the attributes and if they can be promoted to elemental.

[Calluumm]

I left them impure from the original PRs structure. But, the functions largely rely on regamma_p from specialfunctions; that itself is impure elemental, we'd have to change that first if we wanted to change it here.

[jalvesz]

I took a quick glance at the special functions gamma module. It seems the impure attribute was set indeed to enable error stop with a message. This boils down to what should be the design pattern to adopt for the library. Most of those functions seems that would perfectly fit elemental functions. Not that it should be changed but it might be worth discussing.

[Calluumm]

The convention isn't something I'm super familiar with I'm happy to change it to be either.

[jalvesz]

My proposal to discuss would be:

• move numerical kernels to elemental or pure procedures when possible, making sure that NaN is returned for an invalid input.
• Add in the same module an impure version which calls the base kernel but also has an extra intent(out) error argument

idea sketch

module
use, intrinsic :: ieee_arithmetic
implicit none
private

interface my_fun
 module procedure :: my_func_base
 module procedure :: my_func_impure
end interface
public :: my_fun

contains

elemental function my_func_base(x, shape, rate) result(res)
    real, intent(in) :: x, shape, rate
    real :: res
    ! .... crunch numbers
end module

impure elemental function my_func_impure(x, shape, rate, err) result(res)
    real, intent(in) :: x, shape, rate
    integer, intent(out) :: err !> or type(state_type), intent(out) :: err from https://stdlib.fortran-lang.org/page/specs/stdlib_error_state_type.html
    real :: res
    if( ... condition for bad inputs ... ) err = ! some non zero argument or fixed code to retrieve an error message from a static table
   res = my_func_base( x, shape, rate )
end module

Any thoughts @sebastian-mutz, @jvdp1, @perazz ?

[sebastian-mutz]

I like your proposal.

Procedures like distribution functions lend themselves well to be elemental (not "just" pure) and I can see a lot of merit in this – thinking about larger stats/ML applications in an HPC setting, for example. Keeping the core procedures elemental and creating an impure wrapper to handle invalid input and give user feedback (error messages) is the exact approach I took with the FSML distribution module.

To clarify: Are you proposing checking the validity of input in the impure and elemental functions, and handling invalid inputs differently (through NaN returns and error messages, respectively)?

[jalvesz]

Given the elemental procedures have to check in any case to properly return NaN, the impure version could simply catch the NaN and return an appropriate error in that case.

[sebastian-mutz]

That sounds good to me.

[sebastian-mutz]

Hey @Calluumm. Thanks for the PR; nice job on the upgrade and addressing @jalvesz's input.

Note that the original was created before this API discussion and update (#991). To stay consistent with other stats distribution modules, which take the location parameter loc (like scipy, for example), I propose that a loc parameter is added as an input for the various functions in this module. See normal and exponential modules for reference. This may require some more careful thought when updating the core calculations. (Happy to also look at that in detail, of course).

[Calluumm]

@jalvesz
Thanks for the framework, I gave it a go and it tested fine, the above commit is I think how you intended me to try implement it. Do let me know if I've overlooked a detail to your approach

[Calluumm]

To stay consistent with other stats distribution modules, which take the location parameter loc (like scipy, for example), I propose that a loc parameter is added as an input for the various functions in this module. See normal and exponential modules for reference.

@sebastian-mutz looking into this now, are the backwards compatability versions necessary here too? as seen on the exponential/normal modules near the bottom.

[sebastian-mutz]

@sebastian-mutz looking into this now, are the backwards compatability versions necessary here too? as seen on the exponential/normal modules near the bottom.

This was done to not break anyone's work that used the old API. Since the gamma distribution never left the PR stage and never was part of stdlib, you don't have to worry about that.

[Calluumm]

@sebastian-mutz
I've added loc as a paramter but upon updating that in the docs I thought the docs seemed a bit bare compared to the other distribution functions, do you think they need to be more thorough in the descriptions?

[sebastian-mutz]

@sebastian-mutz I've added loc as a paramter but upon updating that in the docs I thought the docs seemed a bit bare compared to the other distribution functions, do you think they need to be more thorough in the descriptions?

Great! That's a good point; it does seem like the gamma dist docs are pretty bare bones in comparison. A more thorough description would certainly be beneficial.

[Calluumm]

Docs should be largely better, templating the existing ones highlighted what gamma was missing