diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/README.md b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/README.md index 9322e3a3c92c..aeedb85d557f 100644 --- a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/README.md +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/README.md @@ -143,6 +143,106 @@ logEachMap( 'x: %0.4f, d1: %0.4f, d2: %0.4f, f(x;d1,d2): %0.4f', x, d1, d2, pdf + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/stats/base/dists/f/pdf.h" +``` + +#### stdlib_base_dists_f_pdf( x, d1, d2 ) + +Evaluates the probability density function (PDF) for an F distribution with numerator degrees of freedom `d1` and denominator degrees of freedom `d2` at a value `x`. + + +```c +double out = stdlib_base_dists_f_pdf( 2.0, 0.5, 1.0 ); +// returns ~0.057 +``` + +The function accepts the following arguments: + +- **x**: `[in] double` input value. +- **d1**: `[in] double` numerator degrees of freedom. +- **d2**: `[in] double` denominator degrees of freedom. + +```c +double stdlib_base_dists_f_pdf( const double x, const double d1, const double d2 ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/stats/base/dists/f/pdf.h" +#include +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + double x; + double d1; + double d2; + double y; + int32_t i; + + for ( i = 0; i < 10; i++ ) { + x = random_uniform( 0.0, 4.0 ); + d1 = random_uniform( 0.0, 10.0 ); + d2 = random_uniform( 0.0, 10.0 ); + y = stdlib_base_dists_f_pdf( x, d1, d2 ); + printf( "x: %lf, d1: %lf, d2: %lf, F(x;d1,d2): %lf\n", x, d1, d2, y ); + } +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/benchmark.native.js new file mode 100644 index 000000000000..a66d0eafbc97 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/benchmark.native.js @@ -0,0 +1,76 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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 resolve = require( 'path' ).resolve; +var tryRequire = require( '@stdlib/utils/try-require' ); +var bench = require( '@stdlib/bench' ); +var Float64Array = require( '@stdlib/array/float64' ); +var uniform = require( '@stdlib/random/base/uniform' ); +var discreteUniform = require( '@stdlib/random/base/discrete-uniform' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var EPS = require( '@stdlib/constants/float64/eps' ); +var format = require( '@stdlib/string/format' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var pdf = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( pdf instanceof Error ) +}; + + +// MAIN // + +bench( format( '%s::native', pkg ), opts, function benchmark( b ) { + var len; + var d1; + var d2; + var x; + var y; + var i; + + len = 100; + x = new Float64Array( len ); + d1 = new Float64Array( len ); + d2 = new Float64Array( len ); + for ( i = 0; i < len; i++ ) { + x[ i ] = uniform( EPS, 100.0 ); + d1[ i ] = discreteUniform( 1, 100 ); + d2[ i ] = discreteUniform( 1, 100 ); + } + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + y = pdf( x[ i % len ], d1[ i % len ], d2[ i % len ]); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + b.pass( 'benchmark finished' ); + b.end(); +}); diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/c/Makefile b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/c/Makefile new file mode 100644 index 000000000000..979768abbcec --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/c/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2026 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. +#/ + +# VARIABLES # + +ifndef VERBOSE + QUIET := @ +else + QUIET := +endif + +# Determine the OS ([1][1], [2][2]). +# +# [1]: https://en.wikipedia.org/wiki/Uname#Examples +# [2]: http://stackoverflow.com/a/27776822/2225624 +OS ?= $(shell uname) +ifneq (, $(findstring MINGW,$(OS))) + OS := WINNT +else +ifneq (, $(findstring MSYS,$(OS))) + OS := WINNT +else +ifneq (, $(findstring CYGWIN,$(OS))) + OS := WINNT +else +ifneq (, $(findstring Windows_NT,$(OS))) + OS := WINNT +endif +endif +endif +endif + +# Define the program used for compiling C source files: +ifdef C_COMPILER + CC := $(C_COMPILER) +else + CC := gcc +endif + +# Define the command-line options when compiling C files: +CFLAGS ?= \ + -std=c99 \ + -O3 \ + -Wall \ + -pedantic + +# Determine whether to generate position independent code ([1][1], [2][2]). +# +# [1]: https://gcc.gnu.org/onlinedocs/gcc/Code-Gen-Options.html#Code-Gen-Options +# [2]: http://stackoverflow.com/questions/5311515/gcc-fpic-option +ifeq ($(OS), WINNT) + fPIC ?= +else + fPIC ?= -fPIC +endif + +# List of includes (e.g., `-I /foo/bar -I /beep/boop/include`): +INCLUDE ?= + +# List of source files: +SOURCE_FILES ?= + +# List of libraries (e.g., `-lopenblas -lpthread`): +LIBRARIES ?= + +# List of library paths (e.g., `-L /foo/bar -L /beep/boop`): +LIBPATH ?= + +# List of C targets: +c_targets := benchmark.out + + +# RULES # + +#/ +# Compiles source files. +# +# @param {string} [C_COMPILER] - C compiler (e.g., `gcc`) +# @param {string} [CFLAGS] - C compiler options +# @param {(string|void)} [fPIC] - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} [INCLUDE] - list of includes (e.g., `-I /foo/bar -I /beep/boop/include`) +# @param {string} [SOURCE_FILES] - list of source files +# @param {string} [LIBPATH] - list of library paths (e.g., `-L /foo/bar -L /beep/boop`) +# @param {string} [LIBRARIES] - list of libraries (e.g., `-lopenblas -lpthread`) +# +# @example +# make +# +# @example +# make all +#/ +all: $(c_targets) + +.PHONY: all + +#/ +# Compiles C source files. +# +# @private +# @param {string} CC - C compiler (e.g., `gcc`) +# @param {string} CFLAGS - C compiler options +# @param {(string|void)} fPIC - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} INCLUDE - list of includes (e.g., `-I /foo/bar`) +# @param {string} SOURCE_FILES - list of source files +# @param {string} LIBPATH - list of library paths (e.g., `-L /foo/bar`) +# @param {string} LIBRARIES - list of libraries (e.g., `-lopenblas`) +#/ +$(c_targets): %.out: %.c + $(QUIET) $(CC) $(CFLAGS) $(fPIC) $(INCLUDE) -o $@ $(SOURCE_FILES) $< $(LIBPATH) -lm $(LIBRARIES) + +#/ +# Runs compiled benchmarks. +# +# @example +# make run +#/ +run: $(c_targets) + $(QUIET) ./$< + +.PHONY: run + +#/ +# Removes generated files. +# +# @example +# make clean +#/ +clean: + $(QUIET) -rm -f *.o *.out + +.PHONY: clean diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/c/benchmark.c b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/c/benchmark.c new file mode 100644 index 000000000000..9025f9b0a82f --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/benchmark/c/benchmark.c @@ -0,0 +1,144 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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. +*/ + +#include "stdlib/stats/base/dists/f/pdf.h" +#include "stdlib/constants/float64/eps.h" +#include "stdlib/math/base/special/ceil.h" +#include +#include +#include +#include +#include + +#define NAME "f-pdf" +#define ITERATIONS 1000000 +#define REPEATS 3 + +/** +* Prints the TAP version. +*/ +static void print_version( void ) { + printf( "TAP version 13\n" ); +} + +/** +* Prints the TAP summary. +* +* @param total total number of tests +* @param passing total number of passing tests +*/ +static void print_summary( int total, int passing ) { + printf( "#\n" ); + printf( "1..%d\n", total ); // TAP plan + printf( "# total %d\n", total ); + printf( "# pass %d\n", passing ); + printf( "#\n" ); + printf( "# ok\n" ); +} + +/** +* Prints benchmarks results. +* +* @param elapsed elapsed time in seconds +*/ +static void print_results( double elapsed ) { + double rate = (double)ITERATIONS / elapsed; + printf( " ---\n" ); + printf( " iterations: %d\n", ITERATIONS ); + printf( " elapsed: %0.9f\n", elapsed ); + printf( " rate: %0.9f\n", rate ); + printf( " ...\n" ); +} + +/** +* Returns a clock time. +* +* @return clock time +*/ +static double tic( void ) { + struct timeval now; + gettimeofday( &now, NULL ); + return (double)now.tv_sec + (double)now.tv_usec/1.0e6; +} + +/** +* Generates a random number on the interval [min,max). +* +* @param min minimum value (inclusive) +* @param max maximum value (exclusive) +* @return random number +*/ +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +/** +* Runs a benchmark. +* +* @return elapsed time in seconds +*/ +static double benchmark( void ) { + double elapsed; + double x[ 100 ]; + double d1[ 100 ]; + double d2[ 100 ]; + double y; + double t; + int i; + + for ( i = 0; i < 100; i++ ) { + x[ i ] = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 100.0 ); + d1[ i ] = stdlib_base_ceil( random_uniform( 1.0, 100.0 ) ); + d2[ i ] = stdlib_base_ceil( random_uniform( 1.0, 100.0 ) ); + } + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_dists_f_pdf( x[ i%100 ], d1[ i%100 ], d2[ i%100 ] ); + if ( y != y ) { + printf( "should not return NaN\n" ); + break; + } + } + elapsed = tic() - t; + if ( y != y ) { + printf( "should not return NaN\n" ); + } + return elapsed; +} + +/** +* Main execution sequence. +*/ +int main( void ) { + double elapsed; + int i; + + // Use the current time to seed the random number generator: + srand( time( NULL ) ); + + print_version(); + for ( i = 0; i < REPEATS; i++ ) { + printf( "# c::%s\n", NAME ); + elapsed = benchmark(); + print_results( elapsed ); + printf( "ok %d benchmark finished\n", i+1 ); + } + print_summary( REPEATS, REPEATS ); +} diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/binding.gyp b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/binding.gyp new file mode 100644 index 000000000000..0d6508a12e99 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/binding.gyp @@ -0,0 +1,170 @@ +# @license Apache-2.0 +# +# Copyright (c) 2026 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. + +# A `.gyp` file for building a Node.js native add-on. +# +# [1]: https://gyp.gsrc.io/docs/InputFormatReference.md +# [2]: https://gyp.gsrc.io/docs/UserDocumentation.md +{ + # List of files to include in this file: + 'includes': [ + './include.gypi', + ], + + # Define variables to be used throughout the configuration for all targets: + 'variables': { + # Target name should match the add-on export name: + 'addon_target_name%': 'addon', + + # Set variables based on the host OS: + 'conditions': [ + [ + 'OS=="win"', + { + # Define the object file suffix: + 'obj': 'obj', + }, + { + # Define the object file suffix: + 'obj': 'o', + } + ], # end condition (OS=="win") + ], # end conditions + }, # end variables + + # Define compile targets: + 'targets': [ + + # Target to generate an add-on: + { + # The target name should match the add-on export name: + 'target_name': '<(addon_target_name)', + + # Define dependencies: + 'dependencies': [], + + # Define directories which contain relevant include headers: + 'include_dirs': [ + # Local include directory: + '<@(include_dirs)', + ], + + # List of source files: + 'sources': [ + '<@(src_files)', + ], + + # Settings which should be applied when a target's object files are used as linker input: + 'link_settings': { + # Define libraries: + 'libraries': [ + '<@(libraries)', + ], + + # Define library directories: + 'library_dirs': [ + '<@(library_dirs)', + ], + }, + + # C/C++ compiler flags: + 'cflags': [ + # Enable commonly used warning options: + '-Wall', + + # Aggressive optimization: + '-O3', + ], + + # C specific compiler flags: + 'cflags_c': [ + # Specify the C standard to which a program is expected to conform: + '-std=c99', + ], + + # C++ specific compiler flags: + 'cflags_cpp': [ + # Specify the C++ standard to which a program is expected to conform: + '-std=c++11', + ], + + # Linker flags: + 'ldflags': [], + + # Apply conditions based on the host OS: + 'conditions': [ + [ + 'OS=="mac"', + { + # Linker flags: + 'ldflags': [ + '-undefined dynamic_lookup', + '-Wl,-no-pie', + '-Wl,-search_paths_first', + ], + }, + ], # end condition (OS=="mac") + [ + 'OS!="win"', + { + # C/C++ flags: + 'cflags': [ + # Generate platform-independent code: + '-fPIC', + ], + }, + ], # end condition (OS!="win") + ], # end conditions + }, # end target <(addon_target_name) + + # Target to copy a generated add-on to a standard location: + { + 'target_name': 'copy_addon', + + # Declare that the output of this target is not linked: + 'type': 'none', + + # Define dependencies: + 'dependencies': [ + # Require that the add-on be generated before building this target: + '<(addon_target_name)', + ], + + # Define a list of actions: + 'actions': [ + { + 'action_name': 'copy_addon', + 'message': 'Copying addon...', + + # Explicitly list the inputs in the command-line invocation below: + 'inputs': [], + + # Declare the expected outputs: + 'outputs': [ + '<(addon_output_dir)/<(addon_target_name).node', + ], + + # Define the command-line invocation: + 'action': [ + 'cp', + '<(PRODUCT_DIR)/<(addon_target_name).node', + '<(addon_output_dir)/<(addon_target_name).node', + ], + }, + ], # end actions + }, # end target copy_addon + ], # end targets +} diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/examples/c/Makefile b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/examples/c/Makefile new file mode 100644 index 000000000000..c8f8e9a1517b --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/examples/c/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2026 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. +#/ + +# VARIABLES # + +ifndef VERBOSE + QUIET := @ +else + QUIET := +endif + +# Determine the OS ([1][1], [2][2]). +# +# [1]: https://en.wikipedia.org/wiki/Uname#Examples +# [2]: http://stackoverflow.com/a/27776822/2225624 +OS ?= $(shell uname) +ifneq (, $(findstring MINGW,$(OS))) + OS := WINNT +else +ifneq (, $(findstring MSYS,$(OS))) + OS := WINNT +else +ifneq (, $(findstring CYGWIN,$(OS))) + OS := WINNT +else +ifneq (, $(findstring Windows_NT,$(OS))) + OS := WINNT +endif +endif +endif +endif + +# Define the program used for compiling C source files: +ifdef C_COMPILER + CC := $(C_COMPILER) +else + CC := gcc +endif + +# Define the command-line options when compiling C files: +CFLAGS ?= \ + -std=c99 \ + -O3 \ + -Wall \ + -pedantic + +# Determine whether to generate position independent code ([1][1], [2][2]). +# +# [1]: https://gcc.gnu.org/onlinedocs/gcc/Code-Gen-Options.html#Code-Gen-Options +# [2]: http://stackoverflow.com/questions/5311515/gcc-fpic-option +ifeq ($(OS), WINNT) + fPIC ?= +else + fPIC ?= -fPIC +endif + +# List of includes (e.g., `-I /foo/bar -I /beep/boop/include`): +INCLUDE ?= + +# List of source files: +SOURCE_FILES ?= + +# List of libraries (e.g., `-lopenblas -lpthread`): +LIBRARIES ?= + +# List of library paths (e.g., `-L /foo/bar -L /beep/boop`): +LIBPATH ?= + +# List of C targets: +c_targets := example.out + + +# RULES # + +#/ +# Compiles source files. +# +# @param {string} [C_COMPILER] - C compiler (e.g., `gcc`) +# @param {string} [CFLAGS] - C compiler options +# @param {(string|void)} [fPIC] - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} [INCLUDE] - list of includes (e.g., `-I /foo/bar -I /beep/boop/include`) +# @param {string} [SOURCE_FILES] - list of source files +# @param {string} [LIBPATH] - list of library paths (e.g., `-L /foo/bar -L /beep/boop`) +# @param {string} [LIBRARIES] - list of libraries (e.g., `-lopenblas -lpthread`) +# +# @example +# make +# +# @example +# make all +#/ +all: $(c_targets) + +.PHONY: all + +#/ +# Compiles C source files. +# +# @private +# @param {string} CC - C compiler (e.g., `gcc`) +# @param {string} CFLAGS - C compiler options +# @param {(string|void)} fPIC - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} INCLUDE - list of includes (e.g., `-I /foo/bar`) +# @param {string} SOURCE_FILES - list of source files +# @param {string} LIBPATH - list of library paths (e.g., `-L /foo/bar`) +# @param {string} LIBRARIES - list of libraries (e.g., `-lopenblas`) +#/ +$(c_targets): %.out: %.c + $(QUIET) $(CC) $(CFLAGS) $(fPIC) $(INCLUDE) -o $@ $(SOURCE_FILES) $< $(LIBPATH) -lm $(LIBRARIES) + +#/ +# Runs compiled examples. +# +# @example +# make run +#/ +run: $(c_targets) + $(QUIET) ./$< + +.PHONY: run + +#/ +# Removes generated files. +# +# @example +# make clean +#/ +clean: + $(QUIET) -rm -f *.o *.out + +.PHONY: clean diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/examples/c/example.c b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/examples/c/example.c new file mode 100644 index 000000000000..8321b7af9285 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/examples/c/example.c @@ -0,0 +1,43 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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. +*/ + +#include "stdlib/stats/base/dists/f/pdf.h" +#include +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + double x; + double d1; + double d2; + double y; + int32_t i; + + for ( i = 0; i < 10; i++ ) { + x = random_uniform( 0.0, 4.0 ); + d1 = random_uniform( 0.0, 10.0 ); + d2 = random_uniform( 0.0, 10.0 ); + y = stdlib_base_dists_f_pdf( x, d1, d2 ); + printf( "x: %lf, d1: %lf, d2: %lf, F(x;d1,d2): %lf\n", x, d1, d2, y ); + } +} diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/include.gypi b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/include.gypi new file mode 100644 index 000000000000..bee8d41a2caf --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/include.gypi @@ -0,0 +1,53 @@ +# @license Apache-2.0 +# +# Copyright (c) 2026 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. + +# A GYP include file for building a Node.js native add-on. +# +# Main documentation: +# +# [1]: https://gyp.gsrc.io/docs/InputFormatReference.md +# [2]: https://gyp.gsrc.io/docs/UserDocumentation.md +{ + # Define variables to be used throughout the configuration for all targets: + 'variables': { + # Source directory: + 'src_dir': './src', + + # Include directories: + 'include_dirs': [ + ' 0 \\) or \\( \operatorname{min}( a, b ) < 1 \\), the two power terms both go in the same direction (towards zero or towards infinity). In this case if either term overflows or underflows, then the product of the two must do so also. Alternatively, if one exponent is less than one, then we can't productively use it to eliminate overflow or underflow from the other term. Problems with spurious overflow/underflow can't be ruled out in this case, but it is _very_ unlikely since one of the power terms will evaluate to a number close to 1. +* +* - If \\( \max( \abs(l1), \abs(l2) ) < 0.5 \\), both exponents are near one and both the exponents are greater than one, and, further, these two power terms tend in opposite directions (one toward zero, the other toward infinity), so we have to combine the terms to avoid any risk of overflow or underflow. We do this by moving one power term inside the other, we have: +* +* ```tex +* (1 + l_1)^a \cdot (1 + l_2)^b \\ +* = ((1 + l_1) \cdot (1 + l_2)^(b/a))^a \\ +* = (1 + l_1 + l_3 + l_1*l_3)^a +* ``` +* +* and +* +* ```tex +* l_3 = (1 + l_2)^(b/a) - 1 \\ +* = \exp((b/a) * \ln(1 + l_2)) - 1 +* ``` +* +* The tricky bit is deciding which term to move inside. By preference we move the larger term inside, so that the size of the largest exponent is reduced. However, that can only be done as long as l3 (see above) is also small. +* +* @private +* @param a - function parameter +* @param b - function parameter +* @param x - function parameter +* @param y - probability equal to `1-x` +* @param normalized - boolean indicating whether to evaluate the power terms of the regularized or non-regularized incomplete beta function +* @returns power terms +*/ +static double ibetaPowerTerms( const double a, const double b, const double x, const double y, const bool normalized ) { + double result; + double smallA; + double ratio; + double agh; + double bgh; + double cgh; + double l1; + double l2; + double l3; + double p1; + double b1; + double b2; + double c; + double l; + + if ( normalized == 0 ) { + // Can we do better here? + return stdlib_base_pow( x, a ) * stdlib_base_pow( y, b ); + } + c = a + b; + + // Combine power terms with Lanczos approximation: + agh = a + STDLIB_CONSTANT_FLOAT64_GAMMA_LANCZOS_G - 0.5; + bgh = b + STDLIB_CONSTANT_FLOAT64_GAMMA_LANCZOS_G - 0.5; + cgh = c + STDLIB_CONSTANT_FLOAT64_GAMMA_LANCZOS_G - 0.5; + result = stdlib_base_gamma_lanczos_sum_expg_scaled( c ); + result /= stdlib_base_gamma_lanczos_sum_expg_scaled( a ) * stdlib_base_gamma_lanczos_sum_expg_scaled( b ); + + // Combine with the leftover terms from the Lanczos approximation: + result *= stdlib_base_sqrt( bgh / STDLIB_CONSTANT_FLOAT64_E ); + result *= stdlib_base_sqrt( agh / cgh ); + + // `l1` and `l2` are the base of the exponents minus one: + l1 = ( ( x * b ) - ( y * agh ) ) / agh; + l2 = ( ( y * a ) - ( x * bgh ) ) / bgh; + if ( stdlib_base_min( stdlib_base_abs(l1), stdlib_base_abs(l2) ) < 0.2 ) { + // When the base of the exponent is very near 1 we get really gross errors unless extra care is taken: + if ( l1 * l2 > 0 || stdlib_base_min( a, b ) < 1 ) { + if ( stdlib_base_abs(l1) < 0.1 ) { + result *= stdlib_base_exp( a * stdlib_base_log1p( l1 ) ); + } else { + result *= stdlib_base_pow( ( x*cgh ) / agh, a ); + } + if ( stdlib_base_abs(l2) < 0.1 ) { + result *= stdlib_base_exp( b * stdlib_base_log1p( l2 ) ); + } else { + result *= stdlib_base_pow((y * cgh) / bgh, b); + } + } + else if ( stdlib_base_max( stdlib_base_abs(l1), stdlib_base_abs(l2) ) < 0.5 ) { + smallA = a < b; + ratio = b / a; + if ( + (smallA && (ratio * l2 < 0.1)) || + (!smallA && (l1 / ratio > 0.1)) + ) { + l3 = stdlib_base_expm1( ratio * stdlib_base_log1p( l2 ) ); + l3 = l1 + l3 + ( l3 * l1 ); + l3 = a * stdlib_base_log1p( l3 ); + result *= stdlib_base_exp( l3 ); + } + else { + l3 = stdlib_base_expm1( stdlib_base_log1p( l1 ) / ratio ); + l3 = l2 + l3 + ( l3 * l2 ); + l3 = b * stdlib_base_log1p( l3 ); + result *= stdlib_base_exp( l3 ); + } + } + else if ( stdlib_base_abs(l1) < stdlib_base_abs(l2) ) { + // First base near 1 only: + l = ( a * stdlib_base_log1p( l1 ) ) + ( b * stdlib_base_ln( ( y*cgh ) / bgh ) ); + if ( l <= STDLIB_CONSTANT_FLOAT64_MIN_LN || l >= STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + l += stdlib_base_ln(result); + if ( l >= STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + return 0.0 / 0.0; + } + result = stdlib_base_exp( l ); + } else { + result *= stdlib_base_exp( l ); + } + } + else { + // Second base near 1 only: + l = ( b * stdlib_base_log1p( l2 ) ) + ( a * stdlib_base_ln( (x*cgh) / agh ) ); + if ( l <= STDLIB_CONSTANT_FLOAT64_MIN_LN || l >= STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + l += stdlib_base_ln(result); + if ( l >= STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + return 0.0 / 0.0; + } + result = stdlib_base_exp( l ); + } else { + result *= stdlib_base_exp( l ); + } + } + } + else { + // General case: + b1 = (x * cgh) / agh; + b2 = (y * cgh) / bgh; + l1 = a * stdlib_base_ln(b1); + l2 = b * stdlib_base_ln(b2); + if ( + l1 >= STDLIB_CONSTANT_FLOAT64_MAX_LN || + l1 <= STDLIB_CONSTANT_FLOAT64_MIN_LN || + l2 >= STDLIB_CONSTANT_FLOAT64_MAX_LN || + l2 <= STDLIB_CONSTANT_FLOAT64_MIN_LN + ) { + // Oops, under/overflow, sidestep if we can: + if ( a < b ) { + p1 = stdlib_base_pow( b2, b / a ); + l3 = a * ( stdlib_base_ln(b1) + stdlib_base_ln(p1) ); + if ( l3 < STDLIB_CONSTANT_FLOAT64_MAX_LN && l3 > STDLIB_CONSTANT_FLOAT64_MIN_LN ) { + result *= stdlib_base_pow( p1 * b1, a ); + } else { + l2 += l1 + stdlib_base_ln(result); + if ( l2 >= STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + return 0.0 / 0.0; + } + result = stdlib_base_exp( l2 ); + } + } + else { + p1 = stdlib_base_pow( b1, a / b ); + l3 = ( stdlib_base_ln(p1) + stdlib_base_ln(b2) ) * b; + if ( l3 < STDLIB_CONSTANT_FLOAT64_MAX_LN && l3 > STDLIB_CONSTANT_FLOAT64_MIN_LN ) { + result *= stdlib_base_pow( p1 * b2, b ); + } else { + l2 += l1 + stdlib_base_ln( result ); + if ( l2 >= STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + return 0.0 / 0.0; + } + result = stdlib_base_exp( l2 ); + } + } + } + else { + // Finally the normal case: + result *= stdlib_base_pow( b1, a ) * stdlib_base_pow( b2, b ); + } + } + return result; +} + +/** +* Computes the partial derivative with respect to x of the incomplete beta function. +* +* @private +* @param x - input value (not equal to 0.0 or 1.0) +* @param a - first parameter +* @param b - second parameter +* @returns value of the partial derivative +*/ +static double ibetaDerivative( const double x, const double a, const double b ) { + double f1; + double y; + + f1 = ibetaPowerTerms( a, b, x, 1.0 - x, true ); + y = ( 1.0 - x ) * x; + f1 /= y; + return f1; +} + + +/** +* Evaluates the probability density function (PDF) for an F distribution with numerator degrees of freedom `d1` and denominator degrees of freedom `d2` at a value `x`. +* +* @param x - input value +* @param d1 - numerator degrees of freedom +* @param d2 - denominator degrees of freedom +* @returns evaluated PDF +* +* @example +* double y = pdf( 2.0, 0.5, 1.0 ); +* // returns ~0.057 +*/ +double stdlib_base_dists_f_pdf( const double x, const double d1, const double d2 ) { + double v1x; + double y; + double z; + if ( + stdlib_base_is_nan( x ) || + stdlib_base_is_nan( d1 ) || + stdlib_base_is_nan( d2 ) || + d1 <= 0.0 || + d2 <= 0.0 + ) { + return 0.0 / 0.0; + } + if ( x < 0.0 || x == STDLIB_CONSTANT_FLOAT64_PINF ) { + return 0.0; + } + if ( x == 0.0 ) { + if ( d1 < 2.0 ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + if ( d1 == 2.0 ) { + return 1.0; + } + return 0.0; + } + v1x = d1 * x; + if ( v1x > d2 ) { + y = ( d2 * d1 ) / ( ( d2 + v1x ) * ( d2 + v1x ) ); + return y * ibetaDerivative( d2 / ( d2+v1x ), d2/2.0, d1/2.0 ); + } + z = d2 + v1x; + y = ((z * d1) - (x * d1 * d1)) / ( z * z ); + return y * ibetaDerivative( v1x / ( d2+v1x ), d1/2.0, d2/2.0 ); +} diff --git a/lib/node_modules/@stdlib/stats/base/dists/f/pdf/test/test.native.js b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/test/test.native.js new file mode 100644 index 000000000000..1bbdb2b2f482 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/f/pdf/test/test.native.js @@ -0,0 +1,325 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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 resolve = require( 'path' ).resolve; +var tape = require( 'tape' ); +var tryRequire = require( '@stdlib/utils/try-require' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var abs = require( '@stdlib/math/base/special/abs' ); +var PINF = require( '@stdlib/constants/float64/pinf' ); +var NINF = require( '@stdlib/constants/float64/ninf' ); +var EPS = require( '@stdlib/constants/float64/eps' ); + + +// VARIABLES // + +var pdf = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( pdf instanceof Error ) +}; + + +// FIXTURES // + +var largeD1 = require( './fixtures/julia/large_d1.json' ); +var largeD2 = require( './fixtures/julia/large_d2.json' ); +var bothLarge = require( './fixtures/julia/both_large.json' ); +var cppLargeD1 = require( './fixtures/cpp/large_d1.json' ); +var cppLargeD2 = require( './fixtures/cpp/large_d2.json' ); +var cppBothLarge = require( './fixtures/cpp/both_large.json' ); + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof pdf, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'if provided `NaN` for any parameter, the function returns `NaN`', opts, function test( t ) { + var y = pdf( NaN, 0.0, 1.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + y = pdf( 0.0, NaN, 1.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + y = pdf( 0.0, 1.0, NaN ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided `+infinity` for `x` and a finite `d1` and `d2`, the function returns `0`', opts, function test( t ) { + var y = pdf( PINF, 1.0, 1.0 ); + t.strictEqual( y, 0.0, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided `-infinity` for `x` and a finite `d1` and `d2`, the function returns `0`', opts, function test( t ) { + var y = pdf( NINF, 1.0, 1.0 ); + t.strictEqual( y, 0.0, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided `d1 <= 0`, the function returns `NaN`', opts, function test( t ) { + var y; + + y = pdf( 2.0, -1.0, 2.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 0.0, -1.0, 2.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, NINF, 1.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, NINF, PINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, NINF, NINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, NINF, NaN ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided `d2 <= 0`, the function returns `NaN`', opts, function test( t ) { + var y; + + y = pdf( 2.0, 2.0, -1.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 0.0, 2.0, -1.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, 1.0, NINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, PINF, NINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, NINF, NINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = pdf( 2.0, NaN, NINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided `d1 < 2.0`, the function returns `+Infinity` when provided `0.0` for `x`', opts, function test( t ) { + var y; + + y = pdf( 0.0, 1.0, 1.0 ); + t.strictEqual( y, PINF, 'returns expected value' ); + + y = pdf( 0.0, 1.5, 1.0 ); + t.strictEqual( y, PINF, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided `d1 = 2.0`, the function returns `1.0` when provided `0.0` for `x`', opts, function test( t ) { + var y; + + y = pdf( 0.0, 2.0, 1.0 ); + t.strictEqual( y, 1.0, 'returns expected value' ); + + y = pdf( 0.0, 2.0, 3.0 ); + t.strictEqual( y, 1.0, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided `d1 > 2.0`, the function returns `0.0` when provided `0.0` for `x`', opts, function test( t ) { + var y; + + y = pdf( 0.0, 3.0, 1.0 ); + t.strictEqual( y, 0.0, 'returns expected value' ); + + y = pdf( 0.0, 3.5, 1.0 ); + t.strictEqual( y, 0.0, 'returns expected value' ); + + t.end(); +}); + +tape( 'the function evaluates the pdf for `x` given large `d1` and `d2` (tested against Julia)', opts, function test( t ) { + var expected; + var delta; + var tol; + var d1; + var d2; + var x; + var y; + var i; + + expected = bothLarge.expected; + x = bothLarge.x; + d1 = bothLarge.d1; + d2 = bothLarge.d2; + for ( i = 0; i < x.length; i++ ) { + y = pdf( x[i], d1[i], d2[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', d1: '+d1[i]+', d2: '+d2[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 2300.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. d1: '+d1[i]+'. d2: '+d2[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the pdf for `x` given large parameter `d1` (tested against Julia)', opts, function test( t ) { + var expected; + var delta; + var tol; + var d1; + var d2; + var x; + var y; + var i; + + expected = largeD1.expected; + x = largeD1.x; + d1 = largeD1.d1; + d2 = largeD1.d2; + for ( i = 0; i < x.length; i++ ) { + y = pdf( x[i], d1[i], d2[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', d1: '+d1[i]+', d2: '+d2[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 150.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. d1: '+d1[i]+'. d2: '+d2[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the pdf for `x` given large parameter `d2` (tested against Julia)', opts, function test( t ) { + var expected; + var delta; + var tol; + var d1; + var d2; + var x; + var y; + var i; + + expected = largeD2.expected; + x = largeD2.x; + d1 = largeD2.d1; + d2 = largeD2.d2; + for ( i = 0; i < x.length; i++ ) { + y = pdf( x[i], d1[i], d2[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', d1: '+d1[i]+', d2: '+d2[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 1800.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. d1: '+d1[i]+'. d2: '+d2[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the pdf for `x` given large `d1` and `d2` (tested against the Boost C++ library', opts, function test( t ) { + var expected; + var delta; + var tol; + var d1; + var d2; + var x; + var y; + var i; + + expected = cppBothLarge.expected; + x = cppBothLarge.x; + d1 = cppBothLarge.d1; + d2 = cppBothLarge.d2; + for ( i = 0; i < x.length; i++ ) { + y = pdf( x[i], d1[i], d2[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', d1: '+d1[i]+', d2: '+d2[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 150.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. d1: '+d1[i]+'. d2: '+d2[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the pdf for `x` given large parameter `d1` (tested against the Boost C++ library)', opts, function test( t ) { + var expected; + var delta; + var tol; + var d1; + var d2; + var x; + var y; + var i; + + expected = cppLargeD1.expected; + x = cppLargeD1.x; + d1 = cppLargeD1.d1; + d2 = cppLargeD1.d2; + for ( i = 0; i < x.length; i++ ) { + y = pdf( x[i], d1[i], d2[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', d1: '+d1[i]+', d2: '+d2[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 100.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. d1: '+d1[i]+'. d2: '+d2[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the pdf for `x` given large parameter `d2` (tested against the Boost C++ library)', opts, function test( t ) { + var expected; + var delta; + var tol; + var d1; + var d2; + var x; + var y; + var i; + + expected = cppLargeD2.expected; + x = cppLargeD2.x; + d1 = cppLargeD2.d1; + d2 = cppLargeD2.d2; + for ( i = 0; i < x.length; i++ ) { + y = pdf( x[i], d1[i], d2[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', d1: '+d1[i]+', d2: '+d2[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 1150.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. d1: '+d1[i]+'. d2: '+d2[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' ); + } + } + t.end(); +});