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

geninvtri_pdf(n, nn, p1, p2, p3, p4)
% geninvtri_pdf.m - evaluates a Generalized Inverse Trinomial Probability Density.
%   See "A Generalization of the Inverse Trinomial", K. Aoyama & K. Shimizu,
%   Keio U., Research Report KSTS/RR-05/003, 06/23/05,
%   iscm.math.um.edu.my/abstracts%5CKunio%20Shimizu.pdf.
%
%       For a proper distribution, p1 + p2 >= p4 + p5 !!!
%
%       Vector forom of PDF!!!
%
%  Created by Jim Huntley,  10/22/08
%

function[pdf] = geninvtri_pdf(n, nn, p1, p2, p3, p4)

%persistent p5 llim lognn logp1 logp2 logp3 logp4 logp5

%if(isempty(p5))
    p5 = 1 - p1 - p2 - p3 - p4;
    llim = 20;
    lognn = log(nn);
    logp1 = log(p1);
    logp2 = log(p2);
    logp3 = log(p3);
    logp4 = log(p4);
    logp5 = log(p5);
%end

sumk = 0;
for jk = 1:n+1
   k = jk - 1;
   suml = 0;
   for jl = 1:llim
       l = jl - 1;
       sumi = 0;
       m = min(nn+k+l,n-k);
       for ji = 1:m+1
           i = ji - 1;
           sumi = sumi + exp(lognn + gammaln(nn+n-i+k+2*l+1) + ...
                  (nn-i+k+l)*logp1 + i*logp2 + (n-i-k)*logp3 - ...
                  (gammaln(nn-i+k+l+1)+gammaln(i+1)+gammaln(n-i-k+1)) - ...
                  log(nn+n-i+k+2*l));
       end
       suml = suml + exp(l*logp5 + log(sumi) - gammaln(l+1));
   end
   sumk = sumk + exp(k*logp4 + log(suml) - gammaln(k+1));
end
pdf = sumk;

return



Contact us