README

libRmath.js

This is a 100% Pure Javascript ( TypeScript ) re-write of Statistical R nmath "core" numerical library found here. This is a manual re-write, "emscripten" was not used.

Summary

libRmath.js port contains all functions implemented in R nmath counterpart:

probability and quantile functions related to 21 distributions.
functions to work with Special functions in mathematics (Bessel,Gamma,Beta).
7 uniform PRNG's. (same output pseudo-random sequence in R for the same initial seeding).
6 normally distributed PRNG's. (same output sequence in R for te same initial seeding).
Function vector (array) arguments follow the R recycling rule.

With this library it becomes trivial to implement hypothesis testing in javascript, calculating p-values and (1 - α) confidence intervals. (ANOVA uses the F-distribution. Tukey HSD uses ptukey function, etc, etc).

All functions in libRmath.so has been re-written to Javascript (Typescript). Use the library with either vanilla Javascript or Typescript. Type definition files are included in the npm package.

Node and Web

The library is an UMD library, it can be used in a web client as in server side node environment.

Installation

serverside

npm install --save lib-r-math.js

web

The module directory contains a minimized bundle for use in html <script> tag. The library is attached to the window.libR object after loading.

<!-- script src="your_server_url/libR.min.js"></script -->
<!-- this example uses unpkg as CDN -->
<script src="https://unpkg.com/lib-r-math.js@1.0.74/dist/lib/libR.min.js">
<script>
  const libR = window.libR;
  //fetch some distribution namespaces
  const { Tukey, Normal, Beta, StudentT, Wilcoxon } = libR;
</script>

Differences with R
Read this first: Helper functions
Uniform Pseudo Random Number Generators
Normal Random Number Generators
Normal and Uniform distributions
- Uniform distribution
- Normal distribution
Other Probability Distributions
Special Functions of Mathematics

Differences with R

Some implementation differences exist with R nmath

PRNG's are not global singletons, but separate object instances and you can have as many as you want. The programmer has the choice of having different deviate generators sharing a common source PRNG.
Wilcoxon Sum Rank functions dwilcox, pwilcox, qwilcox use a fraction of the memory, (R will give memory allocation errors for samples ~1000). The JS solution allocates memory sparsely.

Helper functions for porting `R` programs

Summary

R language operators and function arguments can work with vectorized input. These helper functions are used to mimic this functionality and assist porting of scripts from the R ecosystem using libRmath.js.

`div`

Divides scalar or an array of values with element of the second array or scalar.

Usage:

const libR = require('lib-r-math.js');
const { div } = libR.R;

//1
div(3, 5); //= 3/5
//0.6

div([0, 1, 2, 3], 5);
//[0, 0.2, 0.4, 0.6]

div([10,2,3],[2,4]);// Uses R recycling rules
//[ 5, 0.5, 1.5 ]

`mult`

Multiplies scalar or an array of values with another scalar or array of values. Applies R recycling rules in case of arguments of unequal array length.

Usage:

const libR = require('lib-r-math.js');
const { mult } = libR.R;

//1
mult(3, 5); //= 3*5
//15

mult([0, 1, 2, 3], [5,2]); // R recycling rules apply
//[ 0, 2, 10, 6 ]

`asArray`

Creates a new function from an existing one for it to always return its result as an array.

Usage:


const libR = require('lib-r-math.js');
const { asArray } = libR.R;

const r = asArray(Math.random);

//always returns the result wrapped in an array
r()
//[ 0.39783583929513 ]
r()
//[ 0.04431401890179831 ]
r()
//[ 0.7629304997301447 ]

`sum`

Analog to R's sum function. Calculates the sum of all elements of an array.

const libR = require('lib-r-math.js');
const { sum } = libR.R;

//1
sum(5); //trivial
//5

//2
sum([1, 2, 3, 4]);
//10

`summary`

Gives summary information of numeric data in an array.

typescript decl

declare function summary(data: number[]): ISummary;

//output
interface ISummary {
  N: number; // number of samples in "data"
  mu: number; // mean of "data"
  population: {
    variance: number; // population variance (data is seen as finite population)
    sd: number; // square root of the population variance
  };
  sample: {
    variance: number; // sample variance (data is seen as a small sample from an very large population)
    sd: number; // square root of "sample variance"
  };
  relX; // = x-E(x)
  relX2; // = ( x-E(x) )^2
  stats: {
    min: number; // minimal value from "data"
    '1st Qu.': number; // 1st quantile from "data"
    median: number; // median value from "data
    '3rd Qu.': number; // 3rd quantile from "data"
    max: number; // maximum value in data
  };
}

Usage:

const libR = require('lib-r-math.js');
const { summary } = libR.R;

summary([360, 352, 294, 160, 146, 142, 318, 200, 142, 116])
/*
{ N: 10,
  mu: 223,
  population: { variance: 8447.4, sd: 91.90973833060346 },
  sample: { variance: 9386, sd: 96.88137075826292 },
  relX: [ 137, 129, 71, -63, -77, -81, 95, -23, -81, -107 ],
  relX2: [ 18769, 16641, 5041, 3969, 5929, 6561, 9025, 529, 6561, 11449 ],
  stats: { min: 116, '1st Qu.': 143, median: 180, '3rd Qu.': 312, max: 360 } }
*/

`numberPrecision`

Truncates numbers to a specified significant digits. Takes single numeric value as argument or an array of numbers.

Usage:

const libR = require('lib-r-math.js');

const digits4 = libR.R.numberPrecision(4);

//1 single value
const pr4a = digits4(1.12345678);
//1.123

//2 works with arrays
const pr4b = digits4([0.4553, -2.1243]);
//[ 0.4553, -2.124 ]

`any`

Test a Predicate for each element in an Array. Returns true or false depending on a test function.

Usage:

const libR = require('lib-r-math.js');
const any = libR.R.any;

//1
any([1, 2, 3, 4])(x => x < 2);
//true

//2
any([1, 2, 3, 4])(x => x > 5);
//false

`arrayrify` (DEPRICATED use `multiplex`)

Mimics R vectorized function arguments. Wraps an existing function changing the first first argument to accept both scalar (number) or an array( number[] ).

Note: Only the first argument is vectorized

typescript decl

declare function arrayrify<T, R>(fn: (x: T, ...rest: any[]) => R);

R example

# raise each element by power of 2
c(1,2,3,4)^2
#[1]  1  4  9 16

Javascript equivalent

 const libR = require('lib-r-math.js');
 const { arrayrify } = libR.R;

 // create vectorize "/" operator
 const pow = arrayrify(Math.pow);

 pow(3, 4); // 81

 pow([3, 4, 5], 4); //81 256 625

`each`

Functional analog to Array.prototype.forEach, but can also loop over object properties. The return type can be either an new array or a scalar (see Example).

Example:

const libR = require('lib-r-math.js');
const { each } = libR.R;

each(11)(v => console.log(v * 2)) ;

// single element array result are forced to return scalar
each([3])(v => console.log(v * 2));

each([11, 12])( (v, idx) => console.log({ v, idx}));

//looping over object properties
each({ p:1, name:'myname' })( (value, key) => console.log(`${key}=${value}`))

`flatten` or `c` (alias)

Analog to R's c function. Constructs a final array by (recursively) flattening and merging all of its arguments which can be a combination of scalars and arrays.

Example:

const libR = require('lib-r-math.js');

// optionally rename as `c` to make it look like `R`
const { c } = libR.R;

c(-1, 0, [1], 'r', 'b', [2, 3, [4, 5]]);
// [ -1, 0, 1, 'r', 'b', 2, 3, 4, 5 ]

`map`

Functional analog to Array.prototype.map, but can also loop over object properties. The return type can be either an new array or a scalar (see Example).

Example:

const libR = require('lib-r-math.js');
const { map } = libR.R;

map(11)(v => v * 2);
//22

// single element array result are forced to return scalar
map([3])(v => v * 2);
//6

map([11, 12])( (v, idx) => idx);
// [0, 1]

//looping over object properties
map({ p:1, name:'myname' })( (value, key) => `${key}=${value}`)
//["p=1", "name=myname"]

`selector`

Filter function generator, to be used with Array.prototype.filter to pick elements based on their order (zero based index) in the array. Usually used together with seq to pick items from an array.

NOTE: Always returns an instance of Array.

Example:

