fix: use expm1 in stats/base/dists/lognormal/skewness

`exp(sigma^2) - 1` cancels for small `sigma` and is `0` for
`sigma^2 < 2^-53`. The fixtures are now generated from the textbook
formula in extended precision.

Author: Abhist17

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/lib/main.js

@@ -23,6 +23,7 @@
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 var exp = require( '@stdlib/math/base/special/exp' );
+var expm1 = require( '@stdlib/math/base/special/expm1' );
 
 
 // MAIN //
@@ -55,16 +56,16 @@ var exp = require( '@stdlib/math/base/special/exp' );
 * // returns NaN
 */
 function skewness( mu, sigma ) {
-	var es2;
+	var s2;
 	if (
 		isnan( mu ) ||
 		isnan( sigma ) ||
 		sigma <= 0.0
 	) {
 		return NaN;
 	}
-	es2 = exp( sigma*sigma );
-	return ( es2 + 2.0 ) * sqrt( es2 - 1.0 );
+	s2 = sigma * sigma;
+	return ( exp( s2 ) + 2.0 ) * sqrt( expm1( s2 ) );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/manifest.json

@@ -40,6 +40,7 @@
       "dependencies": [
         "@stdlib/math/base/napi/binary",
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/math/base/special/exp",
         "@stdlib/math/base/special/sqrt"
       ]
@@ -57,6 +58,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/constants/float64/eps",
         "@stdlib/math/base/special/exp",
         "@stdlib/math/base/special/sqrt"
@@ -75,6 +77,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/math/base/special/exp",
         "@stdlib/math/base/special/sqrt"
       ]

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/src/main.c

@@ -20,6 +20,7 @@
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/exp.h"
 #include "stdlib/math/base/special/sqrt.h"
+#include "stdlib/math/base/special/expm1.h"
 
 /**
 * Returns the skewness for a lognormal distribution with location `mu` and scale `sigma`.
@@ -33,14 +34,14 @@
 * // returns ~6.185
 */
 double stdlib_base_dists_lognormal_skewness( const double mu, const double sigma ) {
-	double es2;
+	double s2;
 	if (
 		stdlib_base_is_nan( mu ) ||
 		stdlib_base_is_nan( sigma ) ||
 		sigma <= 0.0
 	) {
 		return 0.0/0.0; // NaN
 	}
-	es2 = stdlib_base_exp( sigma * sigma );
-	return ( es2 + 2.0 ) * stdlib_base_sqrt( es2 - 1.0 );
+	s2 = sigma * sigma;
+	return ( stdlib_base_exp( s2 ) + 2.0 ) * stdlib_base_sqrt( stdlib_base_expm1( s2 ) );
 }

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/test/fixtures/julia/REQUIRE

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

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/test/fixtures/julia/data.json

