Code covered by the BSD License  

Highlights from
ToleranceFactor

ToleranceFactor

by

 

16 May 2009 (Updated )

ToleranceFactor computes the exact tolerance factor for the two-sided tolerance interval

nctcdfVW(x,nu,ncp,tail,method,tol)
function [cdf,z,fun,tail,method,tol] = nctcdfVW(x,nu,ncp,tail,method,tol)
%NCTCDFVW Cummulative distribution function (cdf) of the noncentral
%Student's t-distribution with nu degrees of freedom and noncentrality
%parameter delta. 
%Optionally the upper tail probability is evaluated, 
% i.e. Prob(T_nu(delta) > t). The cdf is computed by direct integration:
%
% If t = 0:
% CDF(t,lowertail) = Prob(T_nu(delta) < t) = normcdf(-delta).
%
% If t > 0:
% CDF(t,lowertail) = Prob(T_nu(delta) < t) = normcdf(-delta) +
% Integral_(-delta)^(+Inf) (1-chi2cdf(nu*(z+delta)^2/t^2)) * normpdf(z) dz.
%
% If t < 0:
% CDF(t,lowertail) = 1 - Prob(T_nu(-delta) < -t).
%
%Syntax:
% cdf = nctcdfVW(t,nu,delta)
% [cdf,z,fun] = nctcdfVW(t,nu,delta,tail,method,tol)
%
%Examples:
% t = 39; nu = 12; delta = 39;
% cdf = nctcdfVW(t,nu,delta)
%
% t = 1000; nu = 3; delta = 1; 
% tail = 'upper'; method = 'Gauss-Kronod';
% [cdf,z,fun] = nctcdfVW(t,nu,delta,tail,method);
% disp(cdf)
% plot(z,fun)
% xlabel('z');ylabel('fun');
% title('(1 - chi2cdf( nu * (z + delta)^2 / t^2 ) ) * normpdf(z)')

%(c) Viktor Witkovsky (witkovsky@savba.sk)
%Ver.: 23-Oct-2009 11:16:29

%% Check the input parameters
if nargin <  3,
    error('VW:nctcdfVW:TooFewInputs','Requires three input arguments.');
end

if nargin < 4
    tail = [];
end

if nargin < 5
    method = [];
end

if nargin < 6
    tol = [];
end

[errorcode x nu ncp] = distchck(3,x,nu,ncp);

if errorcode > 0
    error('VW:nctcdfVW:InputSizeMismatch',...
        'Requires non-scalar arguments to match in size.');
end

if isempty(tail)
    tail = 'best';
end

if isempty(method)
    method = 'Gauss-Kronod';
end

if isempty(tol)
    tol = eps('double');
end

if isnumeric(method) && method == 1
    method = 'Gauss-Kronod';
elseif isnumeric(method) && method ~= 1
    method = 'Gauss-Legendre';
end

%% Initialize the output parameters
if isa(x,'single') || isa(nu,'single') || isa(ncp,'single')
    cdf = zeros(size(x),'single');
else
    cdf = zeros(size(x));
end

z = [];
fun = [];
%% Set the result if x = 0
if x == 0
    if ncp > 0
        cdf = 0.5 * erfc(ncp ./ sqrt(2));
        isLowerTail = true;
    else
        cdf = 0.5 * erfc(-ncp ./ sqrt(2));
        isLowerTail = false;
    end
    switch tail
        case 'lower'
            if ~isLowerTail
                cdf = 1 - cdf;
            end
        case 'upper'
            if isLowerTail
                cdf = 1 - cdf;
            end
        otherwise
            if isLowerTail
                tail = 'lower';
            else
                tail = 'upper';
            end
            warning('VW:nctcdfVW:Tail', ...
                ['The result is the *', tail, '* tail probability !!!']);
    end
    return
end
%% Transform values
if  x < 0
    t = -x;
    delta = -ncp;
    if t > delta
        isLowerTail = true;
        gammaincTail = 'lower';
        cdf0 = 0;
    else
        isLowerTail = false;
        gammaincTail = 'upper';
        cdf0 = 0.5 * erfc(delta ./ sqrt(2));
    end
elseif x > 0
    t = x;
    delta = ncp;
    if t > delta
        isLowerTail = false;
        gammaincTail = 'lower';           
        cdf0 = 0;
    else
        isLowerTail = true;
        gammaincTail = 'upper';
        cdf0 = 0.5 * erfc(delta ./ sqrt(2));
    end
end

%% Set the integraltion limits [A,B]
%TODO: find a better way to set the limits
%MAXUPPER = -norminv(eps(0)) = 38.45

% if nu > 2
%     gammaratio = exp(gammaln((nu-1)/2) - gammaln(nu/2));
%     meanNCT = delta * sqrt(nu/2) * gammaratio;
%     varNCT = nu * (1+ delta^2) / (nu - 2) - delta^2 * nu * gammaratio / 2;
% end

MAXUPPER = 38.45;
UPPER = -norminv(tol);
A = max(-delta,-UPPER);
B = UPPER;

while B < A && B < MAXUPPER
    B = B + UPPER/2;
end

if B < A
    A = B;
end

while A/B > 0.8 && B < MAXUPPER
    B = B + UPPER/A;
end

%% Integrate
switch method
    case 'Gauss-Kronod'
        for i = 1:length(t)
            cdf(i) = quadgk(@(z) Fun(z,t(i),nu(i),delta(i),gammaincTail),...
                A(i),B(i),'AbsTol',tol);
        end
    case 'Gauss-Legendre'
        [bp,wf]=grule(15);
        [z,w] = gweights([A,B],25,bp,wf);
        for i = 1:length(t)
            cdf(i) = w' *  Fun(z,t(i),nu(i),delta(i),gammaincTail);
        end
    otherwise
        [bp,wf]=grule(15);
        [z,w] = gweights([A,B],25,bp,wf);
        for i = 1:length(t)
            cdf(i) = w' *  Fun(z,t(i),nu(i),delta(i),gammaincTail);
        end
end
cdf = cdf + cdf0;
%% Set the result
switch tail
    case 'lower'
        if ~isLowerTail
            cdf = 1 - cdf;
        end
    case 'upper'
        if isLowerTail
            cdf = 1 - cdf;
        end
    otherwise
        if isLowerTail
            tail = 'lower';
        else
            tail = 'upper';
        end
        warning('VW:nctcdfVW:Tail', ...
            ['The result is the *', tail, '* tail probability !!!']);
end

z = linspace(A,B)';
fun = Fun(z,t,nu,delta,gammaincTail);
%% Function Fun (Integrand for the Gauss-Kronod quadrature)
function fun = Fun(z,t,nu,delta,gammaincTail)
%Fun evaluates the Integrand

%(c) Viktor Witkovsky (witkovsky@savba.sk)
%Ver.: 23-Oct-2009 11:16:29

x = nu * (z + delta).^2 / t^2;
fun = gammainc(x/2, nu/2, gammaincTail) .* exp(-0.5 * z.^2) / sqrt(2*pi);

Contact us