const libR = require('lib-r-math.js');
const { selector } = libR.R;

['an', 'array', 'with', 'some', 'elements'].filter(
  selector([0, 2, 3]) // select values at these indexes
);
//[ 'an', 'with', 'some']

['an', 'array', 'with', 'some', 'elements'].filter(
  selector(3) // just one value at position 3
);
//['some']
const seq = libR.R.seq()(); // see "seq" for defaults.

[7, 1, 2, 9, 4, 8, 16].filter(
  selector(
    seq(0, 6, 2) // creates an array [ 0, 2, 4, 6]
  )
);
// returns [7, 2, 4, 16]

`seq`

typescript decl

const seq = (adjustUp = 0) => (adjustDown = adjust) => (
  start: number,
  end: number = 1,
  step: number = 1
) => number[];

R analog to the seq function in R. Generates an array between start and end (inclusive) using step (defaults to 1). The JS implementation ignores the sign of the step argument and only looks at its absolute value.

Like in R, If (end-start)/step is not an exact integer, seq will not overstep the bounds while counting up (or down).

Arguments:

adjustUp: (default 0). If end >= start then adjust value is added to every element in the array.
adjustDown: (default 0). If start >= end then adjustMin value is added to every element in the array.
start: (inclusive) the sequence start value
stop: defaults to 1. (inclusive) the sequence stop value if possible with step
step: defaults to 1, sign is ignored. Sign is inferred from the values of start and stop.

First we look how seq works in R.

seq(1,3,0.5);
#[1] 1.0 1.5 2.0 2.5 3.0

seq(7,-2, -1.3);
#[1]  7.0  5.7  4.4  3.1  1.8  0.5 -0.8

Equivalent in Javascript

const libR = require('lib-r-math.js');

// seqA is a sequence generator
let seqA = libR.R.seq()();

seqA(1, 5);
//[ 1, 2, 3, 4, 5 ]

seqA(5, -3);
//[ 5, 4, 3, 2, 1, 0, -1, -2, -3 ]

seqA(3)
//[3, 2, 1]

//add 1 if stepping upwards, add -2 if stepping downwards
let seqB = libR.R.seq(1)(-2);

seqB(0, 4); //range will be adjusted with '1'
//[ 1, 2, 3, 4]
seqB(6, 5, 0.3); //range will be adjusted with '-2', step
//[4, 3.7, 3.4, 3.1]

`multiplex`

Turns an existing javascript function into one that follows the R argument recycling rule.

Multiplexes the value of several array arguments into one array with the use of a mapping function.

The length of the result is the maximum of the lengths of the parameters. All parameters are recycled to that length.

const libR = require('lib-r-math.js');

const { multiplex, c } = libR.R;

//make the build in Math function follow R-recycling rule
const pow = multiplex(Math.pow);
//
pow([1, 2, 3, 4], 2); //squared
//[ 1, 4, 9, 16 ]

//powers of 2
pow(2, [2, 3, 4])
    //[ 4, 8, 16 ]

//R recycling rule
pow([2, 3], [2, 3, 4, 5]);
//[4, 27, 16, 243]
//4 = 2 ^ 2
//27 = 3 ^ 3
//16 = 2 ^ 4
//243 = 3 ^ 5

`timeseed`

Generates a number based by on your system time clock. Intended use is with PRNG (re)initialization. Its a synchronized function that will wait for some milliseconds before sampling the system clock and returning a result.

Usage:

const libR = require('lib-r-math.js');

const { rng: { timeseed } } = libR;

timeseed();
//2632999169 , based on timestamp

Uniform Pseudo Random Number Generators

Summary

In 'R', the functions that generate random deviates of distributions (Example: Poisson (rpois), Student-t (rt), Normal (rnorm), etc) use uniform PRNG's directly or indirectly (as wrapped in a normal distributed PRNG). This section discusses the uniform distributed PRNG's that have been ported from R to JS.

The 7 Uniform Random Number Generators

All 7 uniform random generators have been ported and tested to yield exactly the same as their R counterpart.

Improvements compared to R

In R it is impossible to use different types of uniform random generators at the same time because of a global shared seed buffer. In our port every random generator has its own buffer and can therefore be used concurrently.

General Usage

All uniform random generator classes have these public methods:

init: set the random generator seed. Same as R set.seed()
seed (read/write property): get/set the current seed values as an array. Same as R .Random.seed.
unif_random: get a random value, same as runif(1) in R

"Mersenne Twister"

From Matsumoto and Nishimura (1998). A twisted GFSR with period 2^19937 - 1 and equi-distribution in 623 consecutive dimensions (over the whole period). The seed is a 624-dimensional set of 32-bit integers plus a current position in that set.

usage example:

const libR = require('lib-r-math.js');
const {
    R: { numberPrecision },
    rng: { MersenneTwister, timeseed }
} = libR;

//helpers
const precision = numberPrecision(9); //9 digits accuracy

//example
const mt = new MersenneTwister(12345); // initialize with seed = 12345

//example
mt.init(timeseed()); // Use seed derived from system clock

//example
mt.init(0); // re-initialize with seed = 0

// show first 8 values of the seed buffer of the mt instance.
mt.seed.slice(0, 8);
/*[ 624,   1280795612,  -169270483,  -442010614,  -603558397,  -222347416,
  1489374793, 865871222 ]
*/

const rmt1 = mt.unif_rand(5);
precision(rmt1);
//[ 0.8966972, 0.265508663, 0.3721239, 0.572853363, 0.90820779 ]

Equivalent in R

RNGkind("Mersenne-Twister")
set.seed(0)

#show first 8 values of the seed buffer
.Random.seed[2:9]
#[1]        624 1280795612 -169270483 -442010614 -603558397 -222347416 1489374793
#[8]  865871222

runif(5)
#[1] 0.8966972 0.2655087 0.3721239 0.5728534
#[5] 0.9082078

"Wichmann-Hill"

The seed, is an integer vector of length 3, where each element is in 1:(p[i] - 1), where p is the length 3 vector of primes, p = (30269, 30307, 30323). The Wichmann–Hill generator has a cycle length of 6.9536e12 = ( 30269 * 30307 * 30323 ), see Applied Statistics (1984) 33, 123 which corrects the original article).

usage example:

const libR = require('lib-r-math.js');
const {
  rng: { WichmannHill, timeseed },
  R: { numberPrecision }
} = libR;

// some helpers

const precision = numberPrecision(9);

// Some options on seeding given below
const wh = new WichmannHill(1234); // initialize seed with 1234 on creation (default 0)
//
wh.init(timeseed()); // re-init seed with a random seed based on timestamp

wh.init(0); // re-init seed to zero
wh.seed; // show seed
//[ 2882, 21792, 10079 ]
const rwh1 = wh.unif_rand(5);
precision(rwh1);
//[ 0.462553151, 0.26582675, 0.57721078, 0.510793206, 0.337560559 ]

in R console:

> RNGkind("Wichmann-Hill")
> set.seed(0)
> runif(5)
[1] 0.4625532 0.2658268 0.5772108 0.5107932
[5] 0.3375606

"Marsaglia-Multicarry"

A multiply-with-carry RNG is used, as recommended by George Marsaglia in his post to the mailing list ‘sci.stat.math’. It has a period of more than 2^60 and has passed all tests (according to Marsaglia). The seed is two integers (all values allowed).

usage example:

const libR = require('lib-r-math.js');
const {
    rng: { MarsagliaMultiCarry, timeseed },
    R: { numberPrecision }
} = libR;

//usefull helpers
const precision = numberPrecision(9); //9 significant digits

// Some options on seeding given below
const mmc = new MarsagliaMultiCarry(1234); // use seed = 1234 on creation

mmc.init(timeseed());
mmc.init(0); // also, defaults to '0' if seed is not specified
mmc.seed;
//[ -835792825, 1280795612 ]

const rmmc = mmc.unif_rand(5);
precision(rmmc);
//[ 0.169153755, 0.53154353, 0.594605297, 0.233315406, 0.45765618 ]

in R console:

> RNGkind("Marsaglia-Multicarry")
> set.seed(0)
# we cannot access the PRNG directly
# we need to use runif wrapper.
> runif(5)
[1] 0.1691538 0.5315435 0.5946053 0.2333154
[5] 0.4576562

"Super Duper"