@@ -1 +1 @@
-{"sigma":[0.6370150969317179,4.553781304799195,1.6897136584317418,0.9193839303313622,3.9822924852920294,4.474403219035906,0.592686297969246,2.7416296765503634,1.7981740492043807,4.4431868848393155,4.148532554898559,3.2913643963196506,4.79141187182625,3.5363073062288897,3.9947975668429283,3.4607327154768663,4.911303012475389,1.153344635169462,2.041732268634454,0.6249634293608552,0.8796871819486618,4.505305536119697,3.581403802763672,1.6901073115863807,4.028177476739651,4.529164600821525,0.6750797988399415,3.4759887479561002,4.087289491837983,4.268967326830709,2.1173795928341654,2.8521120657353394,2.8596338023527537,1.2474558177551986,1.0498505248902934,0.6152484292375637,4.178197027660325,3.2032444242379343,1.9659190803865134,3.7352788529470726,1.6063290781847206,3.134552953557729,2.0466656611158704,4.113090594504293,1.2489460035017719,3.6707704860925916,1.1651214830263745,0.08478662380910329,0.8086809083977753,2.201490352168025,0.07716735510295947,0.7924154094400659,2.4284514514483635,2.9665552495224334,2.9598070594981496,0.2941259897829318,4.582044768029335,3.520641833674695,4.629034796920037,3.418691202355343,4.191087151781913,1.2951657415939288,3.4413387683055277,0.42232175543560957,4.858019378997895,4.929173838427397,1.7758826298953667,2.804153584059721,1.8098649274171819,4.467112160686821,4.860331022607472,1.5303846176156688,3.3970668478318977,0.42105512788462307,2.176553793746152,0.6214910206976132,0.4784317167225749,4.113666856385336,2.594848023703337,0.8649192973711628,4.72548620803634,1.8773557485182912,3.374710792154074,3.515031123243534,3.745189397137043,2.4079901489026243,1.0837585143927597,0.03716708074650854,0.1739616625390661,1.2002513470872567,4.676855220432165,4.150460400523116,3.4830478130236586,3.94324591178333,0.649852688699607,2.4636156066636783,0.9665074266851637,4.3911827694309915,1.0450824557334504,2.7962268422463965],"expected":[2.4764161279229544,3.2277330556233926e13,78.41411831396505,4.989354047981641,2.142834364543264e10,1.1016621049496807e13,2.2193079963557567,78874.65946176853,135.0878068569964,7.256106777708033e12,1.6275040213376935e11,1.1405987408569362e7,9.027272620818556e14,1.4014800568183956e8,2.488700399000577e10,6.34026252147622e7,5.168432885591407e15,9.643355887172216,531.4906307411909,2.40407992534158,4.504816360399443,1.670396050098178e13,2.2682374740483576e8,78.56317312785403,3.719056030555334e10,2.3080240409090246e13,2.71813130362016,7.430993016207147e7,7.63736186961147e10,7.446075027733396e11,846.9190491820564,199237.8287129027,212496.20333464557,13.036437228304512,7.10525978404892,2.3471302422409144,2.3574007130816452e11,4.834144475822791e6,339.59014281690975,1.2277676910782735e9,53.09702552780883,2.516084322903095e6,547.5738920163817,1.0489979785678168e11,13.10144061990895,5.996376126423744e8,10.000976003072989,0.25543049930664,3.7693917873603113,1453.1444090846953,0.23230870448483143,3.6208617761232285,6975.558795448786,540852.0360112038,509365.90523380094,0.9289733039149168,4.754517299567284e13,1.1873261508117507e8,9.100583373659514e13,4.108636626068362e7,2.7714794721910205e11,15.337006309945094,5.1868880449695244e7,1.4118898211118003,2.367304663678963e15,6.728523301425043e15,120.39242273683148,132653.07873580247,143.59838417536739,9.99030848476256e12,2.4484375321981055e15,38.03203775899423,3.2937070724640626e7,1.4067062253649858,1235.2773967088087,2.38358766882956,1.651931570012242,1.0564841463813002e11,24382.897018109004,4.339018056351443,3.522385827430928e14,206.23277592021074,2.624605841897949e7,1.1190626457264058e8,1.3721801267928703e9,6015.848875097824,7.8315707624239295,0.1115911555821832,0.5312441081435969,11.172832080864644,1.7740602176717203e14,1.6670345527408813e11,7.999231946476696e7,1.3470699364061668e10,2.5556175499307954,9023.247830428249,5.649391968525006,3.6426377005462446e12,7.010075787095745,124110.74753325092],"mu":[-3.923554117872933,-2.672795436485885,-2.2122444072050516,-3.568408402715862,-2.834694003325013,-2.725331535073217,-2.696720096474188,-3.013878350382553,-2.831955545999753,-2.7466776644620574,-3.436034558072513,-3.4730286110609487,-3.0895963224853413,-2.821182924768898,-2.536832645622806,-2.726180231911774,-2.296595097554329,-3.8106401272740342,-2.8097696911094654,-2.811903806383317,-2.9102652982258954,-3.4702221056880895,-2.859708351105935,-2.94040597562659,-3.094945663865096,-3.921213930154313,-2.5428422503522072,-2.67379050585622,-3.3236780857123462,-3.296683356727204,-2.602938572351321,-2.0349645805917076,-2.2889456864786712,-3.9223062952962016,-3.13294976267802,-2.233985184170268,-3.569934997641718,-3.437432478800813,-3.0752775803431347,-3.153341799739629,-3.652085368482236,-3.17750789780354,-3.996905481342681,-2.3733415456074605,-3.9783320492367134,-2.9336113423072323,-3.5202708685513393,-2.0313600363807978,-3.1253966454638253,-3.1233368787404405,-3.4910650968663948,-3.2686712579252526,-3.468541520703441,-2.4037422195116056,-3.2912018628252024,-3.67769655578248,-3.9527043446738523,-2.635251157925184,-3.210084574521885,-3.273915042622466,-3.513348213929724,-3.7677657350609195,-2.046113084182904,-2.3043046798878972,-3.6946389627265965,-3.6412390555768788,-2.835260465939693,-3.788745038281978,-3.662923358238609,-2.049567207453098,-2.470588008940093,-2.502434172004819,-3.932048732697886,-2.616659622822441,-3.4497679546515565,-2.479342101997467,-3.264417494325702,-2.4990647808646704,-2.369693967221661,-2.2038070165521706,-2.191769371884672,-2.9470622028168543,-3.156591091420007,-3.484403017140843,-2.7777362689854512,-3.9460852955152887,-2.2544703906931467,-2.2751546746043183,-2.2433126347738788,-3.7413244558148167,-2.1145337639222843,-2.288208248069057,-2.6188399298463394,-2.8250164432507865,-2.2651687786087495,-3.5152560762188796,-2.922009810980355,-3.2420935888015845,-3.7921327960695006,-3.0477125944647567]}
+{"expected":[3.586985256469765e10,373.779791260857,59.50466708200734,3.4415380787977266e9,1.6386381184240403e7,1.0172997381180447e14,5.001805607433798,6.631036737328351e6,9352.628070998855,4.21337594957755,2.798781165421248e12,4922.33120787807,1.5228283842963958e9,2.674299871046978e6,122.13015174419219,1035.9189577933523,130199.87902695996,116.3227861915989,1.708643767226305,450765.0469074596,2.2312711672225504e9,2.7424631486785946,6.918366796398472,3.688357928457058,495372.6978925154,420.34920723270056,1.5177564931265095e6,5.219141379177634e7,0.26650041442215056,1.8846711306881598e15,2.41425814246967e12,1.8622265826813234,7.891027412326184e12,1.7307861623596563,4.249941382202212,2.530032197743837e14,5.460889545110747e9,5.203582270365993e9,1.9574014925502192e6,52803.08181343744,65.96099592541141,1.3902972400189466e6,1.2394530761239788e15,11.82688914346538,456271.4450497873,4.975214320859048,134.76806728901326,7.231802463368179e12,1.3625426264019296,101.45934825209214,0.8810317326118378,14.876345798351975,589.316578179464,84.49936168110384,4.8138671270919913e8,1.6789227630633255e6,2901.6573649244165,299.0562351570923,0.18943479903118687,3.3782410780115786e12,4.5182608698603824e7,14799.38362495365,153.03499783617286,1.4787171213503811,0.2495303450528318,6.8265225902741445e13,2.155542393340309,0.02702997648974367,41.15732905185768,1.7694807396843192e11,23268.73663699636,12282.409756926612,1.4124481980847951e6,1.4488186906704037e10,23870.95098171939,3.425631166668912e11,2665.717083508431,738.8305672334942,0.1005071731034523,1.5348036367493492e15,1.9085528172915957e6,7.10114857582941e11,7.569194152072299e7,5.482172258624223,2.4702398461938917,22.63782079896971,2141.5275575316045,1.49508722630045e8,1.6936777959744817e12,4.198671339857217,47.52472960008406,1.853523042772858e10,5.427293926992595,3.8114737221362585e8,308.8406037116571,4.0175138700804415e15,521.1408363576177,1495.6003784888856,1.738523905776192,5.123230676704532],"mu":[-3.9919212771653934,-2.4675203577738043,-2.3133232397389794,-2.3134108121578683,-3.004026065464043,-2.05837356336741,-2.49812318996603,-2.5306416807810663,-3.447118012337229,-3.418843874046597,-3.1872693466060618,-2.0643024974798223,-3.8012904103193064,-2.4810427668105675,-3.155539828518258,-3.657382153556609,-2.343117512369581,-3.8393558351795316,-2.3845400796305842,-3.805609275297729,-3.2922960276751367,-3.1166519321935797,-2.848917217401585,-2.1723970959391603,-2.007586713345928,-2.8251236917531637,-3.906335860990563,-2.3434668343405507,-3.1926903406977436,-2.3500511833660123,-3.8603944782444453,-2.873041601023419,-2.0905777542698782,-2.057732171465342,-2.7206857427268645,-3.5525749938750093,-2.3131966765397918,-3.259483195538056,-3.8245388126557156,-2.5169496747650673,-2.038229658358553,-3.197245789715201,-2.368380378618896,-3.6546116786562046,-2.841073658719451,-2.5620403262584777,-3.9568635162316195,-3.4985299102384637,-3.7246889282351567,-3.359811144466793,-3.7261182754794495,-2.0314754507295776,-3.968106221129471,-2.8809599067250575,-3.3998650319308092,-3.687300349861946,-3.0386362279036954,-3.6956158078032235,-2.2805054072167215,-2.687022909305322,-2.045240004984451,-2.3672307304768267,-3.404157299923671,-3.8109000246517146,-3.5529467017199456,-2.106985916533536,-3.7092708413996203,-3.8354861620375713,-3.869017239043697,-3.299423136973817,-3.858244815311817,-3.928695208374391,-3.119659522361263,-3.030661288917398,-3.7256252627724735,-3.467099811925748,-2.505505671571521,-2.1926103786909747,-2.7753365145332727,-3.210763768421037,-3.5865867432228002,-2.741752066302074,-2.78318037344663,-2.5557937447205887,-3.941135506328295,-2.331113851823261,-3.218108562714878,-2.9207219197606387,-3.2510886399866488,-3.6659424606050037,-3.007210629822632,-3.8247620482207,-3.456847018723839,-3.5326032057931442,-2.0372619720285288,-2.164365042198029,-2.1704630204476736,-3.733596230484008,-2.2187844516711834,-2.2687255930195844],"sigma":[4.0251842922862195,1.9824315289274075,1.6312485031118258,3.8261541758698088,3.3278559545850865,4.637049572204105,0.9203421311078042,3.235967575342298,2.468472858580671,0.8532791001618828,4.37113296775739,2.379966779875575,3.754449279507997,3.141031711713186,1.7786748659755007,2.149180680030378,2.8019329371854913,1.7691605806017678,0.4909769796864488,2.946008191522815,3.788212983425983,0.6787332811533775,1.0404105101614318,0.799897495837899,2.9566677315027756,2.0024325696405025,3.0803279304375364,3.441939163911001,0.08842847259278197,4.84234996215769,4.35984809828074,0.5236428146453997,4.44947543391722,0.4958033771418563,0.8567081878239957,4.702086297987385,3.8661675295720803,3.862004024162916,3.107735780857265,2.692366923212516,1.6533749042934143,3.070820488868135,4.813415046129358,1.217885663029098,2.9473819022480763,0.9182922951790956,1.797718747356356,4.442935172487339,0.41016779425910754,1.7421229668881004,0.2802921066706343,1.2863644796016904,2.058766822778967,1.7051214804719224,3.650769694141344,3.0912302567264103,2.304442440400722,1.9437913680639918,0.06299878047852536,4.385459025171693,3.4279450582099864,2.529802948273697,1.8219480697776158,0.4384408268751472,0.08284393515215038,4.608283962197117,0.5810243784152835,0.009009565537957953,1.5487821632453476,4.155247454442662,2.5888255648782943,2.505075957360488,3.072535968640209,3.9493965113539984,2.592117454861059,4.207906980061199,2.2920801006895015,2.095523915181539,0.033480485987031905,4.828193509220079,3.105023707888756,4.265262207359698,3.477755394608026,0.9552330710299228,0.6359984375155942,1.4019208131894012,2.2598339046243288,3.542396549374258,4.332660728345871,0.8518903120056542,1.5816044156048348,3.9701334059310374,0.9514402859032423,3.6293887280197143,1.9494231682361685,4.894176108116714,2.038470664032757,2.2058797435512867,0.49748056866430024,0.929537003331069]}

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/test/fixtures/julia/runner.jl

