Code covered by the BSD License

# ToleranceFactor

### Viktor Witkovsky (view profile)

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)
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);```