Marsaglia's famous Super-Duper from the 70's. This is the original version which does not pass the MTUPLE test of the Diehard battery. It has a period of about 4.6*10^18 for most initial seeds. The seed is two integers (all values allowed for the first seed: the second must be odd).

We use the implementation by Reeds et al (1982–84).

usage example:

const libR = require('lib-r-math.js');
const {
    rng: { SuperDuper, timeseed },
    R: { seq, numberPrecision }
} = libR;

//usefull helpers
const precision = numberPrecision(9); //9 significant digits

// Seeding possibilities shown below
const sd = new SuperDuper(1234); // use seed = 1234 on creation
sd.init(timeseed()); // re-initialize with random seed based on timestamp
sd.init(0); // re-initialize with seed = 0.
//
sd.seed;
//[ -835792825, 1280795613 ]

const rsd1 = sd.unif_rand(5);
precision(rsd1);
//[ 0.640403562, 0.592731255, 0.412968712, 0.187729399, 0.267905811 ]

in R console:

> RNGkind("Super-Duper")
> set.seed(0)
> runif(5)
[1] 0.6404036 0.5927313 0.4129687 0.1877294
[5] 0.2679058

"Knuth TAOCP"

An earlier version from Knuth (1997).

The 2002 version was not backwards compatible with the earlier version: the initialization of the GFSR from the seed was altered. R did not allow you to choose consecutive seeds, the reported ‘weakness’, and already scrambled the seeds.

usage example:

const libR = require('lib-r-math.js');
const {
    rng: { KnuthTAOCP, timeseed },
    R: { seq, numberPrecision }
} = libR;

//usefull helpers
const precision = numberPrecision(9); //9 significant digits

// Seeding possibilities shown below
const kn97 = new KnuthTAOCP(1234); // use seed = 1234 on creation
kn97.init(timeseed()); // re-initialize with random seed based on timestamp
kn97.init(0); // re-initialize with seed = 0.

kn97.seed;
// 101 unsigned integer array, only shown the first few values
/*[ 673666444,
  380305043,
  1062889978,
  926003693,
 .
 .]*/

const rkn97 = kn97.unif_rand(5);
// limit precision to 9 digits
precision(rkn97);
//[ 0.627400767, 0.354186672, 0.989893431, 0.862408143, 0.662299205 ]

in R console:

> RNGkind("Knuth-TAOCP")
> set.seed(0)
> runif(5)
[1] 0.6274008 0.3541867 0.9898934 0.8624081 0.6622992

"Knuth TAOCP 2002"

A 32-bit integer GFSR using lagged Fibonacci sequences with subtraction. That is, the recurrence used is

X[j] = (X[j-100] - X[j-37]) mod 2^30

and the ‘seed’ is the set of the 100 last numbers (actually recorded as 101 numbers, the last being a cyclic shift of the buffer). The period is around 2¹²⁹.

usage example:

const libR = require('lib-r-math.js');
const {
    rng: { KnuthTAOCP2002, timeseed },
    R: { numberPrecision }
} = libR;

//some helpers
const precision = numberPrecision(9);

// Seeding possibilities shown below
const kt2002 = new KnuthTAOCP2002(1234); // use seed = 1234 on creation

kt2002.init(timeseed()); //re-initialize with random seed based on timestamp

kt2002.init(0); //re-initialize with seed = 0
kt2002.seed;
// 101 unsigned integer array
//[ 481970911,
//  634898052,
//  994481106,
//  607894626,
//  1044251579,
//  763229919,
//  638368738,
// .
// .]
const kt02 = kt2002.unif_rand(5);
precision(kt02);
//[ 0.195819038, 0.753866884, 0.472411247, 0.193160437, 0.19501841 ]

in R console:

> RNGkind("Knuth-TAOCP-2002")
> set.seed(0)
> runif(5)
[1] 0.1958190 0.7538669 0.4724112 0.1931604
[5] 0.1950184

"L'Ecuyer-CMRG":

A ‘combined multiple-recursive generator’ from L'Ecuyer (1999), each element of which is a feedback multiplicative generator with three integer elements: thus the seed is a (signed) integer vector of length 6. The period is around 2^191.

The 6 elements of the seed are internally regarded as 32-bit unsigned integers. Neither the first three nor the last three should be all zero, and they are limited to less than 4294967087 and 4294944443 respectively.

This is not particularly interesting of itself, but provides the basis for the multiple streams used in package parallel.

usage example:

const libR = require('lib-r-math.js');
const {
    rng: { LecuyerCMRG, timeseed },
    R: { numberPrecision }
} = libR;

//some helpers
const precision = numberPrecision(9);

// Seeding possibilities shown below
const lc = new LecuyerCMRG(1234);

lc.init(timeseed()); //re-initialize with random seed based on timestamp

lc.init(0); //re-initialize with seed = 0
lc.seed;
/*
[ -835792825,
  1280795612,
  -169270483,
  -442010614,
  -603558397,
  -222347416 ]
*/
const lc1 = lc.unif_rand(5);
precision(lc1);
//[ 0.332927492, 0.890352618, 0.163963441, 0.299050824, 0.395239092 ]

in R console:

> RNGkind("L'Ecuyer-CMRG")
> set.seed(0)
> .Random.seed[2:7]  #show the seeds
[1] -835792825 1280795612 -169270483 -442010614
[5] -603558397 -222347416

#pick 6 random numbers
#same numbers as generated in javascript
> runif(5)
[1] 0.3329275 0.8903526 0.1639634 0.2990508
[5] 0.3952391

Normal distributed Random Number Generators

Summary

This section discusses the normal distributed PRNG's that have been ported from R to JS.

All 6 normal random generators have been ported and tested to yield exactly the same as their R counterpart.

General Use

All normal random generator adhere to the same principles:

A constructor that takes an instance of a uniform PRNG as an argument
A function norm_random: get a random value, same as rnorm(1) in R.
A function unif_random: the underlying uniform PRNG.
The default argument for the constructor for normal PRNG is : Mersenne-Twister.
The class instance property rng contains the wrapped uniform PRNG instance.
All PRNG producing normal variates are packaged under the JS library name space rng.normal.

"Ahrens Dieter"

Ahrens, J. H. and Dieter, U. (1973) Extensions of Forsythe's method for random sampling from the normal distribution. Mathematics of Computation 27, 927-937.

example usage:

const libR = require('lib-r-math.js');
const {
    rng: {
        SuperDuper,
        normal: { AhrensDieter }
    },
    R: { numberPrecision }
} = libR;

//helper
const precision = numberPrecision(9);

// explicit use of uniform PRNG
const sd = new SuperDuper(0);
const ad1 = new AhrensDieter(sd);

// At any time reset normal PRNG seed via wrapped uniform.
sd.init(9987);
ad1.rng.init(9987); //same as above

// uses default uniform PRNG "MersenneTwister" with seed 0
const ad2 = new AhrensDieter();

// reference to uniform PRNG under rng property
ad2.rng.init(0);

// bleed the normal PRNG
const rngAd = ad2.norm_rand(5);
precision(rngAd);
//[ -1.17616753, 0.674117732, 1.06414352, -0.143897298, -1.2311498 ]

// its still possible to bleed the uniform PRNG from the normal PRNG
ad2.unif_rand();
//0.2016819310374558

ad2.rng.unif_rand();
//0.8983896849676967

in R console

> RNGkind("Mersenne-Twister",normal.kind="Ahrens-Dieter")
> set.seed(0)
> rnorm(5)
[1] -1.1761675  0.6741177  1.0641435 -0.1438973
[5] -1.2311498
> runif(2)
[1] 0.2016819 0.8983897

Box Muller

Box, G. E. P. and Muller, M. E. (1958) A note on the generation of normal random deviates. Annals of Mathematical Statistics 29, 610–611.

Example usage:

const libR = require('lib-r-math.js');

const {
    rng: {
        SuperDuper,
        normal: { BoxMuller }
    },
    R: { numberPrecision }
} = libR;

// helper
const precision = numberPrecision(9);

const sd = new SuperDuper(0);
const bm1 = new BoxMuller(sd);

// At any time reset normal PRNG seed, with the reference to uniform PRNG
sd.init(0);

// uses default uniform PRNG: MersenneTwister with seed 0
const bm2 = new BoxMuller();