@@ -16,14 +16,15 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import Distributions: skewness, LogNormal
 import JSON
 
 """
 	gen( mu, sigma, name )
 
 Generate fixture data and write to file.
 
+The skewness is evaluated in extended precision from the textbook formula `(exp(sigma^2) + 2) * sqrt(exp(sigma^2) - 1)`, as subtracting `1` from `exp(sigma^2)` in double precision cancels for small `sigma`.
+
 # Arguments
 
 * `mu`: location parameter
@@ -39,9 +40,9 @@ julia> gen( mu, sigma, "data.json" );
 ```
 """
 function gen( mu, sigma, name )
-	z = Array{Float64}( undef, length(mu) );
-	for i in eachindex(mu)
-		z[ i ] = skewness( LogNormal( mu[i], sigma[i] ) );
+	z = setprecision( BigFloat, 2048 ) do
+		es2 = exp.( BigFloat.( sigma ) .^ 2 );
+		Float64.( ( es2 .+ 2 ) .* sqrt.( es2 .- 1 ) )
 	end
 
 	# Store data to be written to file as a collection:
@@ -71,3 +72,8 @@ dir = dirname( file );
 mu = rand( 100 ) .* 2.0 .- 4.0;
 sigma = rand( 100 ) .* 5.0;
 gen( mu, sigma, "data.json" );
