fix: avoid cancellation for large b in stats/base/dists/kumaraswamy/median

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/median/lib/main.js

@@ -21,7 +21,9 @@
 // MODULES //
 
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
+var expm1 = require( '@stdlib/math/base/special/expm1' );
 var pow = require( '@stdlib/math/base/special/pow' );
+var LN2 = require( '@stdlib/constants/float64/ln-two' );
 
 
 // MAIN //
@@ -70,7 +72,8 @@ function median( a, b ) {
 	) {
 		return NaN;
 	}
-	return pow( 1.0 - pow( 2.0, -1.0/b ), 1.0/a );
+	// `1 - 2^(-1/b)` cancels for large `b`; `-expm1(-ln(2)/b)` is the same value without the subtraction:
+	return pow( -expm1( -LN2/b ), 1.0/a );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/median/manifest.json

@@ -40,7 +40,9 @@
       "dependencies": [
         "@stdlib/math/base/napi/binary",
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/pow"
+        "@stdlib/math/base/special/expm1",
+        "@stdlib/math/base/special/pow",
+        "@stdlib/constants/float64/ln-two"
       ]
     },
     {
@@ -56,7 +58,9 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/math/base/special/pow",
+        "@stdlib/constants/float64/ln-two",
         "@stdlib/constants/float64/eps"
       ]
     },
@@ -73,7 +77,9 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/pow"
+        "@stdlib/math/base/special/expm1",
+        "@stdlib/math/base/special/pow",
+        "@stdlib/constants/float64/ln-two"
       ]
     }
   ]

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/median/src/main.c

@@ -18,7 +18,9 @@
 
 #include "stdlib/stats/base/dists/kumaraswamy/median.h"
 #include "stdlib/math/base/assert/is_nan.h"
+#include "stdlib/math/base/special/expm1.h"
 #include "stdlib/math/base/special/pow.h"
+#include "stdlib/constants/float64/ln_two.h"
 
 /**
 * Returns the median of a Kumaraswamy's double bounded distribution.
@@ -39,5 +41,6 @@ double stdlib_base_dists_kumaraswamy_median( const double a, const double b ) {
 	) {
 		return 0.0/0.0; // NaN
 	}
-	return stdlib_base_pow( 1.0 - stdlib_base_pow( 2.0, -1.0/b ), 1.0/a );
+	// `1 - 2^(-1/b)` cancels for large `b`; `-expm1(-ln(2)/b)` is the same value without the subtraction:
+	return stdlib_base_pow( -stdlib_base_expm1( -STDLIB_CONSTANT_FLOAT64_LN2 / b ), 1.0/a );
 }

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/median/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/kumaraswamy/median/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: median, Kumaraswamy
 import JSON
 
 """
 	gen( a, b, name )
 
 Generate fixture data and write to file.
 
+The median is evaluated in extended precision from the textbook formula `(1 - 2^(-1/b))^(1/a)`, as forming `1 - 2^(-1/b)` in double precision cancels for large `b`.
+
 # Arguments
 
 * `a`: first shape parameter
@@ -39,9 +40,8 @@ julia> gen( a, b, \"data.json\" );
 ```
 """
 function gen( a, b, name )
-	z = Array{Float64}( undef, length(a) );
-	for i in eachindex(a)
-		z[ i ] = median( Kumaraswamy( a[ i ], b[ i ] ) );
+	z = setprecision( BigFloat, 2048 ) do
+		Float64.( ( 1 .- BigFloat( 2 ) .^ ( -1 ./ BigFloat.( b ) ) ) .^ ( 1 ./ BigFloat.( a ) ) )
 	end
 
 	# Store data to be written to file as a collection:
@@ -71,3 +71,8 @@ dir = dirname( file );
 a = rand( 1000 );
 b = rand( 1000 );
 gen( a, b, "data.json" );
+
+# Large `b`:
+a = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+b = exp10.( range( 0.0, stop = 15.0, length = 1000 ) );
+gen( a, b, "large_b.json" );

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/median/test/test.js

@@ -31,6 +31,7 @@ var median = require( './../lib' );
 // FIXTURES //
 
 var data = require( './fixtures/julia/data.json' );
+var largeB = require( './fixtures/julia/large_b.json' );
 
 
 // TESTS //
@@ -114,7 +115,24 @@ tape( 'the function returns the median of a Kumaraswamy distribution', function
 	b = data.b;
 	for ( i = 0; i < expected.length; i++ ) {
 		y = median( a[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 9 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 442 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function returns the median of a Kumaraswamy distribution with a large `b`', function test( t ) {
+	var expected;
+	var a;
+	var b;
+	var i;
+	var y;
+
+	expected = largeB.expected;
+	a = largeB.a;
+	b = largeB.b;
+	for ( i = 0; i < expected.length; i++ ) {
+		y = median( a[i], b[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 24 ), true, 'returns expected value' );
 	}
 	t.end();
 });

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/median/test/test.native.js

@@ -32,6 +32,7 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 // FIXTURES //
 
 var data = require( './fixtures/julia/data.json' );
+var largeB = require( './fixtures/julia/large_b.json' );
 
 
 // VARIABLES //
@@ -126,7 +127,24 @@ tape( 'the function returns the median of a Kumaraswamy distribution', opts, fun
 	b = data.b;
 	for ( i = 0; i < expected.length; i++ ) {
 		y = median( a[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 9 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 442 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function returns the median of a Kumaraswamy distribution with a large `b`', opts, function test( t ) {
+	var expected;
+	var a;
+	var b;
+	var i;
+	var y;
+
+	expected = largeB.expected;
+	a = largeB.a;
+	b = largeB.b;
+	for ( i = 0; i < expected.length; i++ ) {
+		y = median( a[i], b[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 24 ), true, 'returns expected value' );
 	}
 	t.end();
 });