FISHER-SNEDECOR F DISTRIBUTION

[shakuntalaParticipant]

% fp.m
% @shortDesc: F Distribution (Fisher-Snedecor). Script sent by M.Urrutia Avisrror MD,PhD
% @fullDesc: Uses an armonic expansion series to compute the exact probability of the F Distribution with v1 and v2 degrees of freedom [P(X)>=x], (upper tail). Ref. Abramowit and Stegun. Handbook of Mathematical Functions. Dover Publications Inc. NY

function fp = fprob(f,v1,v2)
even1 = 0;
even2 = 0;
if round(v1/2)== v1/2
even1 = 1;
elseif round(v2/ 2) == v2/2
even2 = 1;
end
% solve for v1 even or v2 even
if ((even1 == 1) | (even2 == 1))
X = v2 / (v2 + v1 * f);
if even1 == 1
n1 = v2;
n2 = 2;
L = 1;
M = 1;
while n2 <v1
M = M * n1 / n2 * (1 – X);
L = L + M;
n1 = n1 + 2;
n2 = n2 + 2;
end
fp = power(X ,(v2 / 2)) * L;
end
if even2 == 1
n1 = v1;
n2 = 2;
L = 1;
M = 1;
while n2 <v2
M = M * n1 / n2 * X;
L = L + M;
n1 = n1 + 2;
n2 = n2 + 2;
end
fp = 1 – power((1 – X) , (v1 / 2)) * L;
end
% solves for v1 and v2 uneven
elseif (even1 == 0) & (even2 == 0)
X = atan(sqrt((v1 / v2) * f));

n1 = 0;
n2 = 1;
M = cos(X);
L = cos(X);
while n2 <(v2 – 2)
n1 = n1 + 2;
n2 = n2 + 2;
M = M * (n1 / n2 * power(cos(X),2));
L = L + M;
end
A = 2 / pi * (X + sin(X)* L);
if v2 == 1
A =( 2 * X / pi);
end
fact1 = 1;
num1 = (v2 – 1) / 2;
while num1 >0
fact1 = fact1 * num1;
num1 = num1 – 1;
end
fact2 = 1;
num2 = (v2 – 2) / 2;
while num2 >0
fact2 = fact2 * num2;
num2 = num2 – 1;
end
%adjustment for fractionary factorial number
fact2 = fact2 * sqrt(pi) * 0.5%;
R = (2 / sqrt(pi)) * fact1 / fact2;
n1 = v2 + 1;
n2 = 3;
M = 1;
L = 1;
while n2 <v1
M = M * n1 / n2 * power(sin(X), 2);
L = L + M;
n1 = n1 + 2;
n2 = n2 + 2;
end
B = (R * sin(X)/2*power(cos(X),v2)) * L;
if v1 == 1
B = 0;
end
fp = 1 – A + B;
% use fp=1-fp to compute the CDF, P(X)<=x (lower tail)
end
end

[Author]

Posts