fix: evaluate both tails without cancellation in stats/base/dists/cauchy/logcdf

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/lib/factory.js

@@ -24,6 +24,7 @@ var constantFunction = require( '@stdlib/utils/constant-function' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var atan2 = require( '@stdlib/math/base/special/atan2' );
 var ln = require( '@stdlib/math/base/special/ln' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
 
 
 // VARIABLES //
@@ -74,6 +75,14 @@ function factory( x0, gamma ) {
 		if ( isnan( x ) ) {
 			return NaN;
 		}
+		// Below `z = -1`, `0.5 + atan(z)/pi` cancels more and more as `z` goes to `-infinity`; `atan2( gamma, x0-x )/pi` is the same quantity evaluated without the subtraction. Nearer the median the sum does not cancel and is the more accurate of the two...
+		if ( x-x0 < -gamma ) {
+			return ln( ONE_OVER_PI * atan2( gamma, x0-x ) );
+		}
+		// Above `z = 1`, the result is `ln` of a number approaching `1`, so evaluate `ln(1 - atan2(gamma, x-x0)/pi)` with `log1p` instead...
+		if ( x-x0 > gamma ) {
+			return log1p( -ONE_OVER_PI * atan2( gamma, x-x0 ) );
+		}
 		return ln( ( ONE_OVER_PI * atan2( x-x0, gamma ) ) + 0.5 );
 	}
 }

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/lib/main.js

@@ -23,6 +23,7 @@
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var atan2 = require( '@stdlib/math/base/special/atan2' );
 var ln = require( '@stdlib/math/base/special/ln' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
 
 
 // VARIABLES //
@@ -73,6 +74,14 @@ function logcdf( x, x0, gamma ) {
 	) {
 		return NaN;
 	}
+	// Below `z = -1`, `0.5 + atan(z)/pi` cancels more and more as `z` goes to `-infinity`; `atan2( gamma, x0-x )/pi` is the same quantity evaluated without the subtraction. Nearer the median the sum does not cancel and is the more accurate of the two...
+	if ( x-x0 < -gamma ) {
+		return ln( ONE_OVER_PI * atan2( gamma, x0-x ) );
+	}
+	// Above `z = 1`, the result is `ln` of a number approaching `1`, so evaluate `ln(1 - atan2(gamma, x-x0)/pi)` with `log1p` instead...
+	if ( x-x0 > gamma ) {
+		return log1p( -ONE_OVER_PI * atan2( gamma, x-x0 ) );
+	}
 	return ln( ( ONE_OVER_PI * atan2( x-x0, gamma ) ) + 0.5 );
 }
 

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/manifest.json

@@ -43,6 +43,7 @@
         "@stdlib/math/base/assert/is-infinite",
         "@stdlib/math/base/special/atan2",
         "@stdlib/math/base/special/ln",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ninf"
       ]
     },
@@ -63,6 +64,7 @@
         "@stdlib/math/base/assert/is-infinite",
         "@stdlib/math/base/special/atan2",
         "@stdlib/math/base/special/ln",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ninf"
       ]
     },
@@ -83,6 +85,7 @@
         "@stdlib/math/base/assert/is-infinite",
         "@stdlib/math/base/special/atan2",
         "@stdlib/math/base/special/ln",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ninf"
       ]
     }

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/src/main.c

@@ -21,6 +21,7 @@
 #include "stdlib/math/base/assert/is_infinite.h"
 #include "stdlib/math/base/special/atan2.h"
 #include "stdlib/math/base/special/ln.h"
+#include "stdlib/math/base/special/log1p.h"
 #include "stdlib/constants/float64/ninf.h"
 
 static const double ONE_OVER_PI = 0.3183098861837907;
@@ -49,5 +50,13 @@ double stdlib_base_dists_cauchy_logcdf( const double x, const double x0, const d
 	if ( stdlib_base_is_infinite( x ) ) {
 		return ( x < 0.0 ) ? STDLIB_CONSTANT_FLOAT64_NINF : 0.0;
 	}
+	// Below `z = -1`, `0.5 + atan(z)/pi` cancels more and more as `z` goes to `-infinity`; `atan2( gamma, x0-x )/pi` is the same quantity evaluated without the subtraction. Nearer the median the sum does not cancel and is the more accurate of the two...
+	if ( x-x0 < -gamma ) {
+		return stdlib_base_ln( ONE_OVER_PI * stdlib_base_atan2( gamma, x0-x ) );
+	}
+	// Above `z = 1`, the result is `ln` of a number approaching `1`, so evaluate `ln(1 - atan2(gamma, x-x0)/pi)` with `log1p` instead...
+	if ( x-x0 > gamma ) {
+		return stdlib_base_log1p( -ONE_OVER_PI * stdlib_base_atan2( gamma, x-x0 ) );
+	}
 	return stdlib_base_ln( ( ONE_OVER_PI * stdlib_base_atan2( x-x0, gamma ) ) + 0.5 );
 }

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/test/fixtures/julia/runner.jl

