--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/Fcdf.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,48 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Compute cumulative F distribution function
+% 
+% CALL 
+% 
+% p = Fcdf(x,n1,n2);
+% 
+% 
+% INPUT
+% 
+% - x, probability
+% - n1, degree of freedom 1
+% - n2, degree of freedom 2
+% 
+%   References:
+%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
+%      Functions", Government Printing Office, 1964, 26.6.
+% 
+%
+% L Ferraioli 06-12-2010
+%
+% $Id: Fcdf.m,v 1.1 2010/12/06 19:14:36 luigi Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function p = Fcdf(x,n1,n2)
+
+x = n2./(n2+x.*n1);
+p = betainc(x, n2/2, n1/2, 'upper');
+
+% % Compute P when X > 0.
+% k = find(x > 0 & isfinite(n1) & isfinite(n2));
+% if any(k)
+%   k1 = (n2(k) <= x(k).*n1(k));
+%   % use A&S formula 26.6.2 to relate to incomplete beta function
+%   % Also use 26.5.2 to avoid cancellation by subtracting from 1
+%   if any(k1)
+%     kk = k(k1);
+%     xx = n2(kk)./(n2(kk)+x(kk).*n1(kk));
+%     p(kk) = betainc(xx, n2(kk)/2, n1(kk)/2,'upper');
+%   end
+%   if any(~k1)
+%     kk = k(~k1);
+%     num = n1(kk).*x(kk);
+%     xx = num ./ (num+n2(kk));
+%     p(kk) = betainc(xx, n1(kk)/2, n2(kk)/2,'lower');
+%   end
+end