PRBF

Function PRBF

This is a function to find the associated probability of a value (FR) of an

F Ratio

With DA representing Degrees of Freedom Above

And DB representing Degrees of Freedom Below

It returns -999 if any input parameter is less than or equal to zero.

It can also be used to evaluate

Chi-Squared

With DA representing the degrees of freedom

DB set to as large number as possible (999 is probably sufficient)

FR set to the chi-squared value divided by the degrees of freedom DA

Student’s “t”

With DA set to unity (1)

DB set to the degrees of freedom

Fr set to the value of t squared

A normal variate

With DA set to Unity

DB set to as large a number as possible

FR set to the square of the normal variate

(The function GAUS, which is called by PRBF, gives this probability directly.)

It is based on

Jaspen, N., “The calculation of Probabilities corresponding to Values of Z, t, F, and Chi-squared” Education and Psychological Measurement, 1965 Vol 25, pp 887-880

This function is written in ASA FORTRAN IV – but it should be fairly easy to translate into  a more modern language.

Enjoy

 

FUNCTION PRBF(DA,DB,FR)

C PRBF   (SEE 11.8) STATISTICAL FINDS PROBABILITY OF F-RATIO

C DA = DF ABOVE,DB=DF BELOW,FR=F RATIO – PRBF=-999 =OUT OF RANGE

PRBF = -999.

IF(DA.LE.0.0.OR.DB.LE.0.0.OR.FR.LT.0.0) RETURN

PRBF = 1.

IF(FR.EQ.0.) RETURN

LI = 500

IF(DA.EQ.1.0.AND.DB.GT.LI) X = SQRT(FR)

IF(DA.EQ.1.0.AND.DB.GT.LI) GO TO 15

FT = 0.

PRBF = 0.

X = DB/(DB + DA*FR)

NA = DA + 0.5

NB = DB + 0.5

VP = NA + NB – 2

FA = NA

FB = NB

IF(NA.EQ.NA/2*2.AND.NA.LT.LI) GO TO 1

IF(NB.EQ.NB/2*2.AND.NB.LT.LI) GO TO 4

IF(NA+NB.LT.LI) GO TO 6

GO TO 14

C FA EVEN AND LESS THAN LI (THE LIMIT FOR ITERATION)

1 XX = 1.-X

2 FA = FA – 2.

IF(FA.LT.1.) GO TO 3

VP = VP – 2.

FT = XX*VP/FA*(FT+1.)

GO TO 2

3 FT = X**(.5*FB)*(FT+1.)

IF(FT.GT.0.)PRBF = FT

RETURN

C FB EVEN AND LESS THAN LI (THE LIMIT FOR ITERATION)

4 FB = FB – 2.

IF(FB.LT.1.) GO TO 5

VP = VP – 2.

FT = X*VP/FB*(FT+1.)

GO TO 4

5 FT = 1. – (1.-X)**(0.5*FA)*(FT+1.)

IF(FT.GT.0.)PRBF = FT

RETURN

C FA AND FB ODD – FA+FB LESS THAN LIMIT FOR ITERATION

6 TH = ATAN(SQRT(FA*FR/FB))

ST = SIN(TH)

CT = COS(TH)

C2 = CT*CT

S2 = ST*ST

A = 0

B = 0

IF(DB.LE.1.) GO TO 9

7 FB = FB – 2.

IF(FB.LT.2.) GO TO 8

A = C2*(FB-1.)/FB*(1.+A)

GO TO 7

8 A = ST*CT*(1.+A)

9 A = TH + A

IF (DA.LE.1.) GO TO 13

10 FA = FA – 2.

IF (FA.LT.2.) GO TO 11

VP = VP – 2.

B = S2*VP/FA*(1.+B)

GO TO 10

11 GF = 1.

K = DB*0.5

DO 12 I = 1,K

X1 = I

12 GF = X1*GF/(X1-0.5)

B = GF*ST*CT**DB*(1.+B)

13 FT = 1. + 0.636619772368*(B-A)

IF(FT.GT.0.)PRBF = FT

RETURN

C FA AND/OR FB EXCEED LIMITS

14 FA = 2./9./FA

FB = 2./9./FB

CR =FR**0.333333333333

X = (1.-FA +(FB-1.)*CR )/SQRT(FB*CR*CR + FA)

15 FT = GAUS(X)

IF(FT.GT.0.) PRBF = FT

RETURN

END

GAUS

FUNCTION GAUS

 

Function to derive the associated probability of a value on a Normal distribution with Arithmetic mean 0 and Standard Deviation 1

From my book STATCAT – no home should be without one.

This function is based on based on a numerical approximation, since the normal distribution cannot be integrated. (Don’t ask me why, I only work here.)

The numerical approximation was based on an ALGOL program by D. Ibbotsen (ACM Collected Al;gorithms No. 209 , 1963). The original ALGOL formulation has been modified to reduce the rounding errors associated with a negative value of X.(If you estimate the size of the tail of the distribution as, say, 3.141 x 10-9 and subtract this from 1.0 you will get a row of 9s or even 1.0, which doesn’t help much.)

Most computers will not cope with all twelve significant figures, but someone might need them.

This function is written in ASA FORTRAN IV – originally punched on cards, the lines beginning 1,2,3,4,5 are extensions of the previous line, if anyone cares.

 

 

FUNCTION GAUS(X)

C GAUS   (SEE 11.8 ) STATISTICAL FINDS ASS. PROB. OF N VALUE

IF(X.EQ.0.) GAUS = 0.5

IF(X.LT.-6.)GAUS = 0.0

IF(X.GT.6.) GAUS = 1.0

Y =ABS(X)*0.5

IF(X.EQ.0..OR.Y.GT.3.)RETURN

IF(Y.GT.1.) GO TO 15

W = Y*Y

Z =((((((((0.000124818987 *W – 0.001075204047)*W + 0.005198775019)

1  * W    – 0.0191982914)*W + 0.059054035642)*W – 0.151968751364)

2  * W    + 0.319152932694)*W – 0.53192171)*W + 0.797884560593)

3  * Y

IF(X.LE.0.) GAUS = 0.5 – Z

IF(X.GT.0.) GAUS = 0.5 + Z

RETURN

15 Y = Y – 2.

Z = (((((((((((((-0.000045255659*Y + 0.000152529290)*Y

1 – 0.000019538132)*Y – 0.000676904986)*Y + 0.001390604284)*Y

2 – 0.0007948820)*Y – 0.002034254874)*Y + 0.006549791214)*Y

3 – 0.010557625006)*Y + 0.0119447319)*Y – 0.009279453341)*Y

4 + 0.005353579108)*Y – 0.002141268741)*Y + 0.0005353849)*Y

5 – 0.000063342476

IF(X.LE.0.) GAUS = – Z * 0.5

IF(X.GT.0.) GAUS = 1.+ Z*0.5

RETURN

END