F-distribution
In probability theory and statistics, the F-distribution is a continuous probability distribution. The F-distribution arises frequently as the null distribution of a test statistic.
An F-test is any statistical test in which the test statistic has an F-distribution under the null hypothesis. A common use of F-tests is in the analysis of variance (ANOVA).
The cumulative F-distribution can be obtained using the regularized incomplete beta function with the following inputs: X = d1 * F / (d1 * F + d2), a = d1 / 2, and b = d2 / 2, where F is the value obtained in the F-test, d1 is the degrees of freedom for the numerator, and d2 is the degrees of freedom for the denominator.
The source code below applies that method to estimate the probability of the null hypothesis for a given F-value.
[code language="vb"] #COMPILE EXE #REGISTER NONE #DIM ALL ' FUNCTION PfromF(BYVAL F AS EXT,BYVAL dfn AS EXT,BYVAL dfd AS EXT) AS EXT LOCAL ifault AS LONG, Beta, betain AS EXT IF F > 0.0 AND dfn > 0.0 AND dfd > 0.0 THEN CALL IncomplBeta(dfd/(dfd + dfn * F),dfd * 0.5, dfn * 0.5, Beta, betain) FUNCTION = betain ELSE FUNCTION = 1.0 END IF END FUNCTION ' FUNCTION PBMAIN LOCAL Result AS STRING, i1, i2 AS LONG LOCAL P, F, dfn, dfd AS EXT Result$ = INPUTBOX$("Enter F-value, dfn, dfd","Get p-value") i1 = INSTR(Result$, ",") i2 = INSTR(i1 + 1,Result, ",") F = VAL(LEFT$(Result, i1 - 1)) dfn = VAL(MID$(Result, i1 + 1,i2 - i1 - 1)) dfd = VAL(MID$(Result,i2 + 1)) P = PfromF(F, dfn, dfd) MSGBOX "P=" + FORMAT$(P," ##.##################"),,"Result:" END FUNCTION ' SUB IncomplBeta(BYVAL X AS EXT, BYVAL P AS EXT, BYVAL Q AS EXT, _ BYREF beta AS EXT, BYREF betain AS EXT) ' ' Derived from FORTRAN code based on: ' algorithm AS 63 Appl. Statist. (1973), vol.22, no.3 ' Computes incomplete beta function ratio for arguments ' X between zero and one, p (=a) and q (=b) positive. ' Returns log beta and Regularized Incomplete Beta ' LOCAL ns, indx AS LONG LOCAL zero, one, acu AS EXT LOCAL psq, cx, xx, pp, qq, term, ai, rx, temp AS EXT ' define accuracy and initialise zero = 0.0 : one = 1.0 : acu = 1.0E-18 betain = 0.0 : beta = 0.0 ' test for admissibility of arguments IF(p <= zero OR q <= zero) THEN EXIT SUB IF(x < zero OR x > one) THEN EXIT SUB IF(x = zero OR x = one) THEN EXIT SUB ' calculate log of beta by using function GammLn beta = gammln(p) + gammln(q) - gammln(p + q) betain = x ' change tail if necessary psq = p + q cx = one - x IF (p < psq * x) THEN xx = cx cx = x pp = q qq = p indx = 1 ELSE xx = x pp = p qq = q indx = 0 END IF term = one ai = one betain = one ns = qq + cx * psq ' use Soper's reduction formulae. rx = xx / cx temp = qq - ai IF (ns = 0) THEN rx = xx DO term = term * temp * rx / (pp + ai) betain = betain + term temp = ABS(term) IF(temp <= acu AND temp <= acu * betain) THEN EXIT DO ai = ai + one ns = ns - 1 IF (ns >= 0) THEN temp = qq - ai IF (ns = 0) THEN rx = xx ELSE temp = psq psq = psq + one END IF LOOP ' calculate Regularized Incomplete Beta betain = betain * EXP(pp * LOG(xx) + (qq - one) * LOG(cx) - beta) / pp IF indx = 1 THEN betain = one - betain END SUB ' FUNCTION GammLn(BYVAL x AS EXT) AS EXT ' Returns Ln(Gamma()) or 0 on error ' Based on Numerical Recipes gamma.h ' Lanczos, C. 1964, "A Precision Approximation ' of the Gamma Function," SIAM Journal on Numerical ' Analysis, ser. B, vol. 1, pp. 86-96. LOCAL j AS LONG, tmp, y, ser AS EXT DIM cof(0 TO 13) AS LOCAL EXT cof(0) = 57.1562356658629235 cof(1) = -59.5979603554754912 cof(2) = 14.1360979747417471 cof(3) = -0.491913816097620199 cof(4) = 0.339946499848118887e-4 cof(5) = 0.465236289270485756e-4 cof(6) = -0.983744753048795646e-4 cof(7) = 0.158088703224912494e-3 cof(8) = -0.210264441724104883e-3 cof(9) = 0.217439618115212643e-3 cof(10) = -0.164318106536763890e-3 cof(11) = 0.844182239838527433e-4 cof(12) = -0.261908384015814087e-4 cof(13) = 0.368991826595316234e-5 IF x <= 0.0 THEN FUNCTION = 0.0 : EXIT FUNCTION ' Bad argument y = x tmp = x + 5.2421875 tmp = (x + 0.5) * LOG(tmp) - tmp ser = 0.999999999999997092 FOR j = 0 TO 13 y = y + 1 ser = ser + cof(j)/y NEXT j FUNCTION = tmp + LOG(2.5066282746310005 * ser / x) END FUNCTION [/code]