ppnd16

An implementation of the PPND16 function/AS 241 algorithm defined by Michael J. Wichura in Algorithm AS 241: The Percentage Points of the Normal Distribution. The implementation is faithful enough to retain the original's SHOUTINESS and IFAULT error code, so you will probably want to use the normalCdfInv(p, mean, sd) wrapper function instead.

Correctness has been verified by comparing to R's qnorm function, which uses the same algorithm.

Full API

function normalCdfInv(p: number, mean = 0, sd = 1): number
function PPND16(P: number, IFAULT: { value: 0 | 1 }): number

Usage

import { normalCdfInv } from "ppnd16";

// Log a random z-score to the console
// The Number.EPSILON fallback ensures p is greater than 0
console.log(normalCdfInv(Math.random() || Number.EPSILON));

This package can be used in Deno and modern browsers via Skypack.

A Note on Accuracy

Like qnorm, this package is affected by floating-point rounding errors and related issues. This means, among other things, that two of the paper's test cases fail because of a mismatch in the last digit of the result:

// Expected: -0.6744897501960817
// Actual:   -0.6744897501960817
// Result:   PASS
normalCdfInv(0.25);

// Expected: -3.090232306167814
// Actual:   -3.090232306167813
// Result:   FAIL
normalCdfInv(0.001);

// Expected: -9.262340089798408
// Actual:   -9.262340089798405
// Result:   FAIL
normalCdfInv(1e-20);

It's also worth noting that passing p values like 0.000001 and 0.999999 produces asymmetrical results. You can work around this using the following:

import { normalCdfInv as base } from "ppnd16";

// normalCdfInv(0.000001)             = -4.753424308822899
// normalCdfInv(0.999999)             =  4.753424308817089
// normalCdfInv(0.000001, 0, 1, true) =  4.753424308822899
export function normalCdfInv(p, mean = 0, sd = 1, flip = false) {
  return mean + (base(p) * (flip ? -1 : 1) * sd);
}