Question on the Rayleigh MGF formula — looks like a misplaced parenthesis, is a fix PR welcome?

[AdeshDeshmukh]

Hi @Planeshifter,

I'm Adesh. I spend most of my time between competitive programming and contributing to open source, primarily around cloud-native security tooling (like Kubescape and Harbor under CNCF).

Lately I've been working through stdlib's stats/base/dists packages to understand how the distribution functions are implemented across the JS and C layers. Over the past few days I've been tracing the MGF implementations in particular, checking each one against its closed form.

While reading rayleigh/mgf (lib/main.js, lines 73–75, and the same shape in src/main.c, lines 48–49), one parenthesization caught my attention. My understanding of the Rayleigh MGF is 1 + σt·√(π/2)·e^((σt)²/2)·(1 + erf(σt/√2)), but the code computes (1 + σt·e^((σt)²/2)) · √(π/2)·(1 + erf(σt/√2)) — the trailing factor multiplies the leading 1 as well. If I'm reading it right, mgf(0, σ) would return √(π/2) ≈ 1.2533 instead of 1, and an MGF has to equal 1 at t = 0. The JSDoc example (mgf(1.0, 3.0) → ~678.508) also looks consistent with the parenthesized version rather than the closed form, which is what made me pause.

I cross-checked the neighboring implementations — geometric and negative-binomial mgf, and the kurtosis convention across a dozen sibling distributions — so this looks like an isolated slip rather than a pattern. My rough thought was moving the √(π/2)·(1+erf) factor inside so it multiplies only the σt term, in both the JS and C files, plus updating the doc examples. But I'm unsure of your preferred convention for regenerating those.

Rather than opening an issue or PR directly, I wanted to check in to get your direction. Would a small fix PR with updated examples and test fixtures be welcome, or would you prefer a minimal repro here first? I have the repo checked out locally...!

Open to your guidance on the best way forward here, and really appreciate your time..! :)

[Yuno15-bb]

Your reading is right, and it's in the published package too, not just the source tree:

$ node -e "console.log(require('@stdlib/stats-base-dists-rayleigh-mgf')(0.0, 1.0))"
1.2533141373155003

(0.3.1, Node 26.5.0.) M_X(0) = 1 by definition, so that alone settles it.

Confirmed without picking a closed form

To keep the "which parenthesization" question out of the check, I integrated the definition numerically instead — E[e^{tX}] = ∫₀^∞ e^{tx} · (x/σ²) · e^{−x²/2σ²} dx, mpmath quad at 40 digits — and compared it to both candidates:

(t, σ)	package 0.3.1	your form	numerical ∫
(1.0, 3.0)	678.507933	677.004689	677.004689
(1.0, 2.0)	38.6500976	37.2004954	37.2004954
(-1.0, 4.0)	-0.946530144	0.0533904683	0.0533904683
(1.0, 0.5)	2.7152481	1.98200875	1.98200875
(0.5, 0.5)	1.88782736	1.38709316	1.38709316
(0.0, 1.0)	1.25331414	1.0	1.0

The quadrature matches your form to full precision on every case, so the trailing factor multiplying the leading 1 is the whole of it.

There's also a second symptom worth folding into the same PR: two documented examples are negative, which no MGF can be. mgf( -1.0, 4.0 ) // returns ~-0.947 in lib/main.js#L47-L49, and in the ctor package docs/repl.txt#L89-L90, rayleigh.mgf( -0.5 ) → ~-0.91 with σ = 6. Correct: ~0.053 and ~0.086.

On the fixtures — there are none to regenerate

This is the part I'd flag before you start. All three value tests are loops with no reference values, test/test.mgf.js#L65-L71:

// TODO: Add test fixtures (function is not available in R and Julia)
for ( i = 0; i < 100; i++ ) {
    x = randu();
    sigma = randu() * 20.0;
    y = mgf( x, sigma );
    t.strictEqual( !isnan( y ) && isNumber( y ), true, 'returns expected value' );
}

Same TODO and same assertion in test/test.factory.js#L85-L91 and test/test.native.js#L74-L80. And the ctor test checks the method against the same function, ctor/test/test.js#L288-L290:

y = rayleigh.mgf( 0.5 );
t.strictEqual( y, mgf( 0.5, 1.0 ), 'returns expected value' );

So nothing in the suite pins a number — which is also how this survived — and the fix breaks no test. The TODO's reason still holds today: Distributions.jl's rayleigh.jl defines no mgf method, and R has no MGF in base.

A PR would therefore be creating the first fixtures, with a generator that isn't the usual Julia runner. Two ways that still fit the convention next door (rayleigh/cdf/test/fixtures/julia/{runner.jl, small_scale.json, medium_scale.json, large_scale.json}):

• numerical integration of the definition — QuadGK.quadgk in Julia, or mpmath quad. Independent of the closed form, which is what you want when the closed form is the thing under review;
• high-precision evaluation of the closed form (mpmath at 40 dps) as a cross-check on the first.

I ran both and they agree to ~30 digits on the cases above, so either is defensible; using them together makes the fixtures self-checking. Happy to hand over the runner if useful.

The values that move, from grepping the package at 7ff87de:

value	where
~678.508 → ~677.005, ~38.65 → ~37.2, ~-0.947 → ~0.053	mgf/lib/main.js, lib/native.js, lib/index.js, docs/repl.txt, docs/types/index.d.ts (twice), README.md (twice), src/main.c
~2.715 → ~1.982, ~1.888 → ~1.387 (factory)	mgf/lib/factory.js, docs/types/index.d.ts (twice), README.md
~5.586 → ~4.477	ctor/lib/main.js#L378, ctor/README.md#L224
~-0.91 → ~0.086	ctor/docs/repl.txt#L89-L90

The C is a line-for-line mirror (src/main.c#L47-L50), so the same two lines move there.

What I didn't do: I didn't build the native add-on, so the C path is read-only on my side — but the expression is character-for-character the JS one, so I'd expect it to match. And the direction call is @Planeshifter's, I'm not a maintainer; I just had the numbers from checking your claim.

[AdeshDeshmukh]

@Yuno15-bb Thanks for ur POV... just waiting for maintainers judgement! :)

[Yuno15-bb]

Right ! keep me on

[AdeshDeshmukh]

@kgryte PTAL!