fix: avoid overflow in stats/base/dists/poisson/entropy for large lambda

lib/node_modules/@stdlib/stats/base/dists/poisson/entropy/lib/main.js

@@ -20,41 +20,44 @@
 
 // MODULES //
 
-var sumSeries = require( '@stdlib/math/base/tools/sum-series' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
-var factorialln = require( '@stdlib/math/base/special/factorialln' );
-var factorial = require( '@stdlib/math/base/special/factorial' );
+var floor = require( '@stdlib/math/base/special/floor' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
 var exp = require( '@stdlib/math/base/special/exp' );
 var ln = require( '@stdlib/math/base/special/ln' );
+var LN_TWO_PI = require( '@stdlib/constants/float64/ln-two-pi' );
+var EPS = require( '@stdlib/constants/float64/eps' );
+
+
+// VARIABLES //
+
+// Value of `lambda` above which the asymptotic expansion is accurate to within double-precision (relative error < 1e-16):
+var THRESHOLD = 5000.0;
+
+// Relative weight below which the probability mass of a term no longer affects the sum:
+var TOL = EPS * 1.0e-3;
 
 
 // FUNCTIONS //
 
 /**
-* Returns a function to retrieve elements of the series \\( \sum_{k=0}^{\infty} \frac{ \lambda^k \log(k!) }{ k! } \\).
+* Returns `ln(a/b)` for positive `a` and `b`.
+*
+* ## Notes
+*
+* -   When `a` and `b` are within a factor of two of each other, `a-b` is exact (Sterbenz lemma), and `log1p` avoids the cancellation which `ln(a/b)` suffers when `a/b` is close to one.
 *
 * @private
-* @param {NonNegativeNumber} lambda - mean parameter
-* @returns {Function} function to retrieve series elements
+* @param {PositiveNumber} a - numerator
+* @param {PositiveNumber} b - denominator
+* @returns {number} natural logarithm of the ratio
 */
-function seriesClosure( lambda ) {
-	var lk;
-	var k;
-	k = 1;
-	lk = lambda;
-	return seriesElement;
-
-	/**
-	* Returns the current series element.
-	*
-	* @private
-	* @returns {number} series element
-	*/
-	function seriesElement() {
-		k += 1;
-		lk *= lambda;
-		return lk * factorialln( k ) / factorial( k );
+function lnratio( a, b ) {
+	var r = a / b;
+	if ( r > 0.5 && r < 2.0 ) {
+		return log1p( ( a-b ) / b );
 	}
+	return ln( r );
 }
 
 
@@ -63,6 +66,24 @@ function seriesClosure( lambda ) {
 /**
 * Returns the entropy of a Poisson distribution.
 *
+* ## Method
+*
+* For \\(\lambda < 5000\\), the function evaluates the definition \\(H = -\sum_k p_k \ln p_k\\) using weights relative to the mode \\(m = \lfloor \lambda \rfloor\\), \\(w_k = p_k / p_m\\). As the probabilities sum to one, \\(p_m = 1/W\\), where \\(W = \sum_k w_k\\), and thus
+*
+* ```tex
+* H = \ln W - \frac{1}{W} \sum_k w_k \ln w_k
+* ```
+*
+* The logarithms \\(\ln w_k\\) are accumulated using the recurrences \\(w_{k+1} = w_k \lambda/(k+1)\\) and \\(w_{k-1} = w_k k/\lambda\\). This avoids evaluating \\(\lambda^k/k!\\), which overflows for large \\(k\\), and avoids the cancellation between \\(k \ln \lambda\\), \\(\lambda\\), and \\(\ln k!\\), each of which is much larger than the entropy.
+*
+* For \\(\lambda \geq 5000\\), the function uses the asymptotic expansion
+*
+* ```tex
+* H = \frac{1}{2} \ln(2 \pi e \lambda) - \frac{1}{12\lambda} - \frac{1}{24\lambda^2} - \frac{19}{360\lambda^3} + O\left(\lambda^{-4}\right)
+* ```
+*
+* whose truncation error is less than \\(10^{-16}\\) relative to the result.
+*
 * @param {NonNegativeNumber} lambda - mean parameter
 * @returns {PositiveNumber} entropy
 *
@@ -83,18 +104,50 @@ function seriesClosure( lambda ) {
 * // returns NaN
 */
 function entropy( lambda ) {
-	var gen;
-	var out;
+	var sum;
+	var wm1;
+	var lw;
+	var l2;
+	var w;
+	var m;
+	var k;
 	if ( isnan( lambda ) || lambda < 0.0 ) {
 		return NaN;
 	}
 	if ( lambda === 0.0 ) {
 		return 0.0;
 	}
-	gen = seriesClosure( lambda );
-	out = lambda * ( 1.0-ln(lambda) );
-	out += exp( -lambda ) * sumSeries( gen );
-	return out;
+	if ( lambda >= THRESHOLD ) {
+		l2 = lambda * lambda;
+		return ( 0.5 * ( LN_TWO_PI + 1.0 + ln( lambda ) ) ) - ( 1.0 / ( 12.0*lambda ) ) - ( 1.0 / ( 24.0*l2 ) ) - ( 19.0 / ( 360.0*l2*lambda ) ); // eslint-disable-line max-len
+	}
+	m = floor( lambda );
+	wm1 = 0.0; // W - 1
+	sum = 0.0; // Σ w_k ln(w_k)
+
+	// Terms above the mode, using w_{k+1} = w_k * λ/(k+1):
+	lw = 0.0;
+	for ( k = m; ; k++ ) {
+		lw += lnratio( lambda, k+1 );
+		w = exp( lw );
+		wm1 += w;
+		sum += w * lw;
+		if ( w < TOL ) {
+			break;
+		}
+	}
+	// Terms below the mode, using w_{k-1} = w_k * k/λ:
+	lw = 0.0;
+	for ( k = m; k > 0; k-- ) {
+		lw += lnratio( k, lambda );
+		w = exp( lw );
+		wm1 += w;
+		sum += w * lw;
+		if ( w < TOL ) {
+			break;
+		}
+	}
+	return log1p( wm1 ) - ( sum / ( 1.0+wm1 ) );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/poisson/entropy/manifest.json

@@ -40,10 +40,12 @@
       "dependencies": [
         "@stdlib/math/base/napi/unary",
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/factorialln",
-        "@stdlib/math/base/special/factorial",
+        "@stdlib/math/base/special/floor",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/math/base/special/exp",
-        "@stdlib/math/base/special/ln"
+        "@stdlib/math/base/special/ln",
+        "@stdlib/constants/float64/ln-two-pi",
+        "@stdlib/constants/float64/eps"
       ]
     },
     {
@@ -59,10 +61,12 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/factorialln",
-        "@stdlib/math/base/special/factorial",
+        "@stdlib/math/base/special/floor",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/math/base/special/exp",
         "@stdlib/math/base/special/ln",
+        "@stdlib/constants/float64/ln-two-pi",
+        "@stdlib/constants/float64/eps",
         "@stdlib/constants/float64/eps"
       ]
     },
@@ -79,10 +83,12 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/factorialln",
-        "@stdlib/math/base/special/factorial",
+        "@stdlib/math/base/special/floor",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/math/base/special/exp",
-        "@stdlib/math/base/special/ln"
+        "@stdlib/math/base/special/ln",
+        "@stdlib/constants/float64/ln-two-pi",
+        "@stdlib/constants/float64/eps"
       ]
     }
   ]

lib/node_modules/@stdlib/stats/base/dists/poisson/entropy/src/main.c

@@ -18,14 +18,59 @@
 
 #include "stdlib/stats/base/dists/poisson/entropy.h"
 #include "stdlib/math/base/assert/is_nan.h"
-#include "stdlib/math/base/special/factorialln.h"
-#include "stdlib/math/base/special/factorial.h"
+#include "stdlib/math/base/special/floor.h"
+#include "stdlib/math/base/special/log1p.h"
 #include "stdlib/math/base/special/exp.h"
 #include "stdlib/math/base/special/ln.h"
+#include "stdlib/constants/float64/ln_two_pi.h"
+#include "stdlib/constants/float64/eps.h"
+
+// Value of `lambda` above which the asymptotic expansion is accurate to within double-precision (relative error < 1e-16):
+static const double THRESHOLD = 5000.0;
+
+// Relative weight below which the probability mass of a term no longer affects the sum:
+static const double TOL = STDLIB_CONSTANT_FLOAT64_EPS * 1.0e-3;
+
+/**
+* Returns `ln(a/b)` for positive `a` and `b`.
+*
+* ## Notes
+*
+* -   When `a` and `b` are within a factor of two of each other, `a-b` is exact (Sterbenz lemma), and `log1p` avoids the cancellation which `ln(a/b)` suffers when `a/b` is close to one.
+*
+* @param a    numerator
+* @param b    denominator
+* @return     natural logarithm of the ratio
+*/
+static double lnratio( const double a, const double b ) {
+	double r = a / b;
+	if ( r > 0.5 && r < 2.0 ) {
+		return stdlib_base_log1p( ( a-b ) / b );
+	}
+	return stdlib_base_ln( r );
+}
 
 /**
 * Returns the entropy of a Poisson distribution.
 *
+* ## Method
+*
+* For `\lambda < 5000`, the function evaluates the definition `H = -\sum_k p_k \ln p_k` using weights relative to the mode `m = \lfloor \lambda \rfloor`, `w_k = p_k / p_m`. As the probabilities sum to one, `p_m = 1/W`, where `W = \sum_k w_k`, and thus
+*
+* ```tex
+* H = \ln W - \frac{1}{W} \sum_k w_k \ln w_k
+* ```
+*
+* The logarithms `\ln w_k` are accumulated using the recurrences `w_{k+1} = w_k \lambda/(k+1)` and `w_{k-1} = w_k k/\lambda`. This avoids evaluating `\lambda^k/k!`, which overflows for large `k`, and avoids the cancellation between `k \ln \lambda`, `\lambda`, and `\ln k!`, each of which is much larger than the entropy.
+*
+* For `\lambda \geq 5000`, the function uses the asymptotic expansion
+*
+* ```tex
+* H = \frac{1}{2} \ln(2 \pi e \lambda) - \frac{1}{12\lambda} - \frac{1}{24\lambda^2} - \frac{19}{360\lambda^3} + O\left(\lambda^{-4}\right)
+* ```
+*
+* whose truncation error is less than `10^{-16}` relative to the result.
+*
 * @param lambda    mean parameter
 * @return          entropy
 *
@@ -34,34 +79,49 @@
 * // returns ~2.508
 */
 double stdlib_base_dists_poisson_entropy( const double lambda ) {
-	double term;
 	double sum;
-	double out;
-	double lk;
-	int k;
+	double wm1;
+	double lw;
+	double l2;
+	double w;
+	double m;
+	double k;
 
 	if ( stdlib_base_is_nan( lambda ) || lambda < 0.0 ) {
 		return 0.0 / 0.0; // NaN
 	}
 	if ( lambda == 0.0 ) {
 		return 0.0;
 	}
-	out = lambda * ( 1.0 - stdlib_base_ln( lambda ) );
-
-	// Compute: sum_{k=2}^{inf} lambda^k * ln(k!) / k!
-	// k starts at 1, lk = lambda^1, first iteration makes k=2, lk=lambda^2
-	sum = 0.0;
-	lk = lambda; // lambda^1
-	for ( k = 2; k < 1000; k++ ) {
-		lk *= lambda;
-		term = lk * stdlib_base_factorialln( k ) / stdlib_base_factorial( k );
-		sum += term;
+	if ( lambda >= THRESHOLD ) {
+		l2 = lambda * lambda;
+		return ( 0.5 * ( STDLIB_CONSTANT_FLOAT64_LN_TWO_PI + 1.0 + stdlib_base_ln( lambda ) ) ) - ( 1.0 / ( 12.0*lambda ) ) - ( 1.0 / ( 24.0*l2 ) ) - ( 19.0 / ( 360.0*l2*lambda ) );
+	}
+	m = stdlib_base_floor( lambda );
+	wm1 = 0.0; // W - 1
+	sum = 0.0; // Σ w_k ln(w_k)
 
-		// Convergence: check relative size of term
-		if ( term < 1.0e-16 * sum ) {
+	// Terms above the mode, using w_{k+1} = w_k * λ/(k+1):
+	lw = 0.0;
+	for ( k = m; ; k++ ) {
+		lw += lnratio( lambda, k+1.0 );
+		w = stdlib_base_exp( lw );
+		wm1 += w;
+		sum += w * lw;
+		if ( w < TOL ) {
+			break;
+		}
+	}
+	// Terms below the mode, using w_{k-1} = w_k * k/λ:
+	lw = 0.0;
+	for ( k = m; k > 0.0; k-- ) {
+		lw += lnratio( k, lambda );
+		w = stdlib_base_exp( lw );
+		wm1 += w;
+		sum += w * lw;
+		if ( w < TOL ) {
 			break;
 		}
 	}
-	out += stdlib_base_exp( -lambda ) * sum;
-	return out;
+	return stdlib_base_log1p( wm1 ) - ( sum / ( 1.0+wm1 ) );
 }

lib/node_modules/@stdlib/stats/base/dists/poisson/entropy/test/fixtures/julia/REQUIRE

@@ -1,3 +1,3 @@
-Distributions 0.23.8
 julia 1.5
 JSON 0.21
+SpecialFunctions 2.4.0

lib/node_modules/@stdlib/stats/base/dists/poisson/entropy/test/fixtures/julia/large.json

@@ -0,0 +1 @@
+{"expected":[7.636170468247301,7.972537561030521,6.7610278260932,6.917374059516255,6.727728616304167,6.750677706556752,7.244541238871671,6.509536989419537,7.821177896107505,7.269680611911745,7.970309910988387,6.524486872450528,7.21123753063748,6.453430815347789,6.710260269147566,7.591407693115048,7.719570138214441,6.60409923648532,7.275309925886296,6.569444099784174,7.220662662609307,6.195700902330057,6.549083061504953,7.650034426380673,6.943061039654225,5.758624091480383,6.314255620964442,6.854093883650863,7.092309409089726,6.267512362927276,7.234898393348585,7.645304883181197,6.77246469767916,6.690957199656208,6.370033457060966,6.414296508435013,7.346668967419823,6.369298474087861,6.498332839194182,7.529607127859304,7.78693387327993,7.0160844635493245,6.928044277357777,7.139826340202468,6.062510702123153,6.886249261145232,7.1787683992812275,7.110396468179512,6.486337927437986,6.644103590157126,7.433638895070484,6.570855743395127,7.434170413341826,7.797867116769563,7.102138962481776,7.11346739600779,7.962627037349753,6.861612790136011,7.27952065202966,7.193231082968424,7.239885488716427,6.893343446896546,7.485664259242629,7.957547318248099,7.470112622815086,7.475431581084728,7.787995523408879,6.753984879486759,7.6407055038951635,6.089444968267598,6.024594986004413,6.359695233109921,7.13377710402349,7.710607191114309,6.798541224414055,7.564456723270363,7.13716031513703,7.397293948193246,7.121284650541771,6.727903509452129,6.663396008687282,6.975561358374909,7.545200793938018,5.746414776130449,6.33421494448371,7.499186674598618,7.764468119811527,5.97851871881829,7.0678479033698505,6.733128250255487,6.848576538876062,6.072626808270496,6.7933007812668835,7.575245446153155,7.819236694491278,6.661459830661217,6.581854920708985,6.638681032060693,6.631801879734862,7.259650282374132],"lambda":[251315.53326565857,492474.71298286394,43659.771794731394,59687.2581321705,40846.819797937875,42765.2973878073,114829.86335154864,26402.21600383054,363843.52720370906,120750.96191720419,490285.47158018046,27203.5512726292,107430.51508539419,23599.755402424704,39444.41313241093,229794.11735901757,296933.5844250128,31898.89445166201,122118.13198356815,29762.864363011704,109474.81275048848,14094.452562317174,28575.211472947598,258381.49418485272,62833.752321240754,5880.514053612631,17865.799193184277,52591.714641078666,84689.230158319,16271.296811890477,112636.51143423203,255948.96576958572,44669.93950410649,37950.64294699449,19974.229406768616,21823.08400408184,140851.5885728032,19944.889784149582,25817.170369697844,203076.49506414964,339758.77880560927,72714.35888680838,60974.69504154627,93132.39887869543,10798.470393349922,56085.04214822674,100675.8660019734,87808.86722918529,25205.19494818494,34555.946838829004,167611.1040471175,29847.01173730285,167789.37533616516,347269.9280158415,86370.61460762519,88349.83516636364,482809.45586904045,53388.55283910989,123150.88513890035,103630.47313789162,113765.58927913052,56886.46721094308,185990.77378900204,477929.216950017,180294.89658037858,182223.0946912322,340480.95478128234,43049.098297129065,253605.35026226818,11396.11168158752,10009.89674047982,19565.477074651004,92012.4297604441,291658.21154573106,47061.42413460863,217735.68359812972,92637.13470783619,155859.75871650188,89741.99335667417,40861.109896531,35915.33318731494,67053.63737486127,209509.6997880244,5738.662975894696,18593.396811852715,191089.49451916435,324830.7921852627,9128.70049245127,80645.70336010086,41290.324199144525,52014.57337985769,11019.169049149645,46570.7562549823,222484.91377096222,362433.67913489684,35776.52580414858,30510.86832457014,34183.21022702833,33716.12998381677,118352.75659560131]}