// reference to uniform PRNG under rng property
bm2.rng.init(0);
// bleed the normal PRNG
const bm = bm2.norm_rand(5);
precision(bm);
//[ 1.29738758, -0.984378527, -0.732798867, 0.759774198, 1.49998876 ]

// its possible to bleed the uniform PRNG from the normal PRNG
bm2.unif_rand();
//0.8983896849676967

bm2.rng.unif_rand();
//0.944675268605351

R equivalent

> RNGkind("Mersenne-Twister",normal.kind="Box-Muller")
> set.seed(0)
> rnorm(5)
[1]  1.2973876 -0.9843785 -0.7327989  0.7597742
[5]  1.4999888
> runif(2)
[1] 0.8983897 0.9446753

Buggy Kinderman Ramage

Kinderman, A. J. and Ramage, J. G. (1976) Computer generation of normal random variables. Journal of the American Statistical Association 71, 893-896.

The Kinderman-Ramage generator used in versions prior to 1.7.0 (now called "Buggy") had several approximation errors and should only be used for reproduction of old results.

example usage:

const libR = require('lib-r-math.js');
const {
    R: { numberPrecision },
    rng: { SuperDuper, normal: { BuggyKindermanRamage } }
} = libR;
//helper
const precision = numberPrecision(9);

// Possible to arbitraty uniform PRNG source (example: SuperDuper)
const sd = new SuperDuper(0);
const bkm1 = new BuggyKindermanRamage(sd);

// At any time reset normal PRNG seed, with the reference to uniform PRNG
sd.init(0);
bkm1.rng.init(0); //same as above

// uses default uniform PRNG new MersenneTwister with seed 0
const bkm2 = new BuggyKindermanRamage();

// reference to uniform PRNG under rng property
bkm2.rng.init(0);

// bleed the normal PRNG
const bk = bkm2.norm_rand(5);
precision(bk);
//[ 0.3216151, 1.23251561, 0.280369528, -1.17519641, -1.60471361 ]

// its possible to bleed the uniform PRNG from the normal PRNG
bkm2.unif_rand();
//0.17655675252899528

bkm2.rng.unif_rand();
//0.6870228466577828

in R Console

RNGkind("Mersenne-Twister",normal.kind="Buggy Kinderman-Ramage");
set.seed(0);
rnorm(5);
#[1]  0.3216151  1.2325156  0.2803695 -1.1751964
#[5] -1.6047136

runif(2);
#[1] 0.1765568 0.6870228

Inversion

Inverse transform sampling wiki

example usage:

const libR = require('lib-r-math.js');
// Possible to arbitraty uniform PRNG source (example: SuperDuper)
const {
    rng: { SuperDuper, normal: { Inversion } },
    R: { numberPrecision }
} = libR;
//helper
const precision = numberPrecision(9);

const sd = new SuperDuper(0);
const inv1 = new Inversion(sd);
// At any time reset normal PRNG seed, with the reference to uniform PRNG
sd.init(0);

// uses as default uniform PRNG "MersenneTwister" with seed 0;
const inv2 = new Inversion();

// reference to uniform PRNG under rng property
inv2.rng.init(0);

// bleed the normal PRNG
const in2 = inv2.norm_rand(5);
precision(in2);
//[ 1.26295428, -0.326233361, 1.32979926, 1.27242932, 0.414641434 ]

// its possible to bleed the uniform PRNG from the normal PRNG
inv2.unif_rand();
//0.061786270467564464

inv2.rng.unif_rand();
//0.20597457489930093

in R console

> RNGkind("Mersenne-Twister",normal.kind="Inversion")
> set.seed(0)
> rnorm(5)
[1]  1.2629543 -0.3262334  1.3297993  1.2724293
[5]  0.4146414
> runif(2)
[1] 0.06178627 0.20597457

Kinderman Ramage

Kinderman, A. J. and Ramage, J. G. (1976) Computer generation of normal random variables. Journal of the American Statistical Association 71, 893-896.

Non "buggy" version

Example usage:

const libR = require('lib-r-math.js');
const {
    rng: {
        SuperDuper,
        normal: { KindermanRamage }
    },
    R: { numberPrecision }
} = libR;

//helper
const precision = numberPrecision(9);

const sd = new SuperDuper(0);
const km1 = new KindermanRamage(sd);

// At any time reset normal PRNG seed, with the reference to uniform PRNG
sd.init(1234);
km1.rng.init(1234); // same as above

// uses default uniform PRNG MersenneTwister with seed "0"
const km2 = new KindermanRamage();

km2.rng.init(0); // at any time reset PRNG with a new seed.

// bleed the normal PRNG
const k2 = km2.norm_rand(5);
precision(k2);
//[ 0.3216151, 1.23251561, 0.280369528, -1.17519641, -1.60471361 ]

// it's possible to bleed the uniform PRNG from the normal PRNG
km2.unif_rand();
//0.17655675252899528

km2.rng.unif_rand();
//0.6870228466577828

in R console

> RNGkind("Mersenne-Twister",normal.kind="Kinder")
> set.seed(0)
> rnorm(5)
[1]  0.3216151  1.2325156  0.2803695 -1.1751964
[5] -1.6047136
> runif(2)
[1] 0.1765568 0.6870228
>

Normal and Uniform distributions

Summary

In the Section 1 and 2 we discussed uniform and normal PRNG classes. These classes can be used by themselves but mostly intended to be consumed to generate random numbers with a particular distribution (like Uniform, Normal, Gamma, Weibull, Chi-square etc).

Uniform distribution

dunif, qunif, punif, runif

See R doc

These functions are created with the factory method Uniform taking as argument a uniform PRNG. Defaults to Mersenne-Twister.

Usage:

const libR = require('lib-r-math.js');
const {
    Uniform,
    rng: { SuperDuper }
} = libR;

//Create Uniform family of functions using "SuperDuper"
const uni1 = Uniform(new SuperDuper(0));

// Create Uniform family of functions using default "Mersenne-Twister"
const uni2 = Uniform();

// functions exactly named as in `R`
const { runif, dunif, punif, qunif } = uni2;

`dunif`

The density function. See R doc

typescript decl

declare function dunif(
  x: number | number[],
  min = 0,
  max = 1,
  asLog = false
): number | number[];

x: scalar or vector of quantiles
min, max lower and upper limits of the distribution. Must be finite.
asLog if true, results are given as ln.

Example:

const libR = require('lib-r-math.js');
const {
    Uniform,
    R: { numberPrecision, c }
} = libR;

//helper
const precision = numberPrecision(9);

const { runif, dunif, punif, qunif } = Uniform();

const x = [-1, 0, 0.4, 1, 2];

const d0 = dunif(0.5);
//1

const d1 = dunif(x);
// [ 0, 1, 1, 1, 0 ]  Everythin is 1 for inputs between 0 and 1

const d2 = dunif(x, 0, 2);
// [ 0, 0.5, 0.5, 0.5, 0.5 ]

const d3 = dunif(x, 0, 2, true);
precision(d3);
//[ -Infinity,  -0.693147181, -0.693147181, -0.693147181, -0.693147181 ]

`punif`

The probability function. See R doc

typescript decl

declare function punif(
  q: number | number[],
  min = 0,
  max = 1,
  lowerTail = true,
  logP = false
): number | number[];

x: scalar or vector of quantiles
min, max: lower and upper limits of the distribution. Must be finite.
lowerTail: if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP: if true, probabilities p are given as ln(p).

Example:

const libR = require('lib-r-math.js');
const {
    Uniform,
    R: { numberPrecision }
} = libR; // use default Mersenne-Twister PRNG

//helper
const precision = numberPrecision(9);

const { runif, dunif, punif, qunif } = Uniform();

const q = [-2, 0.25, 0.75, 2];

const p1 = punif(0.25);
// 0.25

const p2 = punif(q);
//[ 0, 0.25, 0.75, 1 ]

const p3 = punif(q, 0, 1, false);
precision(p3);
//[ 1, 0.75, 0.25, 0 ]

const p4 = punif(q, 0, 2, false, true);
precision(p4);
//[ 0, -0.133531393, -0.470003629, -Infinity ]

const p5 = punif(q, 0, 2, true, true);
precision(p5);
//[ 0, -0.133531393, -0.470003629, -Infinity ]

`qunif`

The quantile function. See R doc

typescript decl

declare function qunif(
  p: number | number[],
  min = 0,
  max = 1,
  lowerTail = true,
  logP = false
): number | number[];

