Skip to content

Commit 8251feb

Browse files
Abhist17Abhist17
committed
fix: avoid overflow and underflow for scales far from 1 in stats/base/dists/rayleigh/logpdf
For `sigma` outside about [1.5e-154, 1.3e154], `sigma^2` overflows or underflows, and so does `x^2` for large `x`, so the log-PDF came out NaN or a spurious -Infinity where it is finite: `logpdf( 1e200, 1e200 )` is about -461.017 and returned NaN. Keep the existing expression wherever it is safe — it is unchanged, bit for bit, on every point of the existing fixtures — and otherwise evaluate via the ratio `x/sigma`, taking `ln( x/sigma^2 )` directly while that ratio is normal so the two logarithms don't cancel. New `tiny_scale` and `huge_scale` fixtures cover both sides; develop fails all 2,000 of their points.
1 parent 5fc7eb7 commit 8251feb

10 files changed

Lines changed: 240 additions & 10 deletions

File tree

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/lib/factory.js‎

Lines changed: 25 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ var pow = require( '@stdlib/math/base/special/pow' );
2727
var ln = require( '@stdlib/math/base/special/ln' );
2828
var PINF = require( '@stdlib/constants/float64/pinf' );
2929
var NINF = require( '@stdlib/constants/float64/ninf' );
30+
var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
3031

3132

3233
// MAIN //
@@ -46,7 +47,9 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
4647
* // returns ~-0.936
4748
*/
4849
function factory( sigma ) {
50+
var normal;
4951
var s2i;
52+
var lns;
5053
var s2;
5154
if ( isnan( sigma ) || sigma < 0.0 ) {
5255
return constantFunction( NaN );
@@ -56,6 +59,8 @@ function factory( sigma ) {
5659
}
5760
s2 = pow( sigma, 2.0 );
5861
s2i = 1.0 / s2;
62+
normal = ( s2 >= FLOAT64_SMALLEST_NORMAL && s2 < PINF );
63+
lns = ln( sigma );
5964
return logpdf;
6065

6166
/**
@@ -70,13 +75,32 @@ function factory( sigma ) {
7075
* // returns <number>
7176
*/
7277
function logpdf( x ) {
78+
var lq;
79+
var x2;
80+
var r;
81+
var q;
7382
if ( isnan( x ) ) {
7483
return NaN;
7584
}
7685
if ( x < 0.0 || x === PINF ) {
7786
return NINF;
7887
}
79-
return ln( s2i * x ) - (pow( x, 2.0 ) / ( 2.0 * s2 ));
88+
x2 = pow( x, 2.0 );
89+
q = s2i * x;
90+
if ( normal && x2 < PINF && q >= FLOAT64_SMALLEST_NORMAL ) {
91+
return ln( q ) - (x2 / ( 2.0 * s2 ));
92+
}
93+
// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
94+
r = x / sigma;
95+
q = r / sigma;
96+
if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {
97+
lq = ln( q );
98+
} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {
99+
lq = ln( r ) - lns;
100+
} else {
101+
lq = ( ln( x ) - lns ) - lns;
102+
}
103+
return lq - ( 0.5 * r * r );
80104
}
81105
}
82106

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/lib/main.js‎

Lines changed: 26 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ var pow = require( '@stdlib/math/base/special/pow' );
2525
var ln = require( '@stdlib/math/base/special/ln' );
2626
var PINF = require( '@stdlib/constants/float64/pinf' );
2727
var NINF = require( '@stdlib/constants/float64/ninf' );
28+
var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
2829

2930

3031
// MAIN //
@@ -63,7 +64,11 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
6364
*/
6465
function logpdf( x, sigma ) {
6566
var s2i;
67+
var lq;
6668
var s2;
69+
var x2;
70+
var r;
71+
var q;
6772
if (
6873
isnan( x ) ||
6974
isnan( sigma ) ||
@@ -79,7 +84,27 @@ function logpdf( x, sigma ) {
7984
}
8085
s2 = pow( sigma, 2.0 );
8186
s2i = 1.0 / s2;
82-
return ln( s2i * x ) - ( pow( x, 2.0 ) / ( 2.0 * s2 ) );
87+
x2 = pow( x, 2.0 );
88+
q = s2i * x;
89+
if (
90+
s2 >= FLOAT64_SMALLEST_NORMAL &&
91+
s2 < PINF &&
92+
x2 < PINF &&
93+
q >= FLOAT64_SMALLEST_NORMAL
94+
) {
95+
return ln( q ) - ( x2 / ( 2.0 * s2 ) );
96+
}
97+
// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
98+
r = x / sigma;
99+
q = r / sigma;
100+
if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {
101+
lq = ln( q );
102+
} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {
103+
lq = ln( r ) - ln( sigma );
104+
} else {
105+
lq = ( ln( x ) - ln( sigma ) ) - ln( sigma );
106+
}
107+
return lq - ( 0.5 * r * r );
83108
}
84109

85110

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/manifest.json‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,8 @@
4343
"@stdlib/math/base/special/pow",
4444
"@stdlib/math/base/special/ln",
4545
"@stdlib/constants/float64/pinf",
46-
"@stdlib/constants/float64/ninf"
46+
"@stdlib/constants/float64/ninf",
47+
"@stdlib/constants/float64/smallest-normal"
4748
]
4849
},
4950
{
@@ -63,7 +64,8 @@
6364
"@stdlib/math/base/special/pow",
6465
"@stdlib/math/base/special/ln",
6566
"@stdlib/constants/float64/pinf",
66-
"@stdlib/constants/float64/ninf"
67+
"@stdlib/constants/float64/ninf",
68+
"@stdlib/constants/float64/smallest-normal"
6769
]
6870
},
6971
{
@@ -82,7 +84,8 @@
8284
"@stdlib/math/base/special/pow",
8385
"@stdlib/math/base/special/ln",
8486
"@stdlib/constants/float64/pinf",
85-
"@stdlib/constants/float64/ninf"
87+
"@stdlib/constants/float64/ninf",
88+
"@stdlib/constants/float64/smallest-normal"
8689
]
8790
}
8891
]

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/src/main.c‎

