custom_tcdf.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function F = custom_tcdf(x,v)
% custom_tcdf  - CDF of Students t distribution
%
% FORMAT:       p = custom_tcdf(x, v)
%
% Input fields:
%
%       x           t-variate (Student's t has range (-Inf,Inf)
%       v           degrees of freedom (v > 0, non-integer d.f. accepted)
%
% Output fields:
%
%       p           CDF of t distribution with v d.f. at points x
%
% Note: this function has been copied from the SPM2 package published
%       by the Wellcome Department, go here for details:
%       http://www.fil.ion.ucl.ac.uk/spm/

% spm_Tcdf implements the Cumulative Distribution of the Students
% t-distribution.
%
% Definition:
%-----------------------------------------------------------------------
% The CDF F(x) of the Student's t-distribution with v degrees of
% freedom is the probability that a realisation of a t random variable
% X has value less than x; F(x)=Pr{X<x} for X~G(h,c). Student's
% t-distribution is defined for real x and positive integer v (See
% Evans et al., Ch37).
%
% This implementation is not restricted to whole (positive integer) df
% v, rather it will compute for any df v>0.
%
% Variate relationships: (Evans et al., Ch37 & 7)
%-----------------------------------------------------------------------
% The Student's t distribution with 1 degree of freedom is the Standard
% Cauchy distribution, which has a simple closed form CDF.
%
% Algorithm:
%-----------------------------------------------------------------------
% The CDF of the Student's t-distribution with v degrees of freedom
% is related to the incomplete beta function by:
%       Pr(|X|<x) = betainc(v/(v+x^2),v/2,1/2)
% so
%              {     betainc(v/(v+x^2),v/2,1/2) / 2      for x<0
%       F(x) = |   0.5                                   for x=0
%              { 1 - betainc(v/(v+x^2),v/2,1/2) / 2      for x>0
%
% See Abramowitz & Stegun, 26.5.27 & 26.7.1; Press et al., Sec6.4 for
% definitions of the incomplete beta function. The relationship is
% easily verified by substituting for v/(v+x^2) in the integral of the
% incomplete beta function.
%
% MatLab's implementation of the incomplete beta function is used.
%
%
% References:
%-----------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%
%_______________________________________________________________________
% @(#)spm_Tcdf.m	2.2 Andrew Holmes 99/04/26

% Version:  v0.6c
% Build:    7012416
% Date:     Jan-24 2007, 4:10 PM CET
% Author:   Andrew Holmes, SPM2
% URL/Info: http://www.fil.ion.ucl.ac.uk/spm/

% argument check
if nargin < 2 || ...
   ~isa(x, 'double') || ...
   ~isa(v, 'double') || ...
    isempty(v) || ...
    any(isnan(x(:))) || ...
    any(isinf(v(:)) | isnan(v(:)) | v(:) < 0)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(x)
    F = [];
    return;
end
osize = size(x);
if numel(x) == 1 && ...
    numel(v) > 1
    osize = size(v);
end
x = x(:);
v = v(:);
if numel(v) ~= numel(x) && ...
    numel(v) ~= 1 && ...
    numel(x) ~= 1
    error( ...
        'BVQXtools:SizeMismatch', ...
        'Invalid number of elements in v.' ...
    );
end
if numel(x) ~= numel(v)
    if numel(x) == 1
        x = repmat(x, size(v));
    else
        v = repmat(v, size(x));
    end
end

% initialize output
F = zeros(numel(x), 1);

% special case: f is 0.5 when x == 0
F(x == 0) = 0.5;

% special case: Standard Cauchy distribution when v == 1
v1 = v == 1;
F(v1) = 0.5 + atan(x(v1)) / pi;

% compute where defined and no special case
d = find((~v1) & (x ~= 0));
if isempty(d)
    F = reshape(F, osize);
    return
end
xpos = d(x(d) > 0);
if ~isempty(xpos)
    if numel(v) == 1
        F(xpos) = 1 - 0.5 * betainc(v ./ ...
            (v + x(xpos) .^ 2), v / 2, 1 / 2);
    else
        F(xpos) = 1 - 0.5 * betainc(v(xpos) ./ ...
            (v(xpos) + x(xpos) .^ 2), v(xpos) ./ 2, 1 / 2);
    end
end
xpos = find(x(d) < 0);
if ~isempty(xpos)
    if numel(v) == 1
        F(xpos) = 0.5 * betainc(v ./ ...
            (v + x(xpos) .^ 2), v ./ 2, 1 / 2);
    else
        F(xpos) = 0.5 * betainc(v(xpos) ./ ...
            (v(xpos) + x(xpos) .^ 2), v(xpos) ./ 2, 1 / 2);
    end
end
F = reshape(F, osize);