p: scalar or vector of quantiles
min, max lower and upper limits of the distribution. Must be finite.
lowerTail if true (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP if true, probabilities p are given as ln(p).

Example:

const libR = require('lib-r-math.js');
const {
    Uniform,
    R: { numberPrecision, multiplex }
} = libR;

//helper
const precision = numberPrecision(9);

const { runif, dunif, punif, qunif } = Uniform();

const p = [0, 0.1, 0.5, 0.9, 1];

const q1 = qunif(0);
//0

const q2 = qunif(p, -1, 1, false);
//[ 1, 0.8, 0, -0.8, -1 ]

const log = multiplex(Math.log);

const q3 = qunif(log(p), -1, 1, false, true);
//[ 1, 0.8, 0, -0.8, -1 ]

`runif`

Generates random deviates. See R doc

typescript decl

declare function runif(
  n: number = 1,
  min: number = 0,
  max: number = 1
): number | number[];

n: number of deviates. Defaults to 1.
min, max lower and upper limits of the distribution. Must be finite.

Example:

const libR = require('lib-r-math.js');

const {
    Uniform,
    rng: { LecuyerCMRG },
    R: { numberPrecision }
} = libR;

//helper
const _9 = numberPrecision(9);

//explicit PRNG choice, seed = 1234
const lm = new LecuyerCMRG(1234);

const { runif, dunif, punif, qunif } = Uniform(lm);

const r1 = _9(runif(4));
//[ 0.472909817, 0.76978367, 0.216015397, 0.413843973 ]

const r2 = _9(runif(5, -1, 1, true));
//[ 0.122007235, 0.86544455, 0.0295475019, -0.184492403, 0.645749715 ]

Equivalent in R

RNGkind("L'Ecuyer-CMRG");
set.seed(1234);

#1
runif(4);
#[1] 0.4729098 0.7697837 0.2160154 0.4138440

runif(5,-1,1);
#[1]  0.1220072  0.8654446  0.0295475 -0.1844924  0.6457497

Normal distribution

dnorm, qnorm, pnorm, rnorm

Density, distribution function, quantile function and random generation for the normal distribution with mean equal to mean and standard deviation equal to sd. See R doc. These functions are created with the factory method Normal taking as optional argument a normal PRNG (defaults to Inversion.

Usage:

const libR = require('lib-r-math.js');
const {
    Normal,
    rng: {
        SuperDuper,
        normal: { AhrensDieter }
    }
} = libR;

//specify explicit PRNG's
const norm1 = Normal(new AhrensDieter(new SuperDuper(1234)));

//OR just go with defaults: "Inversion" and "Mersenne-Twister".
const norm2 = Normal(); //

//strip and use
const { rnorm, dnorm, pnorm, qnorm } = norm2;

`dnorm`

The density function of the Normal distribution. See R doc

typescript decl

declare function dnorm(
  x: number | number[],
  mu = 0,
  sigma = 1,
  asLog = false
): number | number[];

x:scalar or array of quantiles
mu: mean, default 0.
sigma: standard deviation, default 1.
asLog: give result as ln(..) value

Usage:

const libR = require('lib-r-math.js');
const {
    Normal,
    R: {
        numberPrecision,
        seq: _seq,
        c
    }
} = libR;

//helpers
const seq = _seq()();
const _9 = numberPrecision(9); //9 digits significance

const { rnorm, dnorm, pnorm, qnorm } = Normal();

const d1 = _9(dnorm(0));
//0.39894228

//x=3, µ=4, sd=2
const d2 = _9(dnorm(3, 4, 2));
//0.176032663

const d3 = _9(dnorm(-10));
//7.69459863e-23

//feed it also some *non-numeric*
const x = c(-Infinity, Infinity, NaN, seq(-4, 4));
const d4 = _9(dnorm(x));
/*[
  0,
  0,
  NaN,
  0.000133830226,
  0.00443184841,
  0.0539909665,
  0.241970725,
  0.39894228,
  0.241970725,
  0.0539909665,
  0.00443184841,
  0.000133830226 ]*/

const d5 = _9(dnorm(x, 0, 1, true));
/*[ -Infinity,
    -Infinity,
    NaN,
    -8.91893853,
    -5.41893853,
    -2.91893853,
    -1.41893853,
    -0.918938533,
    -1.41893853,
    -2.91893853,
    -5.41893853,
    -8.91893853 ]*/

Equivalent in R

dnorm(0);
#[1] 0.3989423

dnorm(3, 4, 2);
#[1] 0.1760327

dnorm(-10)
#[1] 7.694599e-23

x = c(-Inf, Inf, NaN, seq(-4, 4));
dnorm(x)
# [1] 0.0000000000 0.0000000000          NaN 0.0001338302 0.0044318484
# [6] 0.0539909665 0.2419707245 0.3989422804 0.2419707245 0.0539909665
#[11] 0.0044318484 0.0001338302

dnorm(x, 0,1, TRUE);
# [1]       -Inf       -Inf        NaN -8.9189385 -5.4189385 -2.9189385
# [7] -1.4189385 -0.9189385 -1.4189385 -2.9189385 -5.4189385 -8.9189385

`pnorm`

The distribution function of the Normal distribution. See R doc

typescript decl

declare function pnorm(
  q: number | number[],
  mu = 0,
  sigma = 1,
  lowerTail = true,
  logP = false
): number | number[];

q:scalar or array of quantiles
mu: mean (default 0)
sigma: standard deviation (default 1)
lowerTail: if true (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP: give result as log value

Usage:

const libR = require('lib-r-math.js');
const {
    Normal,
    R: { numberPrecision, multiplex, seq: _seq }
} = libR;

const { rnorm, dnorm, pnorm, qnorm } = Normal();

// some helpers
const seq = _seq()();
const _9 = numberPrecision(9); //9 digit significance

//data
const q = seq(-1, 1);

const p1 = _9(pnorm(q));
//[ 0.158655254, 0.5, 0.841344746 ]

const p2 = _9(pnorm(q, 0, 1, false));
//[ 0.841344746, 0.5, 0.158655254 ]

const p3 = _9(pnorm(q, 0, 1, false, true));
//[ -0.172753779, -0.693147181, -1.84102165 ]

Equivalent in R

pnorm(-1:1);
#[1] 0.1586553 0.5000000 0.8413447

pnorm(-1:1, lower.tail=FALSE);
#[1] 0.8413447 0.5000000 0.1586553

pnorm(-1:1, log.p= TRUE);
#[1] -0.1727538 -0.6931472 -1.8410216

`qnorm`

The quantile function of the Normal distribution. See R doc

typescript decl

declare function qnorm(
  p: number | number[],
  mu = 0,
  sigma = 1,
  lowerTail = true,
  logP = false
): number | number[];

p: probabilities (scalar or array).
mu: normal mean (default 0).
sigma: standard deviation (default 1).
logP: probabilities are given as ln(p).

Usage:

const libR = require('lib-r-math.js');
const {
    Normal,
    R: { multiplex, seq: _seq, numberPrecision }
} = libR;

//some helpers
const log = multiplex(Math.log);
const _9 = numberPrecision(9);
const seq = _seq()();

const { rnorm, dnorm, pnorm, qnorm } = Normal();

//some data
const p = seq(0, 1, 0.25);
//[0, 0.25, 0.5, 0.75, 1]

const q1 = _9(qnorm(0));
//-Infinity

const q2 = _9(qnorm(p, 0, 2));
//[ -Infinity, -1.3489795, 0, 1.3489795, Infinity ]

const q3 = _9(qnorm(p, 0, 2, false));
//[ Infinity, 1.3489795, 0, -1.3489795, -Infinity ]

//same as q3
const q4 = _9(qnorm(log(p), 0, 2, false, true));
//[ Infinity, 1.3489795, 0, -1.3489795, -Infinity ]

Equivalent in R

p = seq(0, 1, 0.25);

qnorm(0);
#[1] -Inf

qnorm(p, 0, 2);
#[1]     -Inf -1.34898  0.00000  1.34898      Inf

qnorm(p, 0, 2, FALSE);
#[1]      Inf  1.34898  0.00000 -1.34898     -Inf

qnorm(log(p), 0, 2, FALSE, TRUE);
#[1]      Inf  1.34898  0.00000 -1.34898     -Inf

`rnorm`

Generates random normal deviates. See R doc

typescript decl

declare function rnorm(n = 1, mu = 0, sigma = 1): number | number[];

n: number of deviates
mu: mean of the distribution. Defaults to 0.
sigma: standard deviation. Defaults to 1.

Usage:

const libR = require('lib-r-math.js');
const {
    Normal,
    R: { numberPrecision }
} = libR;

//helper
const _9 = numberPrecision(9); // 9 digits

//default Mersenne-Twister/Inversion
const { rnorm, dnorm, pnorm, qnorm } = Normal();

const r1 = _9(rnorm(5));
//[ 1.26295428, -0.326233361, 1.32979926, 1.27242932, 0.414641434 ]

const r2 = _9(rnorm(5, 2, 3));
//[ -2.61985013, -0.785701104, 1.11583866, 1.98269848, 9.21396017 ]

Equivalent in R

RNGkind("Mersenne-Twister",normal.kind="Inversion")
set.seed(0)

rnorm(5)
#[1]  1.2629543 -0.3262334
#[3]  1.3297993  1.2724293
#[5]  0.4146414

rnorm(5,2,3)
#[1] -2.6198501 -0.7857011
#[3]  1.1158387  1.9826985
#[5]  9.2139602

Other Probability Distributions

summary

libRmath.so contains 19 probability distributions (other then Normal and Uniform) with their specific density, quantile and random generators, all are ported and have been verified to yield the same output.

Beta distribution

dbeta, qbeta, pbeta, rbeta

See R doc. See wiki.

These functions are members of an object created by the Beta factory method. The factory method needs an instance of an optional normal PRNG. Various instantiation methods are given below.

Usage:

const libR = require('lib-r-math.js');
const { Beta, rng: { SuperDuper, normal: { BoxMuller } } } = libR;

// explicit use of PRNG's
const explicitB = Beta(new BoxMuller(new SuperDuper(0))); //

// go with defaults 'MersenneTwister" and "Inversion"
const defaultB = Beta();

// Or just go with Default.. defaults to PRNG "Inversion" and "Mersenne-Twister"
const { dbeta, pbeta, qbeta, rbeta } = defaultB;

`dbeta`

The density function of the Beta distribution. See R doc.

$ \frac{\Gamma(a+b)}{Γ(a) Γ(b)} x^{(a-1)}(1-x)^{(b-1)} $

typescript decl

declare function dbeta(
  x: number | number[],
  shape1: number,
  shape2: number,
  ncp = undefined,
  asLog = false
): number | number[];

x: scalar or array of quantiles. 0 <= x <= 1
shape1: non-negative a parameter of the Beta distribution.
shape2: non-negative b parameter of the Beta distribution.
ncp: non centrality parameter. Note: undefined is different then 0
asLog: return result as ln(p)

const libR = require('lib-r-math.js');
const { Beta, R: { numberPrecision } } = libR;

//helpers, 9 digits precision
const _9 = numberPrecision(9);

//just go with Default.. uses Normal(), defaults to PRNG "Inversion" and "Mersenne-Twister"
const { dbeta, pbeta, qbeta, rbeta } = Beta();

//1. ncp argument = 1
const d1 = _9(dbeta(0.4, 2, 2, 1));
//1.28724574

//2., No named arguments in JS, so use undefined to skip
const d2 = _9(dbeta(0.4, 2, 2, undefined, true));
//0.364643114

//3
const d3 = _9(dbeta(0.4, 2, 2, 1, true));
//0.252504851

//4
const d4 = _9(
    dbeta(
        [0, 0.2, 0.4, 0.8, 1, 1.2],
        2,
        2)
);
//[ 0, 0.96, 1.44, 0.96, 0, 0 ]

Equivalent in R

#1
dbeta(0.4,2,2, ncp=1)
#[1] 1.287246

#2
dbeta(0.4,2,2, log = TRUE)
#[1] 0.3646431

#3
dbeta(0.4,2,2, ncp=1, TRUE)
#[1] 0.2525049

#4
dbeta( c(0, 0.2, 0.4, 0.8, 1, 1.2), 2, 2)
#[1] 0.00 0.96 1.44 0.96 0.00 0.00

`pbeta`

The cumulative probability function of the Beta distribution. See R doc.

declare function pbeta(
  q: number | number[],
  shape1: number,
  shape2: number,
  ncp?: number,
  lowerTail: boolean = true,
  logP: boolean = false
): number | number[];

p: quantiles. 0 <= x <= 1
shape1: non-negative a parameter of the Beta distribution.
shape2: non-negative b parameter of the Beta distribution.
ncp: non centrality parameter. Note: undefined is different then 0
lowerTail: if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP: return probabilities as ln(p)

Usage:

const libR = require('lib-r-math.js');
const {
    Beta,
    R: { multiplex, numberPrecision, seq: _seq }
} = libR;

//helpers
// 9 digit precision
const _9 = numberPrecision(9);
const log = multiplex(Math.log);
const seq = _seq()();

//just go with Default.. uses Normal(), defaults to PRNG "Inversion" and "Mersenne-Twister"
const { dbeta, pbeta, qbeta, rbeta } = Beta();
const q = seq(0, 1, 0.2);

//1.
const p1 = _9(pbeta(0.5, 2, 5));
//0.890625

//2.
const p2 = _9(pbeta(0.5, 2, 5, 4));
//0.63923843

//3.
const p3 = _9(pbeta(q, 2, 5, 4));
//[ 0, 0.106517718, 0.438150345, 0.813539396, 0.986024517, 1 ]

//4.
const p4 = _9(pbeta(q, 2, 5, undefined));
//[ 0, 0.345027474, 0.76672, 0.95904, 0.9984, 1 ]

//5. result as as ln(p)
const p5 = _9(pbeta(q, 2, 5, undefined, true, true));
/*[
  -Infinity,     -1.06413123,    -0.265633603,
  -0.0418224949, -0.00160128137,  0
  ]*/

Equivalent in R

q = c(0, 0.2, 0.4, 0.6, 0.8, 1);

#1
pbeta(0.5, 2, 5);
#[1] 0.890625

#2
pbeta(0.5, 2, 5, 4)
#[1] 0.6392384

#3
pbeta(q, 2, 5, 4);
#[1] 0.0000000 0.1061302 0.4381503 0.8135394
#[5] 0.9860245 1.0000000

#4
pbeta(q, 2, 5);
#[1] 0.00000 0.34464 0.76672 0.95904 0.99840 1.00000

#5
pbeta(q, 2, 5, log.p=TRUE)
#[1]         -Inf -1.065254885 -0.265633603
#[4] -0.041822495 -0.001601281  0.000000000

`qbeta`

The quantile function of the Beta distribution. See R doc.

typescript decl

declare function qbeta(
  p: number | number[],
  shape1: number,
  shape2: number,
  ncp = undefined,
  lowerTail = true,
  logP = false
): number | number[];

p: quantiles (scalar or array).
shape1: non-negative a parameter of the Beta distribution.
shape2: non-negative b parameter of the Beta distribution.
ncp: non centrality parameter. Note: undefined is different then 0
lowerTail: if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP: return probabilities as ln(p).

Usage:

const libR = require('lib-r-math.js');
const { Beta, R: { multiplex, numberPrecision } } = libR;

//helpers
const ln = multiplex(Math.log);
const _9 = numberPrecision(9); // 9 digits precision

const { dbeta, pbeta, qbeta, rbeta } = Beta();

//take probabilities in steps of 25%
const p = [0, 0.25, 0.5, 0.75, 1];

//1. always zero, regardless of shape params, because 0 ≤ x ≤ 1.
const q1 = _9(qbeta(0, 99, 66));
//0

//2. 
const q2 = _9(qbeta(p, 4, 5));
//[ 0, 0.329083427, 0.440155205, 0.555486315, 1 ]

//3 ncp = 3
const q3 = _9(qbeta(p, 4, 5, 3));
//[ 0, 0.406861514, 0.521344641, 0.631881288, 1 ]

//4. ncp = undefined, lowerTail = false, logP=false(default)
const q4 = _9(qbeta(p, 4, 5, undefined, false)); //
//[ 1, 0.555486315, 0.440155205, 0.329083427, 0 ]

//5. same as [5] but, logP=true,
const q5 = _9(qbeta(
    ln(p),
    4,
    5,
    undefined,
    false,
    true //p as ln(p)
));
//[ 1, 0.555486315, 0.440155205, 0.329083427, 0 ]

Equivalent in R

p = c(0,.25,.5,.75,1);
#1
qbeta(0,99,66);
#[1] 0

#2
qbeta(p, 4,5);
#[1] 0.0000000 0.3290834 0.4401552 0.5554863 1.0000000

#3
qbeta(p, 4,5,3);
#[1] 0.0000000 0.4068615 0.5213446 0.6318813 1.0000000

#4
qbeta(p, 4,5, lower.tail = FALSE);
#[1] 1.0000000 0.5554863 0.4401552 0.3290834
#[5] 0.0000000

#5
qbeta(  log(p)  ,4,5, lower.tail = FALSE, log.p=TRUE);
#[1] 1.0000000 0.5554863 0.4401552 0.3290834 0.0000000

`rbeta`

Generates random deviates for the Beta distribution. See R doc.

typescript decl

declare function rbeta(
  n: number,
  shape1: number,
  shape2: number,
  ncp = 0 // NOTE: normally the default is "undefined", but not here
): number | number[];

n: number of deviates
shape1: non-negative a parameter of the Beta distribution.
shape2: non-negative b parameter of the Beta distribution.
ncp: non centrality parameter.

const libR = require('lib-r-math.js');
const {
    Beta,
    rng: {
        LecuyerCMRG,
        normal: { Inversion }
    },
    R: { multiplex, numberPrecision }
} = libR;

//helpers
const ln = multiplex(Math.log); //
const _9 = numberPrecision(9);

const lc = new LecuyerCMRG(0);
const { dbeta, pbeta, qbeta, rbeta } = Beta(new Inversion(lc));

//1.
const r1 = _9(rbeta(5, 0.5, 0.5));
//[ 0.800583949,  0.962961579, 0.700710737,  0.169742664, 0.0169845581 ]

//2.
const r2 = _9(rbeta(5, 2, 2, 4));
//[ 0.940977213, 0.803938008, 0.762066155, 0.775315234, 0.0395894783 ]

//3. // re-initialize seed
lc.init(0);

//3
const r3 = _9(rbeta(5, 2, 2));
//[ 0.37955891, 0.240142694, 0.425371111, 0.935280271, 0.636741506 ]

//4.
const r4 = _9(rbeta(5, 2, 2, 5));
//[ 0.532034853, 0.985042931, 0.724819159, 0.67645358, 0.837372377 ]

Same values as in R

Equivalent in R

RNGkind("L'Ecuyer-CMRG", normal.kind ="Inversion")
set.seed(0)

#1
rbeta(5, 0.5, 0.5)
#[1] 0.80058395 0.96296158 0.70071074 0.16974266 0.01698456

#2
rbeta(5, 2, 2, 4)
#[1] 0.94097721 0.80393801 0.76206615 0.77531523 0.03958948

set.seed(0)
#3
rbeta(5, 2, 2);
#[1] 0.3795589 0.2401427 0.4253711 0.9352803 0.6367415

#4
rbeta(5, 2, 2, 5);
#[1] 0.5320349 0.9850429 0.7248192 0.6764536 0.8373724

Binomial distribution

dbinom, qbinom, pbinom, rbinom

See R doc and wiki.

These functions are members of an object created by the Binomial factory method. The factory method needs an instance of a normal PRNG. Various instantiation methods are given below.

Usage:

const libR = require('lib-r-math.js');
const { Binomial, rng: { LecuyerCMRG } } = libR;

// explicit use if PRNG
const lc = new LecuyerCMRG(0);
const explicitB = Binomial(lc);

//default, used "Inversion" and "MersenneTwister"
const defaultB = Binomial();

const { dbinom, pbinom, qbinom, rbinom } = defaultB;

`dbinom`

The density function of the Binomial distribution. See R doc

$p(x) = \frac{n!}{x!(n-x)!} p^{x} (1-p)^{n-x}$

typescript decl

declare function dbinom(
  x: number,
  size: number,
  p: number,
  asLog = false
): number | number[];

x: scalar or array of quantiles.
size: number of trails
p: probability of success.
asLog: return result as ln(p)

Usage:

const libR = require('lib-r-math.js');
const {
    Binomial,
    R: { numberPrecision, seq: _seq }
} = libR;

//helper, 9 digits precision
const _9 = numberPrecision(9);
const seq = _seq()();

//some data
const x = seq(1, 4);

//Binomial()  uses Normal() as default argument,
const { dbinom, pbinom, qbinom, rbinom } = Binomial();

//1. 2 successes out of 4 trials, with success probility 0.3
const d1 = _9(dbinom(2, 4, 0.3));
//0.2646

//2. same as [1], but results as log
const d2 = _9(dbinom(2, 4, 0.3, true));
//-1.32953603

//3. all possibilities out of 4 trials
const d3 = _9(dbinom(x, 4, 0.3));
//[ 0.4116, 0.2646, 0.0756, 0.0081 ]

//4
const d4 = _9(dbinom(x, 4, 0.3, true));
//[ -0.887703275, -1.32953603, -2.582299, -4.81589122 ]

Equivalent in R

#1
dbinom(2,4,0.3)
#[1] 0.2646

#2
dbinom(2,4,0.3, TRUE)
#[1] -1.329536

#3
dbinom(c(1,2,3,4),4,0.3)
#[1] 0.4116 0.2646 0.0756 0.0081

#4
dbinom(c(1,2,3,4),4,0.3, TRUE)
#[1] -0.8877033 -1.3295360 -2.5822990  -4.8158912

`pbinom`

The cumulative probability function of the Binomial distribution. See R doc

declare function pbinom(
  q: number | number[],
  size: number,
  prob: number,
  lowerTail = true,
  logP = false
): number | number[];

q: scalar or array of quantiles.
size: number of trails
prob: probability of success.
lowerTail: if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP: return result as ln(p)

Usage:

const libR = require('lib-r-math.js');
const {
    Binomial,
    R: { numberPrecision, seq: _seq }
} = libR;

//helper, 9 digits precision
const _9 = numberPrecision(9);
const seq = _seq()();

const { dbinom, pbinom, qbinom, rbinom } = Binomial();
//some data

const q = seq(0, 4);
//1.
const p1 = pbinom(4, 4, 0.5);
//1

//2.
const p2 = _9(pbinom(q, 4, 0.5));
//[ 0.0625, 0.3125, 0.6875, 0.9375, 1 ]

//4.
const p3 = _9(pbinom(q, 4, 0.5, false, true));
/*[ -0.0645385211,  -0.374693449,  -1.16315081,  -2.77258872,  -Infinity
]*/

Equivalent in R

q = c(0, 1, 2, 3, 4);
#1
pbinom(4, 4, 0.5)
#[1] 1

#2
pbinom(q, 4, 0.5)
#[1] 0.0625 0.3125 0.6875 0.9375 1.0000

#3
pbinom(q, 4, 0.5, FALSE, TRUE)
#[1] -0.06453852 -0.37469345 -1.16315081
#[4] -2.77258872        -Inf

`qbinom`

The quantile function of the Binomial distribution. See R doc

typescript decl

declare function qbinom(
  p: number | number[],
  size: number,
  prob: number,
  lowerTail = true,
  logP = false
): number | number[];

p: scalar or array of quantiles.
size: number of trails
prob: probability of success.
lowerTail: if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
LogP: return result as ln(p)

Usage:

const libR = require('lib-r-math.js');
const {
    Binomial,
    R: { multiplex, numberPrecision, seq: _seq }
} = libR;

//helpers
const _9 = numberPrecision(9);
const log = multiplex(Math.log);
const seq = _seq()();

const { dbinom, pbinom, qbinom, rbinom } = Binomial();

//data
const p = seq(0, 1, 0.25); //[0, 0.25, 0.5, 0.75, 1];

//1
const q1 = _9(qbinom(0.25, 4, 0.3));
//1

//2.
const q2 = _9(qbinom(p, 40, 0.3));
//[0 10 12 14 40]

//3.
const q3 = _9(qbinom(p, 40, 0.3, false));
//[ 40, 14, 12, 10, 0 ]

//4.  same as 3.
const q4 = _9(qbinom(log(p), 40, 0.3, false, true));
//[ 40, 14, 12, 10, 0 ]

Equivalent in R

p = seq(0,1,0.25); #c(0, 0.25, 0.5, 0.75, 1);

#1
qbinom(.25,4,.3)
#[1] 1

#2
qbinom(p, 40,.3)
#[1]  0 10 12 14 40

#3
qbinom(p, 40,.3, FALSE)
#[1] 40 14 12 10  0

#4
qbinom(log(p), 40,.3, FALSE, TRUE)
#[1] 40 14 12 10  0

`rbinom`

Generates random beta deviates for the Binomial distribution. See R doc.

typescript decl

declare function rbinom(
  n: number,
  size: number,
  prop: number
): number | number[];

n: number of deviates
size: number of trails
prob: probability of success.

Usage:

const libR = require('lib-r-math.js');
const {
    Binomial,
    rng: { KnuthTAOCP2002 }
} = libR;

const kn = new KnuthTAOCP2002(1234);
const { dbinom, pbinom, qbinom, rbinom } = Binomial(kn);

//1.
const r1 = rbinom(2, 40, 0.5);
//[ 24, 19 ]

//2.
const r2 = rbinom(3, 20, 0.5);
//[ 11, 13, 13 ]

//3.
const r3 = rbinom(2, 10, 0.25);
//[ 2, 2 ]

Equivalent in R

RNGkind("Knuth-TAOCP-2002")
set.seed(1234)

#1
rbinom(2, 40, 0.5);
#[1] 24 18

#2
rbinom(3, 20, 0.5);
#[1] 11 13 13

#3
rbinom(2, 10, 0.25);
#[1] 2 2

Negative Binomial distribution

dnbinom, pnbinom, qnbinom, rnbinom.

See [R doc](https: //stat.ethz.ch/R-manual/R-devel/library/stats/html/NegBinomial.html) See wiki

These functions are members of an object created by the NegativeBinomial factory method. This factory method needs an instance of a normal PRNG. Various instantiation methods are given below.

Usage:

const libR = require('lib-r-math.js');
const {
    NegativeBinomial,
    rng: {
        SuperDuper,
        normal: { BoxMuller }
    }
} = libR;

//explicit use PRNG's
const bm = new BoxMuller(new SuperDuper(0));
const explicitNB = NegativeBinomial(bm);

//default uses PRNG "Inverion" and "MersenneTwister"
const defaultNB = NegativeBinomial();

const { dnbinom, pnbinom, qnbinom, rnbinom } = defaultNB;

`dnbinom`

The density function of the Negative Binomial distribution.

$ \frac{Γ(x+n)}{Γ(n) x!} p^{n} (1-p)^{x} $

See [R doc] (https: //stat.ethz.ch/R-manual/R-devel/library/stats/html/NegBinomial.html).

typescript decl

declare function dnbinom(
  x: number | number[],
  size: number,
  prob?: number,
  mu?: number,
  asLog = false
): number | number[];

x: non-negative integer quantiles. Number of failures before reaching size successes.
size: target for number of successful trials, or dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer.
prob: probability of success in each trial. 0 < prob <= 1
mu: alternative parametrization via mean: see ‘Details’ section.
asLog: if true, probabilities p are given as ln(p).

Usage:

const libR = require('lib-r-math.js');
const {
    NegativeBinomial,
    R: { seq: _seq, numberPrecision }
} = libR;

//some helpers
const seq = _seq()();
const _9 = numberPrecision(9);

const { dnbinom, pnbinom, qnbinom, rnbinom } = NegativeBinomial();

//some data
const x = seq(0, 10, 2);

//1.
const d1 = _9(dnbinom(x, 3, 0.5));
//[ 0.125, 0.1875, 0.1171875, 0.0546875, 0.0219726562, 0.00805664062 ]

//2. alternative presentation with `mu` = n*(1-p)/p
const d2 = _9(dnbinom(x, 3, undefined, 3 * (1 - 0.5) / 0.5));
//[ 0.125, 0.1875, 0.1171875, 0.0546875, 0.0219726562, 0.00805664062 ]

//3
const d3 = _9(dnbinom(x, 3, undefined, 3 * (1 - 0.5) / 0.5, true));
/*[ -2.07944154, -1.67397643,  -2.14398006, -2.90612011,
    -3.8179565,  -4.82125861
]*/

Equivalent in R

#1
dnbinom(0:10, size = 3, prob = 0.5);
# [1] 0.125000000 0.187500000 0.187500000 0.156250000 0.117187500 0.082031250
# [7] 0.054687500 0.035156250 0.021972656 0.013427734 0.008056641

#2
dnbinom(0:10, size = 3, mu = 3*(1-0.5)/0.5);
# [1] 0.125000000 0.187500000 0.187500000 0.156250000 0.117187500 0.082031250
# [7] 0.054687500 0.035156250 0.021972656 0.013427734 0.008056641

dnbinom(0:10, size = 3, mu = 3*(1-0.5)/0.5, log=T);
# [1] -2.079442 -1.673976 -1.673976 -1.856298 -2.143980 -2.500655 -2.906120
# [8] -3.347953 -3.817956 -4.310433 -4.821259

`pnbinom`

The gives the cumulative probability function of the Negative Binomial distribution. See [R doc](https: //stat.ethz.ch/R-manual/R-devel/library/stats/html/NegBinomial.html).

typescript decl

 declare function pnbinom(
    q: number | number[],
    size: number,
    prob?: number,
    mu?: number,
    lowerTail = true
    logP =  false
  ): number|number[]

q: non-negative integer quantiles.
size: target for number of successful trials, or dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer.
prob: probability of success in each trial. 0 < prob <= 1
mu: alternative parametrization via mean: see ‘Details’ section.
lowerTail: if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
logP: if true, probabilities p are given as ln(p).

Usage:

const libR = require('lib-r-math.js');
const {
    NegativeBinomial,
    R: { numberPrecision, seq: _seq, c }
} = libR;

//some helpers
const seq = _seq()();
const _9 = numberPrecision(9);

const { dnbinom, pnbinom, qnbinom, rnbinom } = NegativeBinomial();

//some data
const x = c(seq(0, 6), Infinity);
//[ 0, 1, 2, 3, 4, 5, 6, Infinity ]

//1.
const p1 = _9(pnbinom(x, 3, 0.5));
//[ 0.125, 0.3125, 0.5, 0.65625, 0.7734375, 0.85546875, 0.91015625, 1 ]

//2. alternative presentation of 1 with mu = n(1-p)/p
const p2 =

Usage no npm install needed!

README

libRmath.js

Summary

Node and Web

Installation

serverside

web

Table of Contents

Differences with R

Helper functions for porting R programs

Summary

div

mult

asArray

sum

summary

numberPrecision

any

arrayrify (DEPRICATED use multiplex)

R example

Javascript equivalent

each

flatten or c (alias)

map

selector

seq

multiplex

timeseed

Uniform Pseudo Random Number Generators

Summary

The 7 Uniform Random Number Generators

Improvements compared to R

General Usage

"Mersenne Twister"

"Wichmann-Hill"

"Marsaglia-Multicarry"

"Super Duper"

"Knuth TAOCP"

"Knuth TAOCP 2002"

"L'Ecuyer-CMRG":

Normal distributed Random Number Generators

Summary

General Use

"Ahrens Dieter"

Box Muller

Buggy Kinderman Ramage

Inversion

Kinderman Ramage

Normal and Uniform distributions

Summary

Uniform distribution

dunif

punif

qunif

runif

Normal distribution

dnorm

pnorm

qnorm

rnorm

Other Probability Distributions

summary

Beta distribution

dbeta

pbeta

qbeta

rbeta

Binomial distribution

dbinom

pbinom

qbinom

rbinom

Negative Binomial distribution

dnbinom

pnbinom

Helper functions for porting `R` programs

`div`

`mult`

`asArray`

`sum`

`summary`

`numberPrecision`

`any`

`arrayrify` (DEPRICATED use `multiplex`)

`each`

`flatten` or `c` (alias)

`map`

`selector`

`seq`

`multiplex`

`timeseed`

`dunif`

`punif`

`qunif`

`runif`

`dnorm`

`pnorm`

`qnorm`

`rnorm`

`dbeta`

`pbeta`

`qbeta`

`rbeta`

`dbinom`

`pbinom`

`qbinom`

`rbinom`

`dnbinom`

`pnbinom`