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_pdf(x,v,delta)
function y = nc_t_pdf(x,v,delta)
%NC_T_PDF Noncentral T probability density function (pdf).
%   Y = NC_T_PDF(X,V,DELTA) Returns the noncentral T pdf with V 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.     

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

%   B.A. Jones 4-21-95
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.14 $  $Date: 2000/05/26 18:53:08 $


if nargin < 1, 
    error('Requires at least one input argument.');
end

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

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

% if some delta==0, call tcdf for those entries, call nctcdf for other entries.
f = (delta == 0);
if any(f(:))
   y = zeros(size(delta));
   y(f) = t_pdf(x(f),v(f));
   f1 = ~f;
   if any(f1)
      y(f1) = nc_t_pdf(x(f1),v(f1),delta(f1));
   end
   return
end

%   Initialize Y to zero.
y = zeros(size(x));

k = find(v > 0);
if any(k)
   nu = v(k);
   d = delta(k);
   t = x(k);
   lnu = log(nu+t.^2);
   lnu1 = gammaln((nu+1)/2);
   lnu2 = gammaln((nu+2)/2);
   term = exp(-0.5.*(log(pi)+log(nu)) -(d.^2)/2  + lnu1 - gammaln(nu/2) ...
      + ((nu+1)/2).*(log(nu)-lnu));
   
   infsum = zeros(size(k));
   qeps = eps^(3/4);
   tsq = t.^2;
   ratio = tsq .* (d.^2) ./ (nu + tsq);

   % Set up for infinite sum by computing first two terms.
   j = 0;
   t0 = 1;
   t1 = t .* d .* sqrt(2) .* exp(lnu2 - lnu1) ./ sqrt(nu + tsq);
   infsum = infsum + term .* (t0 + t1);
   
   % Sum the series by computing term j from term j-2
   while 1
      t0 = t0 .* ratio .* (nu+j+1) / ((j+1) * (j+2));
      j = j + 1;
      t1 = t1 .* ratio .* (nu+j+1) / ((j+1) * (j+2));
      j = j + 1;
      twoterms = term .* (t0 + t1);
      infsum = infsum + twoterms;

      % Convergence test.
      if all(abs(twoterms)./infsum < qeps), break; end
   end
   y(k) = infsum;
end

p(v <= 0) = NaN;

Contact us