# Addons/stats/distribs

**stats/distribs** - Working with distributions

- Includes verbs for working with the Normal and Uniform distributions
- Other distributions are planned but contributions are very welcome!

Browse history, source and examples using Trac.

## Contents

## Verbs available

### Normal Distribution

**Available in z locale**

`dnorm`v Normal probability density function `pnorm`v Normal cumulative distribution function `pnorm_f`v Normal cumulative distribution function `pnorm_ut`v Upper Tail version of pnorm `qnorm`v Quantile function for Normal distribution `qnorm_ut`v Upper Tail version of qnorm `rnorm`v Random deviates from Normal distribution `tomusigma`v Converts from N[0,1] to N[mu,sigma] `tostd`v Converts from N[mu,sigma] to N[0,1]

**Available in pdistribs locale **

`erf`v Error function `erfc`v Complementary Error function `erfinv`v Inverse of Error function

### Uniform Distribution

**Available in z locale**

`dunif`v Uniform probability density function `punif`v Uniform cumulative distribution function `qunif`v Quantile function for Uniform distribution `runif`v Random deviates from Uniform distribution

## Installation

Use JAL/Package Manager to install the `stats/distribs` addon.

## Usage

Load the `stats/distribs` addon with the following line

load 'stats/distribs'

If you wish to only work with one type of distribution you can load the individual scripts as follows:

load 'stats/distribs/normal'

To see the sampler of usage, open and inspect the test script for the distribution of interest. For example:

- test_normal.ijs script.
- test_uniform.ijs script.

### pnorm01

There are two algorithms here. `pnormh` is more accurate but slower than `pnorm01_f`.

`pnormh` uses built-in primitives and is due originally to Ewart Shaw with some modifications by Roger Hui. It is from from Abramovitz and Stegum 26.2.29 (solved for P) and is also defined in the stats/base/distribution.ijs script.

`pnorm01_f` is coded from a Chebychev expansion in Abramovitz and Stegum 26.2.17. It achieves a maximum absolute error of less than 7.46e_8 over the argument range (_5,5) and less than 0.2 percent relative error.

The `pnorm01` used in the normal script uses the `pnormh` algorithm for `y` values between `_7` and `7` and the `pnorm01_f` algorithm outside that range. This is because the `pnormh` algorithm appears to be unstable outside that range giving spurious results: eg:

0j17 ": ,.pnormh_pnormal_ 10 20 30 40 1.00000000000000020 1.00000000000000040 1.00000000000000090 0.50000000000000000

Note also:

pnormh_pnormal_ _ |NaN error: pnormh_pnormal_ | pnormh_pnormal_ _ pnormh_pnormal_ __ |NaN error: pnormh_pnormal_ | pnormh_pnormal_ __

Whereas:

0j17 ": ,.pnorm01_pnormal_ 10 20 30 40 _ __ 1.00000000000000000 1.00000000000000000 1.00000000000000000 1.00000000000000000 1.00000000000000000 0.00000000000000000

### qnorm01

The explicit `qnorm01` was based on the tacit code on EwartShaw's page EwartShaw/N01CdfInv. An explicit form was developed to improve the performance and ensure the desired results for 0 and 1 i.e.

qnorm01 0 1 __ _

Based on the suggestion of John Randall in this forum thread Fraser Jackson and Ric Sherlock coded the following J version of the algorithm described by P J Acklam. However the algorithm included in the Addon, while slightly slower than `qnorm01_acklam_fast` below, uses the same space and is considerably faster and leaner than `qnorm01_acklam_accurate`. At the same time it is more accurate than either.

qnorm01_acklam_fast=: 3 : 0 z=. ,y s=. ($z)$0 assert. (0<:z) *. z<:1 NB. y outside meaningful bounds a=. _3.969683028665376e01 2.209460984245205e02 _2.759285104469687e02 a=. |.a, 1.383577518672690e02 _3.066479806614716e01 2.506628277459239e00 b=. _5.447609879822406e01 1.615858368580409e02 _1.556989798598866e02 b=. |.b, 6.680131188771972e01 _1.328068155288572e01 1 c=. _7.784894002430293e_03 _3.223964580411365e_01 _2.400758277161838e00 c=. |.c, _2.549732539343734e00 4.374664141464968e00 2.938163982698783e00 d=. 7.784695709041462e_03 3.224671290700398e_01 2.445134137142996e00 d=. |.d, 3.754408661907416e00 1 NB. Define break-points. p_low=. 0.02425 p_high=. 1 - p_low NB. Rational approximation for lower region. v=. (0 < z) *. z < p_low q=. %: _2*^. v#z s=. ((c p. q) % d p. q) (I.v)} s NB. Rational approximation for central region. v=. (p_low <: z) *. z <: p_high q=. (v#z) - 0.5 r=. *:q s=. (q * (a p. r) % b p. r) (I.v)} s NB. Rational approximation for upper region. v=. (p_high < z) *. z < 1 q=. %: _2* ^. 1- v#z s=. (-(c p. q) % d p. q) (I.v)} s NB. equal to 0 or 1 s=. __ (I. z=0)} s s=. _ (I. z=1)} s ($y)$s ) qnorm01_acklam_accurate=: 3 : 0 z=. ,y s=. qnorm01_acklam_fast z NB. Refinement using Halley's rational method v=. (0 < z)*. z < 1 q=. v#s NB. x e=. (v#z) -~ -: erfc -q% %:2 NB. error u=. e * (%:2p1) * ^ (*:q) % 2 NB. f(z)/df(z) NB. q=. q - u NB. Newton's method s=. (q - u% >:q*u%2) (I.v)}s NB. Halley's method ($y)$s )

### rnorm01

Brian Schott, Ric Sherlock and others contributed code to the forum discussion of this function. This explicit version below follows closely the Box-Muller form of Schott and shows the structure of steps in the code more clearly for those not experienced with the tacit form. However the tacit version used in the Addon is leaner.

rnorm01=: 3 : 0 n=. >. -: */y a=. %: _2* ^. runif01 n b=. 2* o. runif01 n r1=. a * 2 o. b r2=. a * 1 o. b y$r1,r2 )

## See Also

- Essays/Normal CDF
- Essays/Extended Precision Functions
- stats/base/distribution.ijs script
- Ewart Shaw's Vector article Hypergeometric Functions and CDFs in J
- Forum thread on generating standard normal variates
- Forum thread on Inverse Normal CDF
- P J Acklam's algorithm for Inverse Normal CDF
- EwartShaw/N01CdfInv

## Contributed by

. Ric Sherlock and Fraser Jackson