diff --git a/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/README.md b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/README.md index 0143acf5f34a..0907e401ea2b 100644 --- a/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/README.md +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/README.md @@ -2,7 +2,7 @@ @license Apache-2.0 -Copyright (c) 2018 The Stdlib Authors. +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. @@ -157,6 +157,105 @@ logEachMap( 'x: %d, r: %0.4f, p: %0.4f, ln(P(X=x;r,p)): %0.4f', x, r, p, logpmf + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/stats/base/dists/negative-binomial/logpmf.h" +``` + +#### stdlib_base_dists_negative_binomial_logpmf( x, r, p ) + +Evaluates the natural logarithm of the probability mass function (PMF) for a negative binomial distribution with number of successes until experiment is stopped `r` and success probability `p`. + +```c +double out = stdlib_base_dists_negative_binomial_logpmf( 5.0, 20.0, 0.8 ); +// returns ~-1.853 +``` + +The function accepts the following arguments: + +- **x**: `[in] double` input value. +- **r**: `[in] double` number of successes until experiment is stopped. +- **p**: `[in] double` success probability. + +```c +double stdlib_base_dists_negative_binomial_logpmf( const double x, const double r, const double p ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/stats/base/dists/negative-binomial/logpmf.h" +#include "stdlib/math/base/special/ceil.h" +#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 r; + double p; + double x; + double y; + int i; + + for ( i = 0; i < 10; i++ ) { + x = stdlib_base_ceil( random_uniform( 0.0, 30.0 ) ); + r = random_uniform( 0.0, 50.0 ); + p = random_uniform( 0.0, 1.0 ); + y = stdlib_base_dists_negative_binomial_logpmf( x, r, p ); + printf( "x: %lf, r: %lf, p: %lf, ln(P(X=x;r,p)): %lf\n", x, r, p, y ); + } +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/benchmark/benchmark.native.js new file mode 100644 index 000000000000..7001a45dfb29 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/benchmark/benchmark.native.js @@ -0,0 +1,72 @@ +/** +* @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 uniform = require( '@stdlib/random/array/uniform' ); +var discreteUniform = require( '@stdlib/random/array/discrete-uniform' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var format = require( '@stdlib/string/format' ); +var EPS = require( '@stdlib/constants/float64/eps' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var logpmf = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( logpmf instanceof Error ) +}; + + +// MAIN // + +bench( format( '%s::native', pkg ), opts, function benchmark( b ) { + var opts; + var r; + var p; + var x; + var y; + var i; + + opts = { + 'dtype': 'float64' + }; + x = discreteUniform( 100, 1, 100, opts ); + r = discreteUniform( 100, 1, 100, opts ); + p = uniform( 100, EPS, 1.0, opts ); + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + y = logpmf( x[ i % x.length ], r[ i % r.length ], p[ i % p.length ] ); + 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/negative-binomial/logpmf/benchmark/c/Makefile b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/benchmark/c/Makefile new file mode 100644 index 000000000000..979768abbcec --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/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/negative-binomial/logpmf/benchmark/c/benchmark.c b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/benchmark/c/benchmark.c new file mode 100644 index 000000000000..d6188624fb8e --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/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/negative-binomial/logpmf.h" +#include "stdlib/constants/float64/eps.h" +#include +#include +#include +#include +#include +#include + +#define NAME "negative-binomial-logpmf" +#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 r[ 100 ]; + double p[ 100 ]; + double x[ 100 ]; + double y; + double t; + int32_t i; + + for ( i = 0; i < 100; i++ ) { + x[ i ] = random_uniform( 1.0, 100.0 ); + r[ i ] = random_uniform( 1.0, 100.0 ); + p[ i ] = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 1.0 ); + } + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_dists_negative_binomial_logpmf( x[ i%100 ], r[ i%100 ], p[ 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/negative-binomial/logpmf/binding.gyp b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/binding.gyp new file mode 100644 index 000000000000..0d6508a12e99 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/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/negative-binomial/logpmf/examples/c/Makefile b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/examples/c/Makefile new file mode 100644 index 000000000000..c8f8e9a1517b --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/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/negative-binomial/logpmf/examples/c/example.c b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/examples/c/example.c new file mode 100644 index 000000000000..76f26c0bee84 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/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/negative-binomial/logpmf.h" +#include "stdlib/math/base/special/ceil.h" +#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 r; + double p; + double x; + double y; + int i; + + for ( i = 0; i < 10; i++ ) { + x = stdlib_base_ceil( random_uniform( 0.0, 30.0 ) ); + r = random_uniform( 0.0, 50.0 ); + p = random_uniform( 0.0, 1.0 ); + y = stdlib_base_dists_negative_binomial_logpmf( x, r, p ); + printf( "x: %lf, r: %lf, p: %lf, ln(P(X=x;r,p)): %lf\n", x, r, p, y ); + } +} diff --git a/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/include.gypi b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/include.gypi new file mode 100644 index 000000000000..bee8d41a2caf --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/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 {NonNegativeNumber} a - function parameter +* @param {NonNegativeNumber} b - function parameter +* @param {Probability} x - function parameter +* @param {Probability} y - probability equal to `1-x` +* @param {boolean} normalized - boolean indicating whether to evaluate the power terms of the regularized or non-regularized incomplete beta function +* @returns {number} 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 *= sqrt( bgh / STDLIB_CONSTANT_FLOAT64_E ); + result *= 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 (0 < x <= 1) +* @param a - first parameter +* @param b - second parameter (must be greater than 1) +* @returns value of the partial derivative +*/ +static double ibetaDerivative( const double x, const double a, const double b ) { + double f1; + double y; + if ( x == 1.0 ) { + return 0.0; + } + // Regular cases: + f1 = ibetaPowerTerms( a, b, x, 1.0 - x, true ); + y = ( 1.0 - x ) * x; + f1 /= y; + return f1; +} + +/** +* Evaluates the natural logarithm of the probability mass function (PMF) for a negative binomial distribution with number of successes until experiment is stopped `r` and success probability `p`. +* +* @param x - input value +* @param r - number of successes until experiment is stopped +* @param p - success probability +* @returns evaluated logPMF +* +* @example +* double y = logpmf( 5.0, 20.0, 0.8 ); +* // returns ~-1.853 +*/ +double stdlib_base_dists_negative_binomial_logpmf( const double x, const double r, const double p ) { + if ( + stdlib_base_is_nan( x ) || + stdlib_base_is_nan( r ) || + stdlib_base_is_nan( p ) || + r <= 0.0 || + p <= 0.0 || + p > 1.0 + ) { + return 0.0 / 0.0; + } + if ( !stdlib_base_is_nonnegative_integer( x ) || p == 0.0 ) { + return STDLIB_CONSTANT_FLOAT64_NINF; + } + return stdlib_base_ln( p ) - stdlib_base_ln( r + x ) + stdlib_base_ln( ibetaDerivative( p, r, x + 1.0 ) ); +} diff --git a/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/test/test.native.js b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/test/test.native.js new file mode 100644 index 000000000000..ba31ce34985c --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/negative-binomial/logpmf/test/test.native.js @@ -0,0 +1,238 @@ +/** +* @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 logpmf = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( logpmf instanceof Error ) +}; + + +// FIXTURES // + +var smallSmall = require( './fixtures/julia/small_small.json' ); +var smallHigh = require( './fixtures/julia/small_high.json' ); +var highSmall = require( './fixtures/julia/high_small.json' ); +var highHigh = require( './fixtures/julia/high_high.json' ); + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof logpmf, '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 = logpmf( NaN, 20.0, 0.5 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + y = logpmf( 0.0, NaN, 0.5 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + y = logpmf( 0.0, 20.0, NaN ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided a negative integer for `x` and a valid `r` and `p`, the function returns `-Infinity`', opts, function test( t ) { + var y = logpmf( NINF, 20, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + y = logpmf( -20.0, 20, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + y = logpmf( -100.0, 20, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + y = logpmf( -1.0, 20, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided a non-integer for `x` and a valid `r` and `p`, the function returns `-Infinity`', opts, function test( t ) { + var y = logpmf( -1.5, 20.0, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + y = logpmf( -0.5, 20.0, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + y = logpmf( 1.5, 20.0, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + y = logpmf( 2.5, 20.0, 0.5 ); + t.strictEqual( y, NINF, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided `r` which is not a positive number, the function returns `NaN`', opts, function test( t ) { + var y; + + y = logpmf( 2.0, -0.5, 0.5 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = logpmf( 0.0, -1.0, 0.5 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = logpmf( 2.0, NINF, 0.5 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + t.end(); +}); + +tape( 'if provided a success probability `p` outside of `[0,1]`, the function returns `NaN`', opts, function test( t ) { + var y; + + y = logpmf( 2.0, 20, -1.0 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = logpmf( 0.0, 20, 1.5 ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = logpmf( 2.0, 20, NINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + y = logpmf( 2.0, 20, PINF ); + t.strictEqual( isnan( y ), true, 'returns expected value' ); + + t.end(); +}); + +tape( 'the function evaluates the logpmf for `x` given large `r` and `p`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var r; + var p; + var x; + var y; + + expected = highHigh.expected; + x = highHigh.x; + r = highHigh.r; + p = highHigh.p; + for ( i = 0; i < x.length; i++ ) { + y = logpmf( x[i], r[i], p[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', r: '+r[i]+', p: '+p[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[i] ); + tol = 80.0 * EPS * abs( expected[i] ); + t.strictEqual( delta <= tol, true, 'within tolerance. x: '+x[i]+'. r: '+r[i]+'. p: '+p[i]+'. y: '+y+'. Expected: '+expected[i]+'. tol: '+tol+'. delta: '+delta+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the logpmf for `x` given large parameter `r` and small `p`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var r; + var p; + var x; + var y; + + expected = highSmall.expected; + x = highSmall.x; + r = highSmall.r; + p = highSmall.p; + for ( i = 0; i < x.length; i++ ) { + y = logpmf( x[i], r[i], p[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', r: '+r[i]+', p: '+p[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[i] ); + tol = 10.0 * EPS * abs( expected[i] ); + t.strictEqual( delta <= tol, true, 'within tolerance. x: '+x[i]+'. r: '+r[i]+'. p: '+p[i]+'. y: '+y+'. Expected: '+expected[i]+'. tol: '+tol+'. delta: '+delta+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the logpmf for `x` given small `r` and large `p`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var r; + var p; + var x; + var y; + + expected = smallHigh.expected; + x = smallHigh.x; + r = smallHigh.r; + p = smallHigh.p; + for ( i = 0; i < x.length; i++ ) { + y = logpmf( x[i], r[i], p[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', r: '+r[i]+', p: '+p[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[i] ); + tol = 45.0 * EPS * abs( expected[i] ); + t.strictEqual( delta <= tol, true, 'within tolerance. x: '+x[i]+'. r: '+r[i]+'. p: '+p[i]+'. y: '+y+'. Expected: '+expected[i]+'. tol: '+tol+'. delta: '+delta+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the logpmf for `x` given small `r` and `p`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var r; + var p; + var x; + var y; + + expected = smallSmall.expected; + x = smallSmall.x; + r = smallSmall.r; + p = smallSmall.p; + for ( i = 0; i < x.length; i++ ) { + y = logpmf( x[i], r[i], p[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+', r: '+r[i]+', p: '+p[i]+', y: '+y+', expected: '+expected[i] ); + } else { + delta = abs( y - expected[i] ); + tol = 10.0 * EPS * abs( expected[i] ); + t.strictEqual( delta <= tol, true, 'within tolerance. x: '+x[i]+'. r: '+r[i]+'. p: '+p[i]+'. y: '+y+'. Expected: '+expected[i]+'. tol: '+tol+'. delta: '+delta+'.' ); + } + } + t.end(); +});