function P = lyngby_cdf_f(x, r1, r2)

% lyngby_cdf_f         - Fisher's F distribution
%
%       function P = lyngby_cdf_f(x, r1, r2)
%
%       Input:  x   Value where the function is evaluated
%               r1  Degrees of freedom
%               r2  Degrees of freedom
%
%       Output: P   Lower tail area
%
%       Lower tail (cumulated) Fisher's F-distribution. The value is 
%       computed by calls to the incomplete beta function.
%
%       Ref: Press et al, Numerical Recipes in C, page 229
%            (Lennart Raade, Beta, 17.2 (p 375) 17.4, page 385)
%
%       See also LYNGBY, LYNGBY_PDF_F, LYNGBY_CDF_BETA,
%                LYNGBY_CDF_CHI2, LYNGBY_CDF_FISHERSG,
%                LYNGBY_CDF_GAUSS, LYNGBY_CDF_T, LYNGBY_CDF_WILRANK.  
%
% $Id: lyngby_cdf_f.m,v 1.2 2003/02/11 19:53:30 fnielsen Exp $


    % Check input arguments
    if nargin < 3, 
      error('This function requires three inputs.'); 
    end
    if any(x < 0)
      error('x cannot be less than zero')
    end
    if any(r1 < 1)
      error('r1 cannot be less than one');
    end
    if any(r2 < 1)
      error('r2 cannot be less than one');
    end

    % Default values
    bBetainc = 1;
    
    
    if bBetainc
      % Evaluate F-distribution with the incomplete beta function
      
      P = zeros(size(x));
      for n = 1:length(x);      
	if x(n) <= 0
	  P(n) = 0;
	else
	  
	  r = r2 ./ (r2+r1 .* x(n));
	  P(n) = 1 - betainc(r, r2/2, r1/2);
	end
      end
      
    else
      % Use an approximation
    
      P = zeros(size(x));
      for n = 1:length(x);      
	if x(n) <= 0
	  P(n) = 0;
	else
	  A = x(n).^(1/3) + 4.5/r1 -  4.5 * x(n).^(1/3)/r2 - 1;
	  B = 4.5 * ( x(n).^(2/3)/r2 + 1/r1);
	  
	  z = A / sqrt(B);

	  if r2 >= 4
	    p = z * (1.5976 + 0.070566*z.^2);
	    P(n) = 1/( 1 + exp(-p) );
	  else
	    x = z + 0.08 * z.^5 / r2.^3;
	    p = x * (1.5976 + 0.070566*x.^2);
	    P(n) = 1/( 1 + exp(-p) );
	  end
	end
      end
      
    end