Lines changed: 26 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
#include "stdlib/math/base/special/ln.h"
2323
#include "stdlib/constants/float64/pinf.h"
2424
#include "stdlib/constants/float64/ninf.h"
25+
#include "stdlib/constants/float64/smallest_normal.h"
2526

2627
/**
2728
* Evaluates the natural logarithm of the probability density function (PDF) for a Rayleigh distribution with scale parameter `sigma` at a value `x`.
@@ -36,7 +37,11 @@
3637
*/
3738
double stdlib_base_dists_rayleigh_logpdf( const double x, const double sigma ) {
3839
double s2i;
40+
double lq;
3941
double s2;
42+
double x2;
43+
double r;
44+
double q;
4045
if (
4146
stdlib_base_is_nan( x ) ||
4247
stdlib_base_is_nan( sigma ) ||
@@ -52,5 +57,25 @@ double stdlib_base_dists_rayleigh_logpdf( const double x, const double sigma ) {
5257
}
5358
s2 = stdlib_base_pow( sigma, 2.0 );
5459
s2i = 1.0 / s2;
55-
return stdlib_base_ln( s2i * x ) - stdlib_base_pow( x, 2.0 ) / ( 2.0 * s2 );
60+
x2 = stdlib_base_pow( x, 2.0 );
61+
q = s2i * x;
62+
if (
63+
s2 >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL &&
64+
s2 < STDLIB_CONSTANT_FLOAT64_PINF &&
65+
x2 < STDLIB_CONSTANT_FLOAT64_PINF &&
66+
q >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL
67+
) {
68+
return stdlib_base_ln( q ) - x2 / ( 2.0 * s2 );
69+
}
70+
// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
71+
r = x / sigma;
72+
q = r / sigma;
73+
if ( q >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL && q < STDLIB_CONSTANT_FLOAT64_PINF ) {
74+
lq = stdlib_base_ln( q );
75+
} else if ( r >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
76+
lq = stdlib_base_ln( r ) - stdlib_base_ln( sigma );
77+
} else {
78+
lq = ( stdlib_base_ln( x ) - stdlib_base_ln( sigma ) ) - stdlib_base_ln( sigma );
79+
}
80+
return lq - ( 0.5 * r * r );
5681
}

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/fixtures/julia/huge_scale.json‎

Lines changed: 1 addition & 0 deletions
Large diffs are not rendered by default.

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/fixtures/julia/runner.jl‎

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,6 @@
1616
# See the License for the specific language governing permissions and
1717
# limitations under the License.
1818

19-
import Distributions: logpdf, Rayleigh
2019
import JSON
2120

2221
"""
@@ -39,9 +38,11 @@ julia> gen( x, sigma, \"data.json\" );
3938
```
4039
"""
4140
function gen( x, sigma, name )
42-
z = Array{Float64}( undef, length(x) );
43-
for i in eachindex(x)
44-
z[ i ] = logpdf( Rayleigh( sigma[i] ), x[i] );
41+
# The log PDF is evaluated in extended precision using the textbook formula `ln(x/sigma^2) - x^2/(2*sigma^2)`, as forming `sigma^2` or `x^2` in double precision overflows or underflows for scale parameters far from `1`...
42+
z = setprecision( BigFloat, 2048 ) do
43+
X = BigFloat.( x );
44+
S = BigFloat.( sigma );
45+
Float64.( log.( X ./ ( S .^ 2 ) ) .- ( ( X .^ 2 ) ./ ( 2 .* ( S .^ 2 ) ) ) );
4546
end
4647

4748
# Store data to be written to file as a collection:
@@ -81,3 +82,13 @@ gen( x, sigma, "medium_scale.json" );
8182
x = rand( 1000 ) .* 20.0;
8283
sigma = rand( 1000 ) .* 5.0;
8384
gen( x, sigma, "small_scale.json" );
85+
86+
# Tiny scale parameter (`sigma^2` underflows):
87+
sigma = exp10.( ( rand( 1000 ) .* 45.0 ) .- 200.0 );
88+
x = sigma .* exp10.( ( rand( 1000 ) .* 4.0 ) .- 3.0 );
89+
gen( x, sigma, "tiny_scale.json" );
90+
91+
# Huge scale parameter (`sigma^2` overflows):
92+
sigma = exp10.( ( rand( 1000 ) .* 45.0 ) .+ 155.0 );
93+
x = sigma .* exp10.( ( rand( 1000 ) .* 4.0 ) .- 3.0 );
94+
gen( x, sigma, "huge_scale.json" );

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/fixtures/julia/tiny_scale.json‎

Lines changed: 1 addition & 0 deletions
Large diffs are not rendered by default.

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/test.factory.js‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,8 @@ var factory = require( './../lib/factory.js' );
3232

3333
var smallScale = require( './fixtures/julia/small_scale.json' );
3434
var mediumScale = require( './fixtures/julia/medium_scale.json' );
35+
var tinyScale = require( './fixtures/julia/tiny_scale.json' );
36+
var hugeScale = require( './fixtures/julia/huge_scale.json' );
3537
var largeScale = require( './fixtures/julia/large_scale.json' );
3638

3739

@@ -187,3 +189,45 @@ tape( 'the created function evaluates the logpdf for `x` given large scale param
187189
}
188190
t.end();
189191
});
192+
193+
tape( 'the created function evaluates the logpdf for `x` given tiny scale parameter `sigma`', function test( t ) {
194+
var expected;
195+
var logpdf;
196+
var sigma;
197+
var i;
198+
var x;
199+
var y;
200+
201+
expected = tinyScale.expected;
202+
x = tinyScale.x;
203+
sigma = tinyScale.sigma;
204+
for ( i = 0; i < x.length; i++ ) {
205+
logpdf = factory( sigma[i] );
206+
y = logpdf( x[i] );
207+
if ( expected[i] !== null ) {
208+
t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
209+
}
210+
}
211+
t.end();
212+
});
213+
214+
tape( 'the created function evaluates the logpdf for `x` given huge scale parameter `sigma`', function test( t ) {
215+
var expected;
216+
var logpdf;
217+
var sigma;
218+
var i;
219+
var x;
220+
var y;
221+
222+
expected = hugeScale.expected;
223+
x = hugeScale.x;
224+
sigma = hugeScale.sigma;
225+
for ( i = 0; i < x.length; i++ ) {
226+
logpdf = factory( sigma[i] );
227+
y = logpdf( x[i] );
228+
if ( expected[i] !== null ) {
229+
t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
230+
}
231+
}
232+
t.end();
233+
});

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/test.logpdf.js‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,8 @@ var logpdf = require( './../lib' );
3333

