[Bug]: rayleigh/mgf returns wrong values due to misplaced parenthesis

Description

I encountered incorrect output from @stdlib/stats/base/dists/rayleigh/mgf while cross-checking the MGF implementations against their closed forms. The Rayleigh implementation multiplies the whole (1 + ...) expression by the √(π/2)·(1 + erf) factor, when that factor should only multiply the σt term.

The current code in lib/main.js (and identically in lib/factory.js and src/main.c) reads as follows:

sigmat = t * sigma;
out = 1.0 + (sigmat * exp( sigmat*sigmat / 2.0 ));
out *= SQRT_HALF_PI * ( erf( sigmat / SQRT2 ) + 1.0 );

The closed form for the Rayleigh MGF is M(t) = 1 + σt·√(π/2)·e^((σt)²/2)·(1 + erf(σt/√2)), so the trailing factor ends up being applied to the leading 1 as well. The quickest way to see the problem is that mgf(0, σ) returns √(π/2) ≈ 1.2533 instead of 1, and every moment-generating function must equal 1 at t = 0. The error also flips signs in places: mgf(-1, 4) currently returns -0.947, whereas the correct value is +0.053.
The JSDoc and README examples (mgf(1.0, 3.0) → ~678.508 and similar) match the buggy formula, so the incorrect values are baked into the documentation as well. I checked the neighboring mgf implementations (geometric, negative-binomial, Poisson, gamma, normal, and others), and this problem looks isolated to Rayleigh.

Proposed Approach

The fix is a small regrouping of the existing expression so that the √(π/2)·(1 + erf) factor multiplies only the σt term, applied in all three places that carry the formula:

sigmat = t * sigma;
out = 1.0 + ( sigmat * SQRT_HALF_PI * exp( sigmat*sigmat / 2.0 ) * ( erf( sigmat / SQRT2 ) + 1.0 ) );

Concretely, I propose updating lib/main.js, lib/factory.js, and src/main.c, along with the corresponding doc examples in README.md, docs/repl.txt, docs/types/index.d.ts, lib/index.js, and lib/native.js. The corrected example outputs, verified against stdlib's own erf implementation, are mgf(1.0, 3.0) → ~677.005, mgf(1.0, 2.0) → ~37.2, mgf(-1.0, 4.0) → ~0.053, mgf.factory(0.5)(1.0) → ~1.982, and mgf.factory(0.5)(0.5) → ~1.387. I have these changes ready locally and fully verified .

I’m happy to prepare a PR for this if the team is aligned , however please let me know if there’s another direction you'd prefer to take...! :)

Related Issues

Related to my earlier discussion #15298

Questions

I have no blocking questions. I would, however, appreciate guidance on whether the pull request should also add numeric fixtures to test.mgf.js and test.factory.js, since those tests currently only assert non-NaN output, which is how this slipped through unnoticed.

Demo

N/A

Reproduction

var mgf = require( '@stdlib/stats/base/dists/rayleigh/mgf' );
mgf( 0.0, 1.0 ); // returns ~1.253 instead of 1
mgf( 1.0, 3.0 ); // returns ~678.508 instead of ~677.005
mgf( -1.0, 4.0 ); // returns ~-0.947 instead of ~0.053

Expected Results

1
677.0046886610472
0.05339046834575567

Actual Results

1.2533141373155003
678.5079332431045
-0.9465301436239729

Version

0.4.1 (develop)

Environments

Node.js

Browser Version

N/A

Node.js / npm Version

