Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor: update stats/base/dists/triangular/mgf implementation and test tolerances #4768

Merged
Merged
Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
51e7260
feat: update fixtures for triangular MGF based on Julia implementation
anandkaranubc Jan 15, 2025
f4234cc
refactor: align main and factory.js with Julia-based triangular MGF i…
anandkaranubc Jan 15, 2025
078ba22
test: adjust tolerance levels in test.native.js for triangular MGF co…
anandkaranubc Jan 15, 2025
263cdaf
refactor: update main.c to align triangular MGF implementation with J…
anandkaranubc Jan 15, 2025
004d18c
refactor: reorganise unused variables in main.c
anandkaranubc Jan 15, 2025
ff1d644
feat: add fixtures for triangular MGF test cases for different condit…
anandkaranubc Jan 16, 2025
aca10b4
feat: add fixtures for triangular MGF test cases for different condit…
anandkaranubc Jan 16, 2025
01ce172
bench: added benchmarks for different cases in triangular mgf impleme…
anandkaranubc Jan 16, 2025
d47eda3
refactor: remove the extra if in main.c
anandkaranubc Jan 16, 2025
7267f00
chore: fix the order of includes in main.c
anandkaranubc Jan 16, 2025
8f01ad8
refactor: apply suggestions from PR review
anandkaranubc Jan 21, 2025
0937887
refactor: add private label to phi2.js
anandkaranubc Jan 22, 2025
2d4f12d
bench: update tolerances for the native tests
anandkaranubc Jan 22, 2025
659542c
bench: move random number generation outside the benchmarking loops a…
anandkaranubc Jan 22, 2025
d36338c
Merge branch 'stdlib-js:develop' into fix/triangular-mgf-divergence
anandkaranubc Mar 1, 2025
63bc5ad
test: update comment about high tolerance issue
anandkaranubc Mar 1, 2025
874afd7
Merge remote-tracking branch 'upstream/develop' into fix/triangular-m…
stdlib-bot Mar 1, 2025
0e26b69
chore: add exports header comment
Planeshifter Mar 6, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@

