function P = lyngby_cdf_t(x,r)

% lyngby_cdf_t         - Student t distribution
%
%       function P = lyngby_cdf_t(x,r)
%
%       Input:  x   Value where the function is evaluated
%               r   Degrees of freedom
%
%       Output: P   Lower tail area
%
%       Lower tail (cumulated) Student t distribution.
%
%       This function uses an approximation to the t distribution. The
%       approximation cannot be used for all values of x and r. For
%       the case where the approximation cannot be used the t
%       distribution is computed via the incomplete Beta function.
%
%	Ref: Lennart Raade, Beta, 17.4, page 385
%            Press, Numerical Recipes, sec 6.4
%
%       See also LYNGBY, LYNGBY_IDF_T, LYNGBY_CDF_F,
%                LYNGBY_CDF_GAUSS. 
%
% $Id: lyngby_cdf_t.m,v 1.4 2003/02/11 19:55:54 fnielsen Exp $


    % By manual testing the approximation in Raade seems not to
    % produce very accurate values for the distribution function when
    % the dofs is small.

    if r < 6
       % This is almost equivalent to equation 6.4.9 in "Numerical
       % Recipes in C". However, the equation in numerical is a
       % two-tailed valued 
     
       P = 1 - 0.5 * betainc(r/(r+x^2), 0.5 * r, 0.5); 
       return
    end
    
    P = zeros(size(x));
    for n = 1:length(x);      
      if x(n) == 0
	P(n) = 0.5;
      else
	A = (x(n).^(2/3)) + (2/9) - ((2*(x(n).^(2/3)))/(9*r)) - 1;
	B = 2 * (1 + (x(n).^(4/3))/r) / 9;
	
	y = A / sqrt(B);

	if y < 0 | ~isreal(x(n))
	  if x(n) > 0
	    P(n) = 1 - 0.5 * betainc(r ./ (r+x(n).^2), 0.5*r, 0.5);  
	  else
	    P(n) = 0.5 * betainc(r ./ (r+x(n).^2), 0.5*r, 0.5);  
	  end
	else
	  p = y * (1.5976 + 0.070566*y.^2);
	  G = 1/( 1 + exp(-p) );
	  
	  P(n) = (1 + G) / 2;
	end
      end
    end
      
      










