stats.base.dists.gamma.cdf() giving unexpected NaN results

[tlby]

Description

Encountered an error when trying to use stats-base-dists-gamma in node 18. Getting unexpected NaN results, seems like an underflow/overflow issue as truncating values sometimes helps.

Related Issues

Related issues # , # , and # .

Questions

No.

Demo

No response

Reproduction

const gamma = require('@stdlib/stats-base-dists-gamma');

function avg(X) {
    return X.reduce((a, x) => a + x, 0) / X.length
}

function fit(X) {
    // https://en.wikipedia.org/wiki/Gamma_distribution#Closed-form_estimators
    const μ = avg(X);
    const θ = avg(X.map((_) => _ * Math.log(_))) - μ * avg(X.map(Math.log));
    const α = μ / θ;
    const β = 1 / θ;
    return { α, β }
}

function main() {
    const X = [28,29,30,31,32];
    const x = 30;
    const p = fit(X);
    for (let n = 2; n < 16; n++) {
        console.log(n, gamma.cdf(
            Number(x.toFixed(n)),
            Number(p.α.toFixed(n)),
            Number(p.β.toFixed(n)),
        ));
    }
    console.log('∞', gamma.cdf(x, p.α, p.β));
}

main();

/* for reference, compare with scipy:
import scipy
import math

def avg(X):
    X = list(X)
    return sum(X) / len(X)

def fit(X):
    μ = avg(X)
    θ = avg(_ * math.log(_) for _ in X) - μ * avg(math.log(_) for _ in X)
    α = μ / θ
    β = 1 / θ
    return { "α": α, "β": β }

def main():
    X = [28,29,30,31,32]
    x = 30;
    p = fit(X)
    θ = 1 / p["β"]
    for n in range(4, 16):
        print(n, scipy.stats.gamma.cdf(
            round(x, n),
            round(p["α"], n),
            0,
            round(θ, n),
        ))
    print('∞', scipy.stats.gamma.cdf(x, p["α"], 0, θ));

main()
*/

Expected Results

4 0.5062803328241918
5 0.5062746887113254
6 0.5062729601819936
7 0.5062729079101176
8 0.5062730531291298
9 0.5062728061540037
10 0.5062727316430549
11 0.5062728060966877
12 0.5062727282243282
13 0.5062728061601118
14 0.5062728061537761
15 0.5062728061540296
∞ 0.5062728061540007

Actual Results

4 0.5062803328241918
5 0.5062746887113254
6 0.5062729601819936
7 0.5062729079101176
8 NaN
9 0.5062728061540037
10 0.5062727316430549
11 NaN
12 0.5062727282243282
13 NaN
14 NaN
15 NaN
∞ NaN

Version

@stdlib/stats-base-dists-gamma 0.0.7

Environments

Node.js

Browser Version

No response

Node.js / npm Version

18

Platform

Linux

Checklist

• Read and understood the Code of Conduct.
• Searched for existing issues and pull requests.

[stdlib-bot]

👋 Hi there! 👋

And thank you for opening your first issue! We will get back to you shortly. 🏃 💨

[Planeshifter]

@tlby Thanks for bringing this to our attention! We have located the source of the problem and will be looking into a solution shortly. Hopefully we have updates by tomorrow.

[Planeshifter]

With the latest changes pushed to the develop branch of this repo, the above code (thanks @tlby for providing a reproducible example!) now returns the expected result:

2 0.5057083630776964
3 0.5062163591098211
4 0.5062803320345663
5 0.5062746876282798
6 0.5062729754839823
7 0.506272811798204
8 0.506272806718444
9 0.5062728061540037
10 0.5062728061464749
11 0.5062728061515548
12 0.5062728061538885
13 0.5062728061540083
14 0.5062728061540018
15 0.5062728061540018
∞ 0.5062728061540018

We are in the midst of planning a new release cycle of the standalone packages and will post an update to this issue once the updates have been released.

const gamma = require('@stdlib/stats/base/dists/gamma');

function avg(X) {
    return X.reduce((a, x) => a + x, 0) / X.length
}

function fit(X) {
    // https://en.wikipedia.org/wiki/Gamma_distribution#Closed-form_estimators
    const μ = avg(X);
    const θ = avg(X.map((_) => _ * Math.log(_))) - μ * avg(X.map(Math.log));
    const α = μ / θ;
    const β = 1 / θ;
    return { α, β }
}

function main() {
    const X = [28,29,30,31,32];
    const x = 30;
    const p = fit(X);
    for (let n = 2; n < 16; n++) {
        console.log(n, gamma.cdf(
            Number(x.toFixed(n)),
            Number(p.α.toFixed(n)),
            Number(p.β.toFixed(n)),
        ));
    }
    console.log('∞', gamma.cdf(x, p.α, p.β));
}

main();

[Planeshifter]

The newest published version of @stdlib/stats-base-dists-gamma contains the fix. Closing this issue now; thanks again for bringing this to our attention!

[tlby]

Sincere thanks!