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

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

@@ -23,7 +23,24 @@
 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 expm1 = require( '@stdlib/math/base/special/expm1' );
+
+
+// FUNCTIONS //
+
+/**
+* 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 ) {
+		return 1;
+	}
+	return ( 2 * ( expm1( x ) - x ) ) / ( x * x );
+}
 
 
 // MAIN //
@@ -45,10 +62,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 +71,6 @@ function factory( a, b, c ) {
 	) {
 		return constantFunction( NaN );
 	}
-	bmc = b - c;
-	bma = b - a;
-	cma = c - a;
 	return mgf;
 
 	/**
@@ -75,19 +85,24 @@ 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 )
+				);
+			}
+			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,24 @@
 
 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 expm1 = require( '@stdlib/math/base/special/expm1' );
+
+
+// FUNCTIONS //
+
+/**
+* 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 ) {
+		return 1;
+	}
+	return ( 2 * ( expm1( x ) - x ) ) / ( x * x );
+}
 
 
 // MAIN //
@@ -73,11 +90,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 +100,21 @@ 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 )
+			);
+		}
+		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/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,22 @@
 
 #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.
+*
+* @private
+* @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 +49,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/runner.jl

@@ -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" );

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

@@ -34,6 +34,10 @@ var factory = require( './../lib/factory.js' );
 var smallRange = require( './fixtures/julia/small_range.json' );
 var mediumRange = require( './fixtures/julia/medium_range.json' );
 var largeRange = require( './fixtures/julia/large_range.json' );
+var aLessCLessB = require( './fixtures/julia/a_less_c_less_b.json' );
+var aLessCEqB = require( './fixtures/julia/a_less_c_eq_b.json' );
+var aEqCLessB = require( './fixtures/julia/a_eq_c_less_b.json' );
+var aEqCEqB = require( './fixtures/julia/a_eq_c_eq_b.json' );
 
 
 // TESTS //
@@ -224,3 +228,127 @@ tape( 'the created function evaluates the MGF for `x` given a large range `b - a
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the MGF for `x` given the case: a < c < b', function test( t ) {
+	var expected;
+	var delta;
+	var mgf;
+	var tol;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = aLessCLessB.expected;
+	x = aLessCLessB.x;
+	a = aLessCLessB.a;
+	b = aLessCLessB.b;
+	c = aLessCLessB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		mgf = factory( a[i], b[i], c[i] );
+		y = mgf( x[i] );
+		if ( y === expected[i] ) {
+			t.equal( y, expected[i], 'x: '+x[i]+', a: '+a[i]+', b: '+b[i]+', c: '+c[i]+', y: '+y+', expected: '+expected[i] );
+		} else {
+			delta = abs( y - expected[ i ] );
+			tol = 1.0 * EPS * abs( expected[ i ] );
+			t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. a: '+a[i]+'. b: '+b[i]+'. c: '+c[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the MGF for `x` given the case: a < c = b', function test( t ) {
+	var expected;
+	var delta;
+	var mgf;
+	var tol;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = aLessCEqB.expected;
+	x = aLessCEqB.x;
+	a = aLessCEqB.a;
+	b = aLessCEqB.b;
+	c = aLessCEqB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		mgf = factory( a[i], b[i], c[i] );
+		y = mgf( x[i] );
+		if ( y === expected[i] ) {
+			t.equal( y, expected[i], 'x: '+x[i]+', a: '+a[i]+', b: '+b[i]+', c: '+c[i]+', y: '+y+', expected: '+expected[i] );
+		} else {
+			delta = abs( y - expected[ i ] );
+			tol = 1.0 * EPS * abs( expected[ i ] );
+			t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. a: '+a[i]+'. b: '+b[i]+'. c: '+c[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the MGF for `x` given the case: a = c < b', function test( t ) {
+	var expected;
+	var delta;
+	var mgf;
+	var tol;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = aEqCLessB.expected;
+	x = aEqCLessB.x;
+	a = aEqCLessB.a;
+	b = aEqCLessB.b;
+	c = aEqCLessB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		mgf = factory( a[i], b[i], c[i] );
+		y = mgf( x[i] );
+		if ( y === expected[i] ) {
+			t.equal( y, expected[i], 'x: '+x[i]+', a: '+a[i]+', b: '+b[i]+', c: '+c[i]+', y: '+y+', expected: '+expected[i] );
+		} else {
+			delta = abs( y - expected[ i ] );
+			tol = 1.0 * EPS * abs( expected[ i ] );
+			t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. a: '+a[i]+'. b: '+b[i]+'. c: '+c[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the MGF for `x` given the case: a = c = b', function test( t ) {
+	var expected;
+	var delta;
+	var mgf;
+	var tol;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = aEqCEqB.expected;
+	x = aEqCEqB.x;
+	a = aEqCEqB.a;
+	b = aEqCEqB.b;
+	c = aEqCEqB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		mgf = factory( a[i], b[i], c[i] );
+		y = mgf( x[i] );
+		if ( y === expected[i] ) {
+			t.equal( y, expected[i], 'x: '+x[i]+', a: '+a[i]+', b: '+b[i]+', c: '+c[i]+', y: '+y+', expected: '+expected[i] );
+		} else {
+			delta = abs( y - expected[ i ] );
+			tol = 1.0 * EPS * abs( expected[ i ] );
+			t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. a: '+a[i]+'. b: '+b[i]+'. c: '+c[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
+		}
+	}
+	t.end();
+});