Node v26.5.0 (the numbers were verified against stdlib's own erf, so the discrepancy lies in the formula rather than the environment).

Platform

macOS (darwin). This is a formula-level bug and is platform-independent.

Checklist

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

👋 Hi there! 👋

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

If you have any further questions while waiting for a response, please join our Zulip community to chat with project maintainers and other community members.

@anandkaranubc

@anandkaranubc Would you mind investigating?

@AdeshDeshmukh

Of course 🫡
Thanks @AdeshDeshmukh for bringing this up!

@AdeshDeshmukh

@AdeshDeshmukh Confirmed, this is a real bug. I integrated the MGF definition directly and your numbers check out on every case. mgf( 0.0, 1.0 ) returning 1.253 instead of 1 is the giveaway. The C implementation has the same two lines, so it's wrong there too. It's been like this since 2017, and the tests only assert that the output is a number, which is how I think it slipped through.

One change to the fix though. Moving the parentheses is mathematically correct but falls apart for negative t. Around σt = -10, 1 + erf( σt/√2 ) rounds to exactly 0 while exp( σ²t²/2 ) is enormous, so you get back 1 when the real answer is ~0.0097. At σt = -40 it's NaN.

Since e^(σ²t²/2)·(1 + erf(σt/√2)) is just erfcx( -σt/√2 ), the whole thing collapses to

out = 1.0 + ( sigmat * SQRT_HALF_PI * erfcx( -sigmat / SQRT2 ) );

Tested this locally in JS and C against ~3000 fixture points with σt in [-20, 20]. No overflow or cancellation issues anywhere, and it drops the exp and erf deps. wald/cdf and normal/logcdf already use erfcx for the same reason.

File list: you have everything except the namespace-level stats/base/dists/rayleigh/docs/types/index.d.ts, which has the same stale values. The ctor package also needs ~5.586 → ~4.477 and ~-0.91 → ~0.086.

Fixtures: yes please. Julia/R/SciPy don't have a Rayleigh MGF, but a Python runner that integrates the definition with scipy.integrate.quad works fine. planck/mgf is a template for Python fixtures. Make sure some points have negative t, otherwise the old code would have passed too.

For reference, here's the change plus the fixture runner: https://gist.github.com/anandkaranubc/18b27f728c9cb96b63c143f24b50c273

@anandkaranubc

Hi @anandkaranubc , thank you for the thorough review ... and for the gist , it made the path very clear...!

I verified everything on my end first. You were right about the naive fix: at σt = −10 it returns exactly 1 (true value ~0.0097, confirmed by integrating the definition directly), and at σt = −40 it goes to NaN. The erfcx formulation holds across the full sweep, so I have implemented your approach as specified — code in all three places, manifest deps, the two missed doc spots, your runner.py verbatim with all three scipy fixture files (~50% negative-t coverage, cross-checked to ~1e-12), and test.mgf.js per your diff.

One honest flag : a single small-scale fixture point sits at 54 ULP against the 50 in your diff — both sides within ~1e-14 of truth, so it reads as fixture quadrature noise. I set that tolerance to 100 ULP to keep the suite green; happy to adjust if you see it differently...!

I will open the PR shortly .
really appreciate you working through this with me...!

fix: resolve incorrect Rayleigh MGF values via erfcx formulation

Sunbelt Computer Software

PL/B Language Development and Support

Navigation Menu

Repository navigation

[Bug]: rayleigh/mgf returns wrong values due to misplaced parenthesis #15456

Description

Proposed Approach

Related Issues

Questions

Demo

Reproduction

Expected Results

Actual Results

Version

Environments

Browser Version

Node.js / npm Version

Platform

Checklist

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Sunbelt Computer Software

PL/B Language Development and Support

[Bug]: rayleigh/mgf returns wrong values due to misplaced parenthesis #15456

Description

Description

Proposed Approach

Related Issues

Questions

Demo

Reproduction

Expected Results

Actual Results

Version

Environments

Browser Version

Node.js / npm Version

Platform

Checklist

Activity

stdlib-bot commented on Sep 23, 2026

kgryte commented on Sep 25, 2026

anandkaranubc commented on Sep 25, 2026

anandkaranubc commented on Sep 27, 2026

anandkaranubc commented on Sep 27, 2026

AdeshDeshmukh commented on Sep 27, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions