Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

nc_t_cdf(x,nu,delta)
function p = nc_t_cdf(x,nu,delta)
%NC_T_CDF Noncentral T cumulative distribution function (cdf).
%   P = NC_T_CDF(X,NU,DELTA) Returns the noncentral T cdf with NU 
%   degrees of freedom and noncentrality parameter, DELTA, at the values 
%   in X.
%
%   The size of Y is the common size of the input arguments. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.     

%   References:
%      [1]  Johnson, Norman, and Kotz, Samuel, "Distributions in
%      Statistics: Continuous Univariate Distributions-2", Wiley
%      1970 p. 205.
%      [2]  Evans, Merran, Hastings, Nicholas and Peacock, Brian,
%      "Statistical Distributions, Second Edition", Wiley
%      1993 pp. 147-148.

%   B.A. Jones 8-22-94
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.13 $  $Date: 2000/05/26 18:53:06 $

if nargin <  3, 
    error('Requires three input arguments.'); 
end

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

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

% Special cases for delta==0 and delta<0.
f0 = (delta == 0);
fn = (delta < 0);
flag0 = any(any(f0));
flagn = any(any(fn));
if (flag0 | flagn)
   fp = ~(f0 | fn);
   p = zeros(size(delta));
   if flag0,        p(f0) = t_cdf(x(f0),nu(f0)); end
   if any(any(fp)), p(fp) = nc_t_cdf(x(fp), nu(fp), delta(fp)); end
   if flagn,        p(fn) = 1 - nc_t_cdf(-x(fn), nu(fn), -delta(fn)); end
   return
end

% Initialize P to zero.
[m,n] = size(x);
p = zeros(m,n);

% Probability that x < 0
p0 = norm_cdf(-delta,0,1);

kneg = find(x < 0);
kzero = find(x == 0);
kpos = 1:m*n;
kpos([kneg(:); kzero(:)]) = [];

%Value passed to Incomplete Beta function.
tmp = (x.^2)./(nu+x.^2);

% Set up for infinite sum.
done = 0;
j = 0;
jk9a  = zeros(size(x));
jk10a = zeros(size(x));
eterm = exp(-0.5*delta.^2);
seps = sqrt(eps);
qeps = eps^(1/4);

% Sum the series.
k0      = find(x~=0);
if any(k0)
   tmp = tmp(k0);
   nu = nu(k0);
   ld = log(0.5*delta(k0).^2);
   jk9 = zeros(size(k0));
   jk10 = zeros(size(k0));
   while 1
      ibeta9 = betainc(tmp,(j+1)/2,nu/2);    % Johnson & Kotz eqn. 9
      ibeta10 = betainc(tmp,(j+1/2),nu/2);    % Johnson & Kotz eqn. 10
      
      %Probability that t lies between 0 and x (x>0)
      nt9   = (exp(0.5*j*ld - gammaln(j/2+1))).*ibeta9;
      jk9   = jk9 + nt9;
      %Probability that t lies between x and -x.
      nt10  = (exp(j*ld - gammaln(j+1))).*ibeta10;
      jk10  = jk10 + nt10; 
      
      % Convergence test.
      if all((nt9./(jk9+qeps)) < seps) & all((nt10./(jk10+qeps)) < seps)
         jk9a(k0) = jk9;
         jk10a(k0) = jk10;
         break;
      end
      j = j+1;
   end
end

% Compute probability for non-negative Xs.
p = p0 + eterm.*jk9a/2;

% Compute probability for negative Xs. P(t < 0) + P(0 < t < |x|) - P(-|x| < t < |x|)
p(kneg) = p(kneg) - eterm(kneg).*jk10a(kneg);


% Return NaN if X is negative or NU is not a positive integer.
p(nu <= 0) = NaN;

return

Contact us