## 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, , 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 ## 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 %}