## Copyright (C) 2015-2016 Nir Krakauer
##
## This function is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## This function is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License, <http://www.gnu.org/licenses/>, for more details.

## -*- texinfo -*-
## @deftypefn {Function File} {} tcdf_integer_df (@var{x}, @var{n})
## For each element of @var{x}, compute the cumulative distribution
## function (CDF) at @var{x} of the t (Student) distribution with
## @var{n} degrees of freedom, i.e., PROB (t(@var{n}) @leq{} @var{x}).
## @var{n} should be a positive integer; @var{x} may be a real-valued array of any size.
##
## The summation method used is different from that of the standard @code{tcdf} function and makes this function more efficient for large @var{x} arrays and small integer @var{n}.
##
## Reference: Christian Walck (2007), Hand-book on Statistical Distributions for Experimentalists, University of Stockholm Internal Report SUF-PFY/96-01, Section 38.12, url: http://www.fysik.su.se/~walck/suf9601.pdf
## @seealso{tcdf}
## @end deftypefn

## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
## Description: CDF of the t distribution with integer degrees of freedom

function cdf = tcdf_integer_df (x, n)

  if (nargin != 2)
    print_usage ();
  endif

  if (!isscalar (n)) || (n <= 0) || (fix (n) != n)
    error ("tcdf_integer_df: N must be a positive integer");
  endif

  if (iscomplex (x) || iscomplex (n))
    error ("tcdf_integer_df: X and N must not be complex");
  endif

  if (n == 1)
    cdf = 0.5 + atan(x)/pi;
  
  elseif (n == 2)
    cdf = 0.5 + x ./ (2 * sqrt(2 + x .^ 2));
    
  else
    xs = x ./ sqrt(n);
    xxf = 1 ./ (1 + xs .^ 2);
    u = s = 1; 

    if mod (n, 2) #n odd

      m = (n - 1) / 2;      
      for i = 2:m
        u .*= (1 - 1/(2*i - 1)) .* xxf;
        s += u;
      endfor
      cdf = 0.5 + (xs .* xxf .* s + atan(xs)) / pi;      
      
    else #n even

      m = n / 2;
      for i = 1:(m - 1)
        u .*= (1 - 1/(2*i)) .* xxf;
        s += u;
      endfor
      cdf = 0.5 + (xs .* sqrt(xxf) .* s) / 2;

    endif

  endif

endfunction


%!shared x,y
%! x = [-Inf 0 1 Inf];
%! y = [0 1/2 3/4 1];
%!assert (tcdf_integer_df (x, 1), y)

%% Test input validation
%!error tcdf_integer_df ()
%!error tcdf_integer_df (1)
%!error tcdf_integer_df (1, 0)
%!error tcdf_integer_df (1,2,3)
%!error tcdf_integer_df (ones (3), ones (2))
%!error tcdf_integer_df (ones (2), ones (2))
%!error tcdf_integer_df (i, 2)
%!error tcdf_integer_df (2, i)

## check accuracy for small positive and large negative values of x
%!shared x, tol
%! tol = 20 * eps;
%! x = [logspace(-10, -16, 7) -logspace(1, 5, 5)];
%!assert (tcdf_integer_df (x, 2), tcdf (x, 2), tol)
%!assert (tcdf_integer_df (x, 3), tcdf (x, 3), tol)
%!assert (tcdf_integer_df (x, 19), tcdf (x, 19), tol)
%!assert (tcdf_integer_df (x, 20), tcdf (x, 20), tol)
%!assert (tcdf_integer_df (x, 999), tcdf (x, 999), tol)
%!assert (tcdf_integer_df (x, 1000), tcdf (x, 1000), tol)

%{
#compare time needed with the built-in function
N = 3000;
df = 10000;
x = randn(N, "single");
tic; c1 = tcdf (x, df); t1 = toc; #0.8 s with N=1E3, 7.8 s with N=3E3
tic; c2 = tcdf_integer_df (x, df); t2 = toc; #0.08 s with N=1E3, 1.0 s with N=3E3
%}