@@ -16,7 +16,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import Distributions: logcdf, Cauchy
 import JSON
 
 """
@@ -41,9 +40,9 @@ julia> gen( x, x0, gamma, \"data.json\" );
 ```
 """
 function gen( x, x0, gamma, name )
-	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logcdf( Cauchy( x0[i], gamma[i] ), x[i] );
+	# The log CDF is evaluated in extended precision, as `0.5 + atan(z)/pi` cancels in double precision left of the median (Distributions.jl's `logcdf` evaluates it in `Float64` and is off by hundreds of ULP there)...
+	z = setprecision( BigFloat, 2048 ) do
+		Float64.( log.( 0.5 .+ ( atan.( ( BigFloat.( x ) .- BigFloat.( x0 ) ) ./ BigFloat.( gamma ) ) ./ BigFloat( pi ) ) ) );
 	end
 
 	# Store data to be written to file as a collection:
@@ -87,3 +86,10 @@ x = rand( 1000 ) .* 2.0;
 x0 = rand( 1000 ) .* 1.0;
 gamma = rand( 1000 ) .* 50.0;
 gen( x, x0, gamma, "large_gamma.json" );
+
+# Far left tail, many orders of magnitude below the median:
+z = .-exp10.( ( rand( 1000 ) .* 290.0 ) .+ 6.0 );
+x0 = ( rand( 1000 ) .* 20.0 ) .- 10.0;
+gamma = ( rand( 1000 ) .* 10.0 ) .+ 0.1;
+x = ( z .* gamma ) .+ x0;
+gen( x, x0, gamma, "left_tail.json" );

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/test/test.factory.js

@@ -32,6 +32,7 @@ var factory = require( './../lib/factory.js' );
 
 var largeGamma = require( './fixtures/julia/large_gamma.json' );
 var negativeMedian = require( './fixtures/julia/negative_median.json' );
+var leftTail = require( './fixtures/julia/left_tail.json' );
 var positiveMedian = require( './fixtures/julia/positive_median.json' );
 
 
@@ -197,3 +198,24 @@ tape( 'the created function evaluates the logcdf for `x` given `x0` and `gamma`
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the logcdf for `x` given `x0` and `gamma` far into the left tail', function test( t ) {
+	var expected;
+	var logcdf;
+	var gamma;
+	var x0;
+	var i;
+	var x;
+	var y;
+
+	expected = leftTail.expected;
+	x = leftTail.x;
+	x0 = leftTail.x0;
+	gamma = leftTail.gamma;
+	for ( i = 0; i < x.length; i++ ) {
+		logcdf = factory( x0[i], gamma[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/test/test.logcdf.js

@@ -32,6 +32,7 @@ var logcdf = require( './../lib' );
 
 var largeGamma = require( './fixtures/julia/large_gamma.json' );
 var negativeMedian = require( './fixtures/julia/negative_median.json' );
+var leftTail = require( './fixtures/julia/left_tail.json' );
 var positiveMedian = require( './fixtures/julia/positive_median.json' );
 
 
@@ -148,3 +149,22 @@ tape( 'the function evaluates the logcdf for `x` given `x0` and `gamma` (`x0 > 0
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logcdf for `x` given `x0` and `gamma` far into the left tail', function test( t ) {
+	var expected;
+	var gamma;
+	var x0;
+	var x;
+	var y;
+	var i;
+
+	expected = leftTail.expected;
+	x = leftTail.x;
+	x0 = leftTail.x0;
+	gamma = leftTail.gamma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], x0[i], gamma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});

lib/node_modules/@stdlib/stats/base/dists/cauchy/logcdf/test/test.native.js

@@ -33,6 +33,7 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 
 var largeGamma = require( './fixtures/julia/large_gamma.json' );
 var negativeMedian = require( './fixtures/julia/negative_median.json' );
+var leftTail = require( './fixtures/julia/left_tail.json' );
 var positiveMedian = require( './fixtures/julia/positive_median.json' );
 
 
@@ -157,3 +158,22 @@ tape( 'the function evaluates the logcdf for `x` given `x0` and `gamma` (`x0 > 0
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logcdf for `x` given `x0` and `gamma` far into the left tail', opts, function test( t ) {
+	var expected;
+	var gamma;
+	var x0;
+	var x;
+	var y;
+	var i;
+
+	expected = leftTail.expected;
+	x = leftTail.x;
+	x0 = leftTail.x0;
+	gamma = leftTail.gamma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], x0[i], gamma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});