var bench = require( '@stdlib/bench' );
var Float64Array = require( '@stdlib/array/float64' );
var randu = require( '@stdlib/random/base/randu' );
var uniform = require( '@stdlib/random/base/uniform' );
var isnan = require( '@stdlib/math/base/assert/is-nan' );
var EPS = require( '@stdlib/constants/float64/eps' );
var pkg = require( './../package.json' ).name;
Expand All @@ -47,10 +47,10 @@ bench( pkg, function benchmark( b ) {
mode = new Float64Array( len );

for ( i = 0; i < len; i++ ) {
t[ i ] = randu() * 5.0;
min[ i ] = randu() * 10.0;
max[ i ] = min[ i ] + ( randu() * 40.0 ) + EPS;
mode[ i ] = min[ i ] + ( ( max[ i ] - min[ i ] ) * randu() );
t[ i ] = uniform( 0.0, 5.0 );
min[ i ] = uniform( 0.0, 10.0 );
max[ i ] = uniform( min[ i ] + EPS, min[ i ] + 40.0 );
mode[ i ] = uniform( min[ i ], max[ i ]);
}

b.tic();
Expand All @@ -73,6 +73,7 @@ bench( pkg+':factory', function benchmark( b ) {
var mode;
var min;
var max;
var len;
var t;
var y;
var i;
Expand All @@ -81,11 +82,15 @@ bench( pkg+':factory', function benchmark( b ) {
max = 1.5;
mode = 0.5;
mymgf = mgf.factory( min, max, mode );
len = 100;
t = new Float64Array( len );
for ( i = 0; i < len; i++ ) {
t[ i ] = uniform( 0.0, 5.0 );
}

b.tic();
for ( i = 0; i < b.iterations; i++ ) {
t = randu() * 5.0;
y = mymgf( t );
y = mymgf( t[ i % len ]);
if ( isnan( y ) ) {
b.fail( 'should not return NaN' );
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
var resolve = require( 'path' ).resolve;
var bench = require( '@stdlib/bench' );
var Float64Array = require( '@stdlib/array/float64' );
var randu = require( '@stdlib/random/base/randu' );
var uniform = require( '@stdlib/random/base/uniform' );
var isnan = require( '@stdlib/math/base/assert/is-nan' );
var tryRequire = require( '@stdlib/utils/try-require' );
var EPS = require( '@stdlib/constants/float64/eps' );
Expand Down Expand Up @@ -56,10 +56,10 @@ bench( pkg+'::native', opts, function benchmark( b ) {
mode = new Float64Array( len );

for ( i = 0; i < len; i++ ) {
t[ i ] = randu() * 5.0;
min[ i ] = randu() * 10.0;
max[ i ] = min[ i ] + ( randu() * 40.0 ) + EPS;
mode[ i ] = min[ i ] + ( ( max[ i ] - min[ i ] ) * randu() );
t[ i ] = uniform( 0.0, 5.0 );
min[ i ] = uniform( 0.0, 10.0 );
max[ i ] = uniform( min[ i ] + EPS, min[ i ] + 40.0 );
mode[ i ] = uniform( min[ i ], max[ i ]);
}

b.tic();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
var constantFunction = require( '@stdlib/utils/constant-function' );
var isnan = require( '@stdlib/math/base/assert/is-nan' );
var exp = require( '@stdlib/math/base/special/exp' );
var pow = require( '@stdlib/math/base/special/pow' );
var phi2 = require( './phi2.js' );


// MAIN //
Expand All @@ -45,10 +45,6 @@ var pow = require( '@stdlib/math/base/special/pow' );
* // returns ~10.205
*/
function factory( a, b, c ) {
var bmc;
var bma;
var cma;

if (
isnan( a ) ||
isnan( b ) ||
Expand All @@ -58,9 +54,6 @@ function factory( a, b, c ) {
) {
return constantFunction( NaN );
}
bmc = b - c;
bma = b - a;
cma = c - a;
return mgf;

/**
Expand All @@ -75,19 +68,19 @@ function factory( a, b, c ) {
* // returns <number>
*/
function mgf( t ) {
var ret;

if ( isnan( t ) ) {
return NaN;
}
if ( t === 0.0 ) {
return 1.0;
if ( a < c ) {
if ( c < b ) {
return exp( c*t ) * ( ( (c-a)*phi2( (a-c)*t ) ) + ( (b-c)*phi2( (b-c)*t ) ) ) / ( b-a ); // eslint-disable-line max-len
}
return exp( c*t ) * phi2( ( a-c ) * t );
}
if ( c < b ) {
return exp( c*t ) * phi2( ( b-c ) * t );
}
ret = ( bmc * exp( a*t ) ) - ( bma * exp( c*t ) );
ret += cma * exp( b*t );
ret *= 2.0;
ret /= bma * cma * bmc * pow( t, 2.0 );
return ret;
return exp( c*t );
}
}

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@

var isnan = require( '@stdlib/math/base/assert/is-nan' );
var exp = require( '@stdlib/math/base/special/exp' );
var pow = require( '@stdlib/math/base/special/pow' );
var phi2 = require( './phi2.js' );


// MAIN //
Expand Down Expand Up @@ -73,11 +73,6 @@ var pow = require( '@stdlib/math/base/special/pow' );
* // returns NaN
*/
function mgf( t, a, b, c ) {
var bmc;
var bma;
var cma;
var ret;

if (
isnan( t ) ||
isnan( a ) ||
Expand All @@ -88,17 +83,16 @@ function mgf( t, a, b, c ) {
) {
return NaN;
}
if ( t === 0.0 ) {
return 1.0;
if ( a < c ) {
if ( c < b ) {
return exp( c*t ) * ( ( (c-a)*phi2( (a-c)*t ) ) + ( (b-c)*phi2( (b-c)*t ) ) ) / ( b-a ); // eslint-disable-line max-len
}
return exp( c*t ) * phi2( ( a-c ) * t );
}
if ( c < b ) {
return exp( c*t ) * phi2( ( b-c ) * t );
}
bmc = b - c;
bma = b - a;
cma = c - a;
ret = ( bmc * exp( a*t ) ) - ( bma * exp( c*t ) );
ret += cma * exp( b*t );
ret *= 2.0;
ret /= bma * cma * bmc * pow( t, 2.0 );
return ret;
return exp( c*t );
}


Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
/**
* @license Apache-2.0
*
* Copyright (c) 2025 The Stdlib Authors.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

'use strict';

// MODULES //

var expm1 = require( '@stdlib/math/base/special/expm1' );


// MAIN //

/**
* Helper function for repeated computation in the MGF formula.
*
* @private
* @param {number} x - input value
* @returns {number} evaluated result
*/
function phi2( x ) {
if ( x === 0.0 ) {
return 1.0;
}
return ( 2.0 * ( expm1( x ) - x ) ) / ( x*x );
}

module.exports = phi2;
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@
"@stdlib/math/base/napi/quaternary",
"@stdlib/math/base/assert/is-nan",
"@stdlib/math/base/special/exp",
"@stdlib/math/base/special/pow"
"@stdlib/math/base/special/expm1"
]
},
{
Expand Down
51 changes: 27 additions & 24 deletions lib/node_modules/@stdlib/stats/base/dists/triangular/mgf/src/main.c
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,21 @@

#include "stdlib/stats/base/dists/triangular/mgf.h"
#include "stdlib/math/base/assert/is_nan.h"
#include "stdlib/math/base/special/expm1.h"
#include "stdlib/math/base/special/exp.h"
#include "stdlib/math/base/special/pow.h"

/**
* Helper function for repeated computation in the MGF formula.
*
* @param x input value
* @return evaluated result
*/
static double phi2( const double x ) {
if ( x == 0.0 ) {
return 1.0;
}
return ( 2.0 * ( stdlib_base_expm1( x ) - x ) ) / ( x*x );
}

/**
* Evaluates the moment-generating function (MGF) for a triangular distribution with lower limit `a`, upper limit `b`, and mode `c` at a value `t`.
Expand All @@ -35,29 +48,19 @@
* // returns ~1.021
*/
double stdlib_base_dists_triangular_mgf( const double t, const double a, const double b, const double c ) {
double bmc;
double bma;
double cma;
double ret;
if (
stdlib_base_is_nan( t ) ||
stdlib_base_is_nan( a ) ||
stdlib_base_is_nan( b ) ||
stdlib_base_is_nan( c ) ||
a > c ||
c > b
) {
return 0.0/0.0; // NaN
if ( stdlib_base_is_nan( t ) || stdlib_base_is_nan( a ) || stdlib_base_is_nan( b ) || stdlib_base_is_nan( c ) || a > c || c > b ) {
return 0.0 / 0.0; // NaN
}
if ( t == 0.0 ) {
return 1.0;
if ( a < c ) {
if ( c < b ) {
double term1 = ( c-a ) * phi2( ( a-c ) * t );
double term2 = ( b-c ) * phi2( ( b-c ) * t );
return stdlib_base_exp( c*t ) * ( term1+term2 ) / ( b-a );
}
return stdlib_base_exp( c*t ) * phi2( ( a-c ) * t );
}
if ( c < b ) {
return stdlib_base_exp( c*t ) * phi2( ( b-c ) * t );
}
bmc = b - c;
bma = b - a;
cma = c - a;
ret = ( bmc * stdlib_base_exp( a*t ) ) - ( bma * stdlib_base_exp( c*t ) );
ret += cma * stdlib_base_exp( b*t );
ret *= 2.0;
ret /= bma * cma * bmc * stdlib_base_pow( t, 2.0 );
return ret;
return stdlib_base_exp( c*t );
}

Large diffs are not rendered by default.

Large diffs are not rendered by default.

Large diffs are not rendered by default.

Large diffs are not rendered by default.

Large diffs are not rendered by default.

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
Expand Up @@ -93,3 +93,15 @@ b = ( rand( 1000 ) .* 80.0 ) .+ a;
c = a .+ ( b .- a ) .* rand();
x = rand( 1000 );
gen( x, a, b, c, "large_range.json" );

# Case: a < c < b
gen( x, a, b, c, "a_less_c_less_b.json" );

# Case: a < c = b
gen( x, a, b, b, "a_less_c_eq_b.json" );

# Case: a = c < b
gen( x, a, b, a, "a_eq_c_less_b.json" );

# Case: a = c = b
gen( x, a, a, a, "a_eq_c_eq_b.json" );

Large diffs are not rendered by default.

Loading