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

connegbin_pdf(n, k, P)
% connegbin_pdf.m - evaluates a Condensed Negative Binomial Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%   Wiley, p.247, 2005.  
%
%   Created by  J. Huntley,  12/07/06.
%

function [pdf] = connegbin_pdf(n, k, P)

%persistent logk logP logopp

%if(isempty(logk))
    logk = gammaln(k);
    logP = log(P);
    logopp = log(1+P);
%end

% Initializations.
twox = 2*n;

% Calculate PDF.
if(n == 0)
    a0 = 1 / (1+P)^k;
    a1 = k * P / (1+P)^(k+1);
    pdf = a0 + 0.5 * a1;
elseif(n > 0)
    %a2xm1 = factorial(k+twox-2) * P^(twox-1) / (factorial(twox-1) * factorial(k-1) ...
    %        * (1 + P)^(k+twox-1));
    %a2x = factorial(k+twox-1) * P^twox / (factorial(twox) * factorial(k-1) ...
    %        * (1 + P)^(k+twox));
    %a2xp1 = factorial(k+twox) * P^(twox+1) / (factorial(twox+1) * factorial(k-1) ...
    %        * (1 + P)^(k+twox+1));
    a2xm1 = exp(gammaln(k+twox-1) + (twox-1)*logP - (gammaln(twox) + logk ...
            + (k+twox-1)*logopp));
    a2x = exp(gammaln(k+twox) + twox*logP - (gammaln(twox+1) + logk ...
            + (k+twox)*logopp));
    a2xp1 = exp(gammaln(k+twox+1) + (twox+1)*logP - (gammaln(twox+2) + logk ...
            + (k+twox+1)*logopp));
    pdf = 0.5 * a2xm1 + a2x + 0.5 * a2xp1;
end

return


    

Contact us