Numeric Integration & Probability
Vector 1:2, page 117
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 117
At school we learned a lot of clever tricks to integrate a wide class of functions \(f(x)\) to get a definite integral
Unfortunately, there is an even wider class of functions which occur frequently in all branches of science and engineering for which no algebraic solution can be found in terms of known simple functions. This is certainly the case in statistics, where if \(f(x)\) is a probability density function the above integral represents the probability of an observation falling in the interval \((a,b)\). In such cases one must resort to numerical integration techniques to get numerical answers. Often the results are tabulated for users to look up, such as the normal, \(t\), chisquared, and \(F\) distributions which are so commonly used in statistics and found in e.g. Biometrika Tables by Pearson and Hartley (1958).
A lot of ingenious methods were developed to get good approximations with as little calculation as possible, using analytical techniques to place error bounds on the results. Biometrika and other tables were computed to very good accuracy, with great skill and dedication, using the famous Brunsviga hand calculating machines. With modern computers, especially using APL, it is so easy to calculate the function over a large grid of points that simple techniques will be perfectly adequate for most purposes, and accuracy may be controlled simply by repeating the calculation over a finer grid to see if results differ by less than some prescribed accuracy; if so, stop, otherwise make the grid finer again.
The simplest well known technique, the trapezium rule, approximates the function by series of straight lines. Very nearly as simple, and usually quite a bit more accurate for the same grid size is Simpson’s rule. If the interval \((a,b)\) is divided up into \(2M\) intervals with division points at
the function \(f(x)\) is approximated by quadratic curves, rather than straight lines, over each pair of intervals \((x_0,x_2)\) \((x_2,x_4)\ldots\) . This leads to the approximate integral
The APL function Y NUMINT X enables one very easily to carry out numerical integration using Simpson’s rule for any well behaved function (one which is unbounded, for example, would need special treatment). The full specification of the right argument is a three element vector X←A,B,N. The result is then an N-vector containing the integral of the function over each of N equal subintervals into which the interval (A,B) is divided; e.g. if X←1 5 2 the result consists of the integral over (1,3) and the integral over (3,5). If N is omitted, so that X←A,B, it is automatically assumed that N←1 and the result is simply the integral over the whole range (A,B). If X contains a single item only, X←B we assume A←0 and N←1 so that the result is the integral over (0,B).
The function that is to be integrated is specified by the left argument Y in the form of a string vector. The usual form taken by Y will be the name of a function with an optional left argument to specify the shape if required. The function is expecting a vector as its right argument and returns as a result a vector of the same shape, containing the value of the function at each point of the input vector. The string is catenated with an appropriate 'GRID' in line [21] and executed, ⍎ ; this is also done with a different value for GRID in line [23]. In general one thus has a statement of the form 'A functionname' NUMINT X, with A a shape parameter vector.
2 Examples Using the Normal Distribution¶
The probability density function, pdf, of the standard normal distribution is given by
This is calculated for each element of the vector X by the function NORMALPDF X. As an example the probability attached to the interval (1.5,2.4) say is
while
is the probability over the interval (0, 1.96). The probabilities for the three intervals (0,1), (1,2), (2,3) can be found from
For the normal distribution the probability of the random variable X being less than a given value \(x > 0\). can be written as \(P(0 < X < x) + \tfrac12\). This can be computed for each \(x = 1, 2, 3, 4\) say by
which agrees exactly with results in Table 1 of the Biometrika Tables.
The accuracy of integration is controlled simultaneously for all intervals over which the integral is required by a global variable ∆INTERR which was set at 1E¯7 for all the above examples. Suppose now we go right up into the tail and require the probability for each of the intervals (14,16), (16,18) and (18,20)
Now these are very much smaller than the accuracy setting, so we should not be surprised if these were not too good. If now we let
That’s much better, these are known to be correct to four figures.
In principle it would be a good idea to make ∆INTERR smaller still to cope with some of the very small values. However, with ∆INTERR←1E¯60 this soon leads to WS FULL! The reason is that at this level of accuracy the grid has to be very fine indeed over intervals for which the integral is relatively large. The moral of this is you should not simultaneously use intervals over which the values of the integral vary by a very large amount, and that ∆INTERR should not be made much smaller than you really need it. In isolation, with ∆INTERR←1E¯80
works perfectly well. The WS FULL problem previously was obviously due to the first interval.
These calculations were done on a PRIME 750. If you are working on a much smaller machine you would run out of space quicker. In that case you may prefer to take a different approach and approximate the integral by a suitable algebraic function. In the case of normal distribution this is done in Hawkes (1982) and is implemented in the function NORMTAIL X. This is known to be accurate to at least four significant figures, so for example compute
(¯1↓Z) - 1↓Z← NORMTAIL 14 16 18 20
7.793E¯45 6.388E¯58 9.741E¯73
)DIGITS 7
1 - NORMTAIL ⍳4
0.8413451 0.9772489 0.9986502 0.9999683
These should be compared with results (2) and (1) respectively which are accurate to the number of figures shown. The approximation in NORMTAIL is good enough for most purposes and requires less workspace than using NUMINT (although it is not faster).
3 Some More Distributions¶
The shape of the probability density of student’s \(t\) distribution depends on a parameter known as the degrees of freedom, DF. This is calculated at a vector of points X by the function DF TPDF X. Then for the \(t\) distribution with 5 degrees of freedom say we can find
which agrees exactly with the results in Biometrika Table 9. Again exact agreement with Biometrika Table 7 is found from the chisquared distribution with 12 degrees of freedom.
which are the probabilities of exceeding 6, 12, 18, and 24 respectively.
The chisquared distribution is a special case \((\alpha = \mathit{DF}/2,\ \lambda = \tfrac12)\) of the gamma distribution which has pdf
This is computed in A GAMMAPDF X where A is a two element vector containing values of \(\alpha\) and \(\lambda\).
Note that if \(\alpha<1\), you need to be careful not to include zero in any interval to be integrated over because the function tends to infinity. If \(\alpha = 1\) we have the well-known negative exponential distribution
Of course, we do not really need to use numerical integration on this one, but it can be done very simply, e.g.
'0.5×*-0.5×' NUMINT 0 8 4
0.63212 0.23254 0.08555 0.03147
'1 0.5 GAMMAPDF' NUMINT 0 8 4
0.63212 0.23254 0.08555 0.03147
of course gives exactly the same result. The former version, however, is interesting because it is done with only raw APL in the left argument.
Our range of basic probability distributions is completed by the \(F\) distribution and the Beta distribution, each having a left argument which is a vector of 2 elements.
in agreement with p.44 of Soper (1921), evalutes the integral
for \(a = 325\), \(b = 577\). Again note that there may be difficulties at \(x = 0\) and/or \(x = 1\) if \(a\) and/or \(b\) is less than unity
Of course \(a\), \(b\) are not restricted to integers
agrees with the result on p.17 of Soper.
References¶
E. S. PEARSON and O. HARTLEY (1958), Biometrika Tables for Statisticians Vol.1, Cambridge University Press.
A. G. HAWKES (1982), Approximating the normal tail, The Statistician, 31, 231-236.
H. E. SOPER (1921), Tracts for computers No. VII, Cambridge University Press.
Function listings¶
∇ R←Y NUMINT X;N;RANGE;RT;DX;GRID;FO;FE;M;S
[1] ⍝ NUMERICAL INTEGRATION. <Y> IS A STRING OF FORM 'A FNNAME'
[2] ⍝ THE GLOBAL <∆INTERR> SPECIFIES ACCURACY
[3] ⍝ RESULT IS INTEGRAL OVER <N> SUBINTERVALS OF X[1],X[2]
[4] XCHECK:→((1≥⍴⍴X)∧(3≥⍴,X))/XOPTIONS
[5] XMESS:'X SHOULD BE VECTOR OF SHAPE 3. X[3] IS NUMBER OF '
[6] 'INTERVALS INTO WHICH (X[1],X[2]) DIVIDED'
[7] 'DEFAULT X[3] IS 1 IF NOT SPECIFIED'
[8] →0,0⍴⎕←'IF ONLY SCALAR X GIVEN THEN IT IS TAKEN AS X←0,X,1'
[9] XOPTIONS:→(1<⍴,X)/TWO
[10] X←0,X
[11] TWO:→(3=⍴X)/XCHECK2
[12] X←X,1
[13] XCHECK2:→(X[3]≠1⌈⌊X[3])/XMESS
[14] M←1
[15] R←(N←X[3])⍴0
[16] RANGE←X[2]-X[1]
[17] LOOP:RT←R
[18] M←M×2
[19] DX←RANGE÷S←M×N
[20] GRID←X[1]+DX×0,⍳S
[21] FO←⍎Y,' GRID'
[22] GRID←X[1]+DX×(⍳S)-0.5
[23] FE←⍎Y,' GRID'
[24] R←(DX÷6)×+/(N,M)⍴(4×FE)+(1↓FO)+¯1↓FO
[25] →(∆INTERR<⌈/|R-RT)/LOOP
∇
∇ R←NORMALPDF X
[1] R←(÷(○2)*0.5)×*-0.5×X*2
∇
∇ R←DF TPDF X
[1] R←((!0.5×DF-1)÷(!¯1+DF÷2)×(○DF)*0.5)÷(1+(X*2)÷DF)*0.5×1+DF
∇
∇ R←DF FPDF X
[1] R←(!¯1+0.5×+/DF)×((÷/DF)*0.5×1↑DF)÷×/!¯1+DF÷2
[2] R←R×(X*¯1+0.5×1↑DF)÷(1+X×÷/DF)*+/DF÷2
∇
∇ R←A BETAPDF X
[1] ⍝ PDF OF BETA DISTRIBUTION
[2] R←((!1++/A)÷×/!A)×(X*1↑A)×(1-X)*1↓A←A-1
∇
∇ R←A GAMMAPDF X;L
[1] ⍝ <A> IS ALPHA,LAMDA
[2] R←(X*A)×(*-L×X)÷(!A←1↑A-1)÷(L←1↓A)*1↑A
∇
∇ R←DF CHISQPDF X
[1] R←((DF÷2),0.5)GAMMAPDF X
∇
∇ R←NORMTAIL X;Z;T;Q1;Q2
[1] ⍝ UPPER TAIL N(0,1) FOR ARRAY <X>. HAWKES 1982-STATISTICIAN
[2] Z←|X
[3] T←Z×Z≤2.2
[4] T←T×1+T×Tׯ7.5166E¯3+T×T×3.174E¯4-2.965E¯6×T×T
[5] Q1←0.5×1-(0⌈1-*-T×T÷○0.5)*0.5
[6] Q2←(1.184+Z)×(*-Z×Z÷2)÷((1.209+Z×1.176+Z)×(○2)*0.5)
[7] T←Z≤2.2
[8] R←(T×Q1)+(1-T)×Q2
[9] T←X≥0
[10] R←(T×R)+(1-T)×(1-R)
∇