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

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