+
+# Small `sigma` (above `1e-154`, where `sigma^2` underflows):
+mu = ( rand( 1000 ) .* 4.0 ) .- 2.0;
+sigma = exp10.( range( -150.0, stop = -1.0, length = 1000 ) );
+gen( mu, sigma, "small_sigma.json" );

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/test/fixtures/julia/small_sigma.json

@@ -31,6 +31,7 @@ var skewness = require( './../lib' );
 // FIXTURES //
 
 var data = require( './fixtures/julia/data.json' );
+var smallSigma = require( './fixtures/julia/small_sigma.json' );
 
 
 // TESTS //
@@ -86,7 +87,26 @@ tape( 'the function returns the skewness of a lognormal distribution', function
 	for ( i = 0; i < mu.length; i++ ) {
 		y = skewness( mu[i], sigma[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 22 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function returns the skewness of a lognormal distribution with a small `sigma`', function test( t ) {
+	var expected;
+	var sigma;
+	var mu;
+	var y;
+	var i;
+
+	expected = smallSigma.expected;
+	mu = smallSigma.mu;
+	sigma = smallSigma.sigma;
+	for ( i = 0; i < mu.length; i++ ) {
+		y = skewness( mu[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 		}
 	}
 	t.end();

lib/node_modules/@stdlib/stats/base/dists/lognormal/skewness/test/test.js

@@ -32,6 +32,7 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 // FIXTURES //
 
 var data = require( './fixtures/julia/data.json' );
+var smallSigma = require( './fixtures/julia/small_sigma.json' );
 
 
 // VARIABLES //
@@ -98,7 +99,26 @@ tape( 'the function returns the skewness of a lognormal distribution', opts, fun
 	for ( i = 0; i < mu.length; i++ ) {
 		y = skewness( mu[i], sigma[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 22 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function returns the skewness of a lognormal distribution with a small `sigma`', opts, function test( t ) {
+	var expected;
+	var sigma;
+	var mu;
+	var y;
+	var i;
+
+	expected = smallSigma.expected;
+	mu = smallSigma.mu;
+	sigma = smallSigma.sigma;
+	for ( i = 0; i < mu.length; i++ ) {
+		y = skewness( mu[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 		}
 	}
 	t.end();