3434
var smallScale = require( './fixtures/julia/small_scale.json' );
3535
var mediumScale = require( './fixtures/julia/medium_scale.json' );
36+
var tinyScale = require( './fixtures/julia/tiny_scale.json' );
37+
var hugeScale = require( './fixtures/julia/huge_scale.json' );
3638
var largeScale = require( './fixtures/julia/large_scale.json' );
3739

3840

@@ -180,3 +182,57 @@ tape( 'the function evaluates the logpdf for `x` given large scale parameter `si
180182
}
181183
t.end();
182184
});
185+
186+
tape( 'the function evaluates the logpdf for `x` given tiny scale parameter `sigma`', function test( t ) {
187+
var expected;
188+
var delta;
189+
var sigma;
190+
var tol;
191+
var i;
192+
var x;
193+
var y;
194+
195+
expected = tinyScale.expected;
196+
x = tinyScale.x;
197+
sigma = tinyScale.sigma;
198+
for ( i = 0; i < x.length; i++ ) {
199+
y = logpdf( x[i], sigma[i] );
200+
if ( expected[i] !== null ) {
201+
if ( y === expected[i] ) {
202+
t.strictEqual( y, expected[i], 'x: '+x[i]+', sigma: '+sigma[i]+', y: '+y+', expected: '+expected[i] );
203+
} else {
204+
delta = abs( y - expected[ i ] );
205+
tol = 3.0 * EPS * abs( expected[ i ] );
206+
t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. sigma: '+sigma[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
207+
}
208+
}
209+
}
210+
t.end();
211+
});
212+
213+
tape( 'the function evaluates the logpdf for `x` given huge scale parameter `sigma`', function test( t ) {
214+
var expected;
215+
var delta;
216+
var sigma;
217+
var tol;
218+
var i;
219+
var x;
220+
var y;
221+
222+
expected = hugeScale.expected;
223+
x = hugeScale.x;
224+
sigma = hugeScale.sigma;
225+
for ( i = 0; i < x.length; i++ ) {
226+
y = logpdf( x[i], sigma[i] );
227+
if ( expected[i] !== null ) {
228+
if ( y === expected[i] ) {
229+
t.strictEqual( y, expected[i], 'x: '+x[i]+', sigma: '+sigma[i]+', y: '+y+', expected: '+expected[i] );
230+
} else {
231+
delta = abs( y - expected[ i ] );
232+
tol = 3.0 * EPS * abs( expected[ i ] );
233+
t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. sigma: '+sigma[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
234+
}
235+
}
236+
}
237+
t.end();
238+
});

