fix: rewrite Wald mode formula in numerically stable rationalized form

Multiplying numerator and denominator by the conjugate avoids
catastrophic cancellation in sqrt(4*lambda^2 + 9*mu^2) - 3*mu when
lambda << mu, dropping max ULP error against the Julia fixtures from
~874 to ~298. Tighten the test tolerance accordingly.

Author: Planeshifter

lib/node_modules/@stdlib/stats/base/dists/wald/mode/lib/main.js

@@ -54,7 +54,6 @@ var sqrt = require( '@stdlib/math/base/special/sqrt' );
 * // returns NaN
 */
 function mode( mu, lambda ) {
-	var numerator;
 	if (
 		isnan( mu ) ||
 		isnan( lambda ) ||
@@ -63,8 +62,7 @@ function mode( mu, lambda ) {
 	) {
 		return NaN;
 	}
-	numerator = sqrt( 4.0*lambda*lambda + 9.0*mu*mu ) - 3.0*mu;
-	return ( mu*numerator ) / ( 2.0*lambda );
+	return ( 2.0*mu*lambda ) / ( sqrt( 4.0*lambda*lambda + 9.0*mu*mu ) + 3.0*mu );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/wald/mode/src/main.c

@@ -32,7 +32,6 @@
 * // returns ~0.325
 */
 double stdlib_base_dists_wald_mode( const double mu, const double lambda ) {
-	double numerator;
 	if (
 		stdlib_base_is_nan( mu ) ||
 		stdlib_base_is_nan( lambda ) ||
@@ -41,6 +40,5 @@ double stdlib_base_dists_wald_mode( const double mu, const double lambda ) {
 	) {
 		return 0.0/0.0; // NaN
 	}
-	numerator = stdlib_base_sqrt( 4.0*lambda*lambda + 9.0*mu*mu ) - 3.0*mu;
-	return mu * numerator / ( 2.0*lambda );
+	return ( 2.0*mu*lambda ) / ( stdlib_base_sqrt( 4.0*lambda*lambda + 9.0*mu*mu ) + 3.0*mu );
 }

lib/node_modules/@stdlib/stats/base/dists/wald/mode/test/test.js

@@ -117,7 +117,7 @@ tape( 'the function returns the mode of a Wald distribution', function test( t )
 		if ( y === expected[i] ) {
 			t.strictEqual( y, expected[i], 'mu:'+mu[i]+', lambda: '+lambda[i]+', y: '+y+', expected: '+expected[i] );
 		} else {
-			t.ok( isAlmostSameValue( y, expected[i], 1000 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
+			t.ok( isAlmostSameValue( y, expected[i], 500 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
 		}
 	}
 	t.end();

lib/node_modules/@stdlib/stats/base/dists/wald/mode/test/test.native.js

@@ -126,7 +126,7 @@ tape( 'the function returns the mode of a Wald distribution', opts, function tes
 		if ( y === expected[i] ) {
 			t.strictEqual( y, expected[i], 'mu:'+mu[i]+', lambda: '+lambda[i]+', y: '+y+', expected: '+expected[i] );
 		} else {
-			t.ok( isAlmostSameValue( y, expected[i], 1000 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
+			t.ok( isAlmostSameValue( y, expected[i], 500 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
 		}
 	}
 	t.end();