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 //
@@ -44,23 +61,16 @@ var pow = require( '@stdlib/math/base/special/pow' );
 * y = mgf( 2.0 );
 * // returns ~10.205
 */
-function factory( a, b, c ) {
-	var bmc;
-	var bma;
-	var cma;
-
+function factory(a, b, c) {
 	if (
-		isnan( a ) ||
-		isnan( b ) ||
-		isnan( c ) ||
+		isnan(a) ||
+		isnan(b) ||
+		isnan(c) ||
 		a > c ||
 		c > b
 	) {
-		return constantFunction( NaN );
+		return constantFunction(NaN);
 	}
-	bmc = b - c;
-	bma = b - a;
-	cma = c - a;
 	return mgf;
 
 	/**
@@ -74,20 +84,28 @@ function factory( a, b, c ) {
 	* var y = mgf( 0.5 );
 	* // returns <number>
 	*/
-	function mgf( t ) {
-		var ret;
-
-		if ( isnan( t ) ) {
+	function mgf(t) {
+		if (isnan(t)) {
 			return NaN;
 		}
-		if ( t === 0.0 ) {
+		if (t === 0.0) {
 			return 1.0;
 		}
-		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;
+		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);
+		}
+		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 //
@@ -72,33 +89,35 @@ var pow = require( '@stdlib/math/base/special/pow' );
 * var y = mgf( 0.5, 1.0, 0.0, 1.5 );
 * // returns NaN
 */
-function mgf( t, a, b, c ) {
-	var bmc;
-	var bma;
-	var cma;
-	var ret;
-
+function mgf(t, a, b, c) {
 	if (
-		isnan( t ) ||
-		isnan( a ) ||
-		isnan( b ) ||
-		isnan( c ) ||
+		isnan(t) ||
+		isnan(a) ||
+		isnan(b) ||
+		isnan(c) ||
 		a > c ||
 		c > b
 	) {
 		return NaN;
 	}
-	if ( t === 0.0 ) {
+	if (t === 0.0) {
 		return 1.0;
 	}
-	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;
+	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);
+	}
+	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

@@ -16,10 +16,24 @@
 * limitations under the License.
 */
 
-#include "stdlib/stats/base/dists/triangular/mgf.h"
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/exp.h"
-#include "stdlib/math/base/special/pow.h"
+#include "stdlib/math/base/special/expm1.h"
+#include "stdlib/stats/base/dists/triangular/mgf.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,22 @@
 * // 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;
 	}
-	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;
+	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 );
+	}
+	return stdlib_base_exp( c * t );
 }

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

@@ -103,9 +103,7 @@ tape( 'the function evaluates the MGF for `x` given a small range `b - a`', opts
 			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 ] );
-
-			// FIXME: the tolerance for the C implementation is widely different from the JavaScript implementation, and it is not clear why.
-			tol = 1.0e7 * EPS * abs( expected[ i ] );
+			tol = 4.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+'.' );
 		}
 	}
@@ -134,9 +132,7 @@ tape( 'the function evaluates the MGF for `x` given a medium range `b - a`', opt
 			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 ] );
-
-			// FIXME: the tolerance for the C implementation is widely different from the JavaScript implementation, and it is not clear why.
-			tol = 1.0e7 * EPS * abs( expected[ i ] );
+			tol = 4.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+'.' );
 		}
 	}
@@ -165,9 +161,7 @@ tape( 'the function evaluates the MGF for `x` given a large range `b - a`', opts
 			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 ] );
-
-			// FIXME: the tolerance for the C implementation is widely different from the JavaScript implementation, and it is not clear why.
-			tol = 1.0e7 * EPS * abs( expected[ i ] );
+			tol = 4.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+'.' );
 		}
 	}