‎lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/test.native.js‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,8 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
3333

3434
var smallScale = require( './fixtures/julia/small_scale.json' );
3535
var mediumScale = require( './fixtures/julia/medium_scale.json' );
36+
var tinyScale = require( './fixtures/julia/tiny_scale.json' );
37+
var hugeScale = require( './fixtures/julia/huge_scale.json' );
3638
var largeScale = require( './fixtures/julia/large_scale.json' );
3739

3840

@@ -164,3 +166,41 @@ tape( 'the function evaluates the logpdf for `x` given large scale parameter `si
164166
}
165167
t.end();
166168
});
169+
170+
tape( 'the function evaluates the logpdf for `x` given tiny scale parameter `sigma`', opts, function test( t ) {
171+
var expected;
172+
var sigma;
173+
var i;
174+
var x;
175+
var y;
176+
177+
expected = tinyScale.expected;
178+
x = tinyScale.x;
179+
sigma = tinyScale.sigma;
180+
for ( i = 0; i < x.length; i++ ) {
181+
y = logpdf( x[i], sigma[i] );
182+
if ( expected[i] !== null ) {
183+
t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
184+
}
185+
}
186+
t.end();
187+
});
188+
189+
tape( 'the function evaluates the logpdf for `x` given huge scale parameter `sigma`', opts, function test( t ) {
190+
var expected;
191+
var sigma;
192+
var i;
193+
var x;
194+
var y;
195+
196+
expected = hugeScale.expected;
197+
x = hugeScale.x;
198+
sigma = hugeScale.sigma;
199+
for ( i = 0; i < x.length; i++ ) {
200+
y = logpdf( x[i], sigma[i] );
201+
if ( expected[i] !== null ) {
202+
t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
203+
}
204+
}
205+
t.end();
206+
});

0 commit comments

Comments
 (0)