Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

generates random variates from over 870 univariate distributions

nc_f_cdf(x,nu1,nu2,delta)
function p = nc_f_cdf(x,nu1,nu2,delta)
%NC_F_CDF Noncentral F cumulative distribution function (cdf).
%   P = NC_F_CDF(X,NU1,NU2,DELTA) Returns the noncentral F cdf with numerator 
%   degrees of freedom (df), NU1, denominator df, NU2, and noncentrality
%   parameter, DELTA, at the values in X.
%
%   The size of P 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. 192.
%      [2]  Evans, Merran, Hastings, Nicholas and Peacock, Brian,
%      "Statistical Distributions, Second Edition", Wiley
%      1993 p. 73-74.

%   B.A. Jones 2-7-95, ZP You 9-17-98
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.12 $  $Date: 2000/05/26 18:53:04 $

if nargin <  4,
    error('Requires four input arguments.'); 
end

[errorcode x nu1 nu2 delta] = distchck(4,x,nu1,nu2,delta);

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

[m,n] = size(x);

% Initialize p to zero.
y = zeros(m,n);
p = y;

k = (nu1 <= 0 | nu2 <= 0 | x < 0 | delta < 0);
k1 = ~k;
% set Y to NaN for indices where X is negative or NU1 or NU2 are not positive.
p(k) = NaN;
if ~any(k1), return; end

if any(k) % reset variable so that we don't have to deal with bad indices.
   y = y(k1);
   x = x(k1);
   nu1 = nu1(k1);
   nu2 = nu2(k1);
   delta = delta(k1);
end

% to simply computation, pre-divide nu1,nu2 and delta
nu1 = nu1/2; nu2 = nu2/2; delta = delta/2;


%Value passed to Beta distribution function.
tmp = nu1.*x./(nu2+nu1.*x);
j = 0;
epsq = eps^(1/4);
epsr = sqrt(eps);

% Sum the series.
j0 = floor(delta);
for j=j0:-1:0
   deltay = poiss_pdf(j,delta).*beta_cdf(tmp,j+nu1,nu2);
   y = y + deltay;
   
   % Convergence test:  change must be small and not increasing
   if all(deltay(:)./(y(:)+epsq)<epsr), break; end
end

for j=j0+1:j0+1000    % so that we don't get into an infinite loop
   deltay = poiss_pdf(j,delta).*beta_cdf(tmp,j+nu1,nu2);
   y = y + deltay;
   
   % Convergence test:  change must be small and not increasing
   if all(deltay(:)./(y(:)+epsq)<epsr), break; end
end

if (j == j0+1000) 
  warning('Failed to converge, use the result cautiously.\n');
end

p(k1) = y;

return

Contact us