refactor: update stats/base/dists/triangular/mgf implementation and test tolerances

PR-URL: #4768
Closes: #4704

Co-authored-by: Philipp Burckhardt <pburckhardt@outlook.com>
Co-authored-by: stdlib-bot <noreply@stdlib.io>
Reviewed-by: Philipp Burckhardt <pburckhardt@outlook.com>
Reviewed-by: Athan Reines <kgryte@gmail.com>
Signed-off-by: Philipp Burckhardt <pburckhardt@outlook.com>

Author: 3 people

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/benchmark/benchmark.js

@@ -22,7 +22,7 @@
 
 var bench = require( '@stdlib/bench' );
 var Float64Array = require( '@stdlib/array/float64' );
-var randu = require( '@stdlib/random/base/randu' );
+var uniform = require( '@stdlib/random/base/uniform' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var EPS = require( '@stdlib/constants/float64/eps' );
 var pkg = require( './../package.json' ).name;
@@ -47,10 +47,10 @@ bench( pkg, function benchmark( b ) {
 	mode = new Float64Array( len );
 
 	for ( i = 0; i < len; i++ ) {
-		t[ i ] = randu() * 5.0;
-		min[ i ] = randu() * 10.0;
-		max[ i ] = min[ i ] + ( randu() * 40.0 ) + EPS;
-		mode[ i ] = min[ i ] + ( ( max[ i ] - min[ i ] ) * randu() );
+		t[ i ] = uniform( 0.0, 5.0 );
+		min[ i ] = uniform( 0.0, 10.0 );
+		max[ i ] = uniform( min[ i ] + EPS, min[ i ] + 40.0 );
+		mode[ i ] = uniform( min[ i ], max[ i ]);
 	}
 
 	b.tic();
@@ -73,6 +73,7 @@ bench( pkg+':factory', function benchmark( b ) {
 	var mode;
 	var min;
 	var max;
+	var len;
 	var t;
 	var y;
 	var i;
@@ -81,11 +82,15 @@ bench( pkg+':factory', function benchmark( b ) {
 	max = 1.5;
 	mode = 0.5;
 	mymgf = mgf.factory( min, max, mode );
+	len = 100;
+	t = new Float64Array( len );
+	for ( i = 0; i < len; i++ ) {
+		t[ i ] = uniform( 0.0, 5.0 );
+	}
 
 	b.tic();
 	for ( i = 0; i < b.iterations; i++ ) {
-		t = randu() * 5.0;
-		y = mymgf( t );
+		y = mymgf( t[ i % len ]);
 		if ( isnan( y ) ) {
 			b.fail( 'should not return NaN' );
 		}

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/benchmark/benchmark.native.js

@@ -23,7 +23,7 @@
 var resolve = require( 'path' ).resolve;
 var bench = require( '@stdlib/bench' );
 var Float64Array = require( '@stdlib/array/float64' );
-var randu = require( '@stdlib/random/base/randu' );
+var uniform = require( '@stdlib/random/base/uniform' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var tryRequire = require( '@stdlib/utils/try-require' );
 var EPS = require( '@stdlib/constants/float64/eps' );
@@ -56,10 +56,10 @@ bench( pkg+'::native', opts, function benchmark( b ) {
 	mode = new Float64Array( len );
 
 	for ( i = 0; i < len; i++ ) {
-		t[ i ] = randu() * 5.0;
-		min[ i ] = randu() * 10.0;
-		max[ i ] = min[ i ] + ( randu() * 40.0 ) + EPS;
-		mode[ i ] = min[ i ] + ( ( max[ i ] - min[ i ] ) * randu() );
+		t[ i ] = uniform( 0.0, 5.0 );
+		min[ i ] = uniform( 0.0, 10.0 );
+		max[ i ] = uniform( min[ i ] + EPS, min[ i ] + 40.0 );
+		mode[ i ] = uniform( min[ i ], max[ i ]);
 	}
 
 	b.tic();

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/lib/factory.js

@@ -23,7 +23,7 @@
 var constantFunction = require( '@stdlib/utils/constant-function' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var exp = require( '@stdlib/math/base/special/exp' );
-var pow = require( '@stdlib/math/base/special/pow' );
+var phi2 = require( './phi2.js' );
 
 
 // MAIN //
@@ -45,10 +45,6 @@ var pow = require( '@stdlib/math/base/special/pow' );
 * // returns ~10.205
 */
 function factory( a, b, c ) {
-	var bmc;
-	var bma;
-	var cma;
-
 	if (
 		isnan( a ) ||
 		isnan( b ) ||
@@ -58,9 +54,6 @@ function factory( a, b, c ) {
 	) {
 		return constantFunction( NaN );
 	}
-	bmc = b - c;
-	bma = b - a;
-	cma = c - a;
 	return mgf;
 
 	/**
@@ -75,19 +68,19 @@ function factory( a, b, c ) {
 	* // returns <number>
 	*/
 	function mgf( t ) {
-		var ret;
-
 		if ( isnan( t ) ) {
 			return NaN;
 		}
-		if ( t === 0.0 ) {
-			return 1.0;
+		if ( a < c ) {
+			if ( c < b ) {
+				return exp( c*t ) * ( ( (c-a)*phi2( (a-c)*t ) ) + ( (b-c)*phi2( (b-c)*t ) ) ) / ( b-a ); // eslint-disable-line max-len
+			}
+			return exp( c*t ) * phi2( ( a-c ) * t );
+		}
+		if ( c < b ) {
+			return exp( c*t ) * phi2( ( b-c ) * t );
 		}
-		ret = ( bmc * exp( a*t ) ) - ( bma * exp( c*t ) );
-		ret += cma * exp( b*t );
-		ret *= 2.0;
-		ret /= bma * cma * bmc * pow( t, 2.0 );
-		return ret;
+		return exp( c*t );
 	}
 }
 

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/lib/main.js

@@ -22,7 +22,7 @@
 
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var exp = require( '@stdlib/math/base/special/exp' );
-var pow = require( '@stdlib/math/base/special/pow' );
+var phi2 = require( './phi2.js' );
 
 
 // MAIN //
@@ -73,11 +73,6 @@ var pow = require( '@stdlib/math/base/special/pow' );
 * // returns NaN
 */
 function mgf( t, a, b, c ) {
-	var bmc;
-	var bma;
-	var cma;
-	var ret;
-
 	if (
 		isnan( t ) ||
 		isnan( a ) ||
@@ -88,17 +83,16 @@ function mgf( t, a, b, c ) {
 	) {
 		return NaN;
 	}
-	if ( t === 0.0 ) {
-		return 1.0;
+	if ( a < c ) {
+		if ( c < b ) {
+			return exp( c*t ) * ( ( (c-a)*phi2( (a-c)*t ) ) + ( (b-c)*phi2( (b-c)*t ) ) ) / ( b-a ); // eslint-disable-line max-len
+		}
+		return exp( c*t ) * phi2( ( a-c ) * t );
+	}
+	if ( c < b ) {
+		return exp( c*t ) * phi2( ( b-c ) * t );
 	}
-	bmc = b - c;
-	bma = b - a;
-	cma = c - a;
-	ret = ( bmc * exp( a*t ) ) - ( bma * exp( c*t ) );
-	ret += cma * exp( b*t );
-	ret *= 2.0;
-	ret /= bma * cma * bmc * pow( t, 2.0 );
-	return ret;
+	return exp( c*t );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/lib/phi2.js

@@ -0,0 +1,45 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2025 The Stdlib Authors.
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*    http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+'use strict';
+
+// MODULES //
+
+var expm1 = require( '@stdlib/math/base/special/expm1' );
+
+
+// MAIN //
+
+/**
+* Helper function for repeated computation in the MGF formula.
+*
+* @private
+* @param {number} x - input value
+* @returns {number} evaluated result
+*/
+function phi2( x ) {
+	if ( x === 0.0 ) {
+		return 1.0;
+	}
+	return ( 2.0 * ( expm1( x ) - x ) ) / ( x*x );
+}
+
+
+// EXPORTS //
+
+module.exports = phi2;

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/manifest.json

@@ -41,7 +41,7 @@
         "@stdlib/math/base/napi/quaternary",
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/exp",
-        "@stdlib/math/base/special/pow"
+        "@stdlib/math/base/special/expm1"
       ]
     },
     {

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/src/main.c

@@ -18,8 +18,21 @@
 
 #include "stdlib/stats/base/dists/triangular/mgf.h"
 #include "stdlib/math/base/assert/is_nan.h"
+#include "stdlib/math/base/special/expm1.h"
 #include "stdlib/math/base/special/exp.h"
-#include "stdlib/math/base/special/pow.h"
+
+/**
+* Helper function for repeated computation in the MGF formula.
+*
+* @param x   input value
+* @return    evaluated result
+*/
+static double phi2( const double x ) {
+	if ( x == 0.0 ) {
+		return 1.0;
+	}
+	return ( 2.0 * ( stdlib_base_expm1( x ) - x ) ) / ( x*x );
+}
 
 /**
 * Evaluates the moment-generating function (MGF) for a triangular distribution with lower limit `a`, upper limit `b`, and mode `c` at a value `t`.
@@ -35,29 +48,19 @@
 * // returns ~1.021
 */
 double stdlib_base_dists_triangular_mgf( const double t, const double a, const double b, const double c ) {
-	double bmc;
-	double bma;
-	double cma;
-	double ret;
-	if (
-		stdlib_base_is_nan( t ) ||
-		stdlib_base_is_nan( a ) ||
-		stdlib_base_is_nan( b ) ||
-		stdlib_base_is_nan( c ) ||
-		a > c ||
-		c > b
-	) {
-		return 0.0/0.0; // NaN
+	if ( stdlib_base_is_nan( t ) || stdlib_base_is_nan( a ) || stdlib_base_is_nan( b ) || stdlib_base_is_nan( c ) || a > c || c > b ) {
+		return 0.0 / 0.0; // NaN
 	}
-	if ( t == 0.0 ) {
-		return 1.0;
+	if ( a < c ) {
+		if ( c < b ) {
+			double term1 = ( c-a ) * phi2( ( a-c ) * t );
+			double term2 = ( b-c ) * phi2( ( b-c ) * t );
+			return stdlib_base_exp( c*t ) * ( term1+term2 ) / ( b-a );
+		}
+		return stdlib_base_exp( c*t ) * phi2( ( a-c ) * t );
+	}
+	if ( c < b ) {
+		return stdlib_base_exp( c*t ) * phi2( ( b-c ) * t );
 	}
-	bmc = b - c;
-	bma = b - a;
-	cma = c - a;
-	ret = ( bmc * stdlib_base_exp( a*t ) ) - ( bma * stdlib_base_exp( c*t ) );
-	ret += cma * stdlib_base_exp( b*t );
-	ret *= 2.0;
-	ret /= bma * cma * bmc * stdlib_base_pow( t, 2.0 );
-	return ret;
+	return stdlib_base_exp( c*t );
 }

lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/test/fixtures/julia/a_eq_c_eq_b.json

@@ -93,3 +93,15 @@ b = ( rand( 1000 ) .* 80.0 ) .+ a;
 c = a .+ ( b .- a ) .* rand();
 x = rand( 1000 );
 gen( x, a, b, c, "large_range.json" );
+
+# Case: a < c < b
+gen( x, a, b, c, "a_less_c_less_b.json" );
+
+# Case: a < c = b
+gen( x, a, b, b, "a_less_c_eq_b.json" );
+
+# Case: a = c < b
+gen( x, a, b, a, "a_eq_c_less_b.json" );
+
+# Case: a = c = b
+gen( x, a, a, a, "a_eq_c_eq_b.json" );