Calculating Probabilities
Vector 12:2, page 109
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 109
Alan Sykes rightly makes the point in Calculating Probabilities for Elementary Distributions (Vector Vol.12 No.1) that anyone who is in the business of doing serious probability calculations should regard a function LNNCR which computes the logarithm of nCr as essential. The function given by Alan can be appreciably shortened by using "each". First a basic function for scalar n and r is:
which can be generalised to Alan’s LNNCR, where the right argument is an array of any shape, by:
If only a small number of values are required, which is often likely to be the case in practice, this is the simplest way to go about things.
However, as Alan points out, if the right argument is an array with a largish number of items, an efficiency gain is achieved by computing many nCrs, and then selecting the ones which are actually wanted. An extension of lnncr to do this is
Alan says that extending this to array left arguments as for the "shriek" primitive is tricky and calls for a volunteer. I suggest:
Unless Alan has some subtlety in mind which is escaping me, it seems I have just volunteered! Turning to HYPERGEOM, the function which evaluates hypergeometric probabilities, I prefer the form:
This is not only shorter than Alan’s version, but also emphasises the essential symmetry of the calculation, which is not apparent in the mathematical notation, namely that the n,N-0,n which makes up the left argument to lnncr refer to the overall sizes of population and sample, while r,m-0,r refer to the respective numbers of "successes".
Continuing Alan’s illustration of the National Lottery odds, the inverse probabilities of gaining 3, 4, 5, 6 successes, that is the odds expressed as "1 in …" are given by
Next, Poisson probabilities and cumulative probabilities are given by the following chains of functions:
[0] Z←N pois MU
[1] Z←*(+/(⍟MU)-⍟⍳N)-MU
[0] Z←N pois1 MU
[1] Z←N pois¨MU
[0] Z←N POIS MU
[1] Z←(⊂N)pois1¨MU
8 4⍕⊃0 1 2 3 4 5 POIS 1 2
.3679 .3679 .1839 .0613 .0153 .0031
.1353 .2707 .2707 .1804 .0902 .0361
[0] Z←N poiscdf MU
[1] Z←+/((↑N)+0,⍳-/-N)pois¨MU
[0] Z←N poiscdf1 MU
[1] Z←+\((0,1+¯1⌽1+N),¨N)poiscdf¨MU
[0] Z←N POISCDF MU
[1] Z←(⊂N)poiscdf1¨MU
1000 POISCDF 1000
0.50840937
8 4⍕⊃100 120 140 POISCDF 100 120
.5266 .9773 .9999
.0347 .5243 .9669
The functions LNNCR1, POIS and POISCDF all work with lower rank arguments, so that the subsidiary functions are not required at user level.
Sooner or later, as the size of the integer arguments increases, even the above functions will hit limits, either domain errors, if even the logarithms exceed the largest representable number, or wsfull if the number of integers in the numerator and denominator of the recursive formula become too large.
The first problem is addressed by increasing the base value of the logarithms (the value e is, after all, arbitrary).
The second problem is addressed by splitting the recursively defined fraction into blocks of a size which nearly fills the workspace, and using iteration. (On my computer 10000 is a desirable size.)
When computing very large binomial coefficients, it might well be the case that the number of decimal digits in the result was the primary object of interest, in which case it is appropriate to use:
and the iterative function indicated above is:
[0] Z←n lcr r;t;u;i
[1] Z←n lnncr10 10000⌊ t←r⌊n-r ⋄ i←0
[2] L1:→(t<u←10000×i←i+1)/0
[3] Z←Z++/(10⍟1+n-u)-10⍟u+⍳10000 ⋄ →L1
As an example, the number of decimal digits in 1000000 C 500000 is: