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

pfq(a,b,z,d)
function [y, tt, nterms] = pfq(a,b,z,d)
%PFQ Generalized Hypergeometric Function pFq
%   [Y,TT,NTERMS] = PFQ(A,B,Z,D) returns the Generalized Hypergeometric
%   Function pFq with numerator parameters A, denominator parameters B,
%   evaluated at the values in Z, with respect to derivatives D.
%
%   PFQ(A,B,Z) returns the values evaluated at Z
%
%   Inputs:
%   A : real or complex valued vector (or empty [])
%   B : real or complex valued vector (or empty [])
%   Z : real or complex valued M-by-N-by-... array
%   D : scalar or vector of real non-negative integers
%
%   Outputs:
%   Y      : Evaluation of pFq at values in Z, with respect to derivative D
%   TT     : Estimated Truth Table on values given
%              +1 : Terms were convergent, answer valid close to precision
%               0 : Uncharacterized, use answers with caution
%              -1 : Terms were divergent, use answers with severe caution
%   NTERMS : Number of terms used in the calculation for each output
%
%   If Z is M-by-N-by-... array and D is a vector 1-by-Q then the outputs
%   will be an array of size M-by-N-by-...-by-Q
%
%   Examples:
%   %0F1: Confluent Hypergeometric Function
%   pfq( [] , 1 , 1.5 )
%   pfq( [] , 1+1i , 1.5-1i )
%   pfq( [] , 1 , -2 )
%   %1F1: Kummer Confluent Hypergeometric Function
%   pfq( 1 , 2 , 3 )
%   pfq( 2+1i , 2, 0.5 )
%   pfq( 10 , 1/3 , -1 )
%   %2F1: Hypergeometric Function
%   pfq( [1/3 1/3] , 2/3 , 1/2 )
%   pfq( [2+1i -1i] , 3/4 , 0.5-0.5i )
%   %PFQ: Generalized Hypergeometric Function
%   pfq( [1 1] , [3 3 3] , 2 )
%   pfq( [1i 1i 1i] , [2 2 2] , -1i )
%   pfq( [1 2 3 4] , [5 6 7] , [.1 .3 .5] )
%
%   Vectorization with derivatives:
%   Evaluate pFq above, at [.1 .3 .5], but with respect to the 0th, 1st and
%   2nd derivative. Show all outputs
%   [y , tt , nterms] = pfq( [1 2 3 4] , [5 6 7] , [.1 .3 .5] , 0:2 )
%
%   Finite degree polynomial, with derivatives
%   x = [.1 .3 .5];
%   y = pfq( [1 -2] , -2 , x , [0:4] )
%   [1+x+x.^2; 1+2.*x] %Theoretical answer
%

%
%   Mike Sheppard
%   MIT Lincoln Laboratory
%   michael.sheppard@ll.mit.edu
%   Last Modified 25-Jul-2011

%
%   Comments, suggestions, or improvements are welcome
%


%Check inputs
if nargin < 3
    error('pfq:TooFewInputs',...
        'Requires at least three input argument.');
end
if nargin == 3
    d=0; %Default zero'th derivative
end
if ~isvector(a) && ~isempty(a)
    error('pfq:A must be either a vector or empty');
end
if ~isvector(b) && ~isempty(b)
    error('pfq:B must be either a vector or empty');
end
if ~(isvector(d) && all(d==round(d)) && all(d>=0))
    error('pfq:D must be either a scalar or vector of real non-negative integers');
end


%Tolerance
tol_exp=50;
%Tolerance Exponent, base 10
%   Terms are added until the norm of each term, for all inputs and
%   derivatives, are ...
%   1. Less than 10^(-TOL_EXP) in absolute value; converging
%   2. Greater than 10^(TOL_EXP) in absolute value; diverging
%   Note: It may be oscillatory and an exact answer exists, however this
%   will stop when the terms exceed 10^(TOL_EXP) in absolute value and
%   yields -1 for Estimated Truth Table, above.

%Maximum number of steps
maxsteps=1e4;
    



%
%   The definition of pFq by a formal power series is:
%   pFq = sum( [C(a,n) / C(b,n)] * [z^n/n!] , n=0, ... , infinity )
%      with
%   C(v,n)=prod( [gamma(v(j)+n) / gamma(v(j)] , j=1, ... ,length(v) )
%
%   The derivatives are calculated equivalently, with [C(a,n)/C(b,n)]
%   term remaining the same and using d'th derivative of
%   z^k/k! -> z^(k-d) / (k-d)!
%
%   Terms are added until all the terms, for all inputs and
%   derivatives, exceed TOL given below; or if expansion terminates
%   with a finite polynomial or maximum allowed number of steps
%


% INITIALIZE PARAMETERS
%Make A and B row vectors for consistency
a=(a(:)).'; b=(b(:)).'; %non-conjugate transpose of column vector
%Keep track of which A's and B's are negative integers
anegint=(a<0)&(a==round(a)); %Negative integer values of A
bnegint=(b<0)&(b==round(b)); %Negative integer values of B
%If any A's is a negative integer a polynomial is produce,
%find stopping point
if any(anegint)
    minnega=min(-a(anegint)); %stopping point
else
    minnega=Inf; %no stopping
end

sz=size(z); ld=length(d); z=z(:); lnz=log(z); lz=length(lnz);
y=zeros([lz ld]); tol_high=10^abs(tol_exp); tol_low=tol_high^-1;
n=0; keepgoing=1;
indxz=true(lz,ld); %Logical index of which Z, and which derivatives, to still work on
tt=zeros(lz,ld); %Estimated Truth Table
nterms=ones(lz,ld); %Number of terms, minimum number is one


while keepgoing
    
    %Compute the log of product of terms
    %C(v,n)=prod( [gamma(v(j)+n) / gamma(v(j)] , j=1, ... ,length(v) )
    %Returns log of product, by summing the individual log of terms
    %Subfunction watches out for negative integers and does those
    %separately
    lnnum=lnprodterms(a,anegint,n);
    lnden=lnprodterms(b,bnegint,n);
    term=zeros([lz ld]); %set to zero
    
    %Loop through each derivative given
    for dn=1:ld
        %Use only Z values that need more computation, for this derivative
        indxzd=indxz(:,dn); z_temp=z(indxzd); findxzd=find(indxzd);
        if ~isempty(findxzd)
            %For each derivative, check to see if term in original expansion is
            %valid part of derivative. The indexing is for the actual function;
            %and does not change for derivatives, but coefficient and exponent
            %do change accordingly
            %d'th derivative of (coef) * z^n / n! -> (coef) * z^(n-d) / (n-d)!
            if n>=d(dn)
                %Rewrite in log space, using d'th derivative
                lnzn=(n-d(dn)).*lnz(indxzd); lnfact=lngamma(n-d(dn)+1);
                %Natural log of term for all z for specific derivative
                lnterm=lnnum-lnden-lnfact+lnzn;
                %Check special case if z==0. If so, and if n==d(dn) then term
                %is coef as computed, else zero. If z~=0 then coef is correct.
                ztz=(z_temp==0); indx0=findxzd(ztz); indx1=findxzd(~ztz);
                %keyboard
                term(indx0,dn)=(n==d(dn)).*exp(lnnum-lnden-lnfact);
                term(indx1,dn)=exp(lnterm(~ztz));
                %Record number of steps
                %Term of original series - #derivative + one
                %as started with zero'th term (z^0)
                nterms(findxzd,dn)=n-d(dn)+1;
                %Terminate if all values are outside range of tolerance;
                %either diverging or converging, keep all those within range to
                %keep on going with terms
                normM=abs(term(indxzd,dn));
                lo_tol=(normM<tol_low); hi_tol=(normM(:,1)>tol_high);
                if any(lo_tol | hi_tol)
                    %Add either +1 or -1 to Estimated Truth Table
                    tt(findxzd(lo_tol),dn)=1; tt(findxzd(hi_tol),dn)=-1;
                    %Delete those points from further calculations
                    indxz(findxzd(lo_tol | hi_tol),dn)=false;
                end
            else
                %If term in original expansion is not a valid part of the
                %derivative then just add zero
                term(:,dn)=0;
            end
        end  %~isempty(findxzd)
    end
    
    %Add the new term
    y=y+term;
    
    %Break out of loop if n is equal to minimum of any A's that are
    %negative integers (finite degree polynomial)
    if n==minnega
        keepgoing=0;
    end
    
    %If all new terms for all values exceed the tolerance; stop.
    if all(indxz(:)==false)
        keepgoing=0; %Stop, all values exceed tolerance
    else
        n=n+1; %Keep going to the next term
    end
    
    %Break out of loop after MAXSTEPS terms, if haven't already
    %This is usually the case when z is close to 1
    if n>maxsteps
        keepgoing=0;
        %Record all remaining ones as uncharacterized, 
        %and maxsteps number of steps. 
        nterms(indxz)=maxsteps;
        tt(indxz)=0; %uncharacterized
    end
    
    
end


%For values of z==0, results are accurate for any derivative
if any(z==0)
    k0=(z==0); tt(k0,:)=1; 
end

%For finite polynomials the results are accurate for any derivative
%except if output exceeded given tolerance
if isfinite(minnega)
    tt(tt~=-1)=1; %keep -1 the same, change all others to 1
end

%For those with MAXSTEPS, change TT to 0 for uncharacterized
tt(nterms==maxsteps)=0;

%Reshape final result to proper dimensions from input
y=reshape(y(:),[sz ld]);
tt=reshape(tt(:),[sz ld]);
nterms=reshape(nterms(:),[sz ld]);
if isvector(z)
    if sz(1)==1 && ld>1
        y=(squeeze(y))';
        tt=(squeeze(tt))';
        nterms=(squeeze(nterms))';
    else
        y=squeeze(y);
        tt=squeeze(tt);
        nterms=squeeze(nterms);
    end
end


%Fix round-off
lim=100*eps;
%Make real part zero, keep imaginary
kr=abs(real(y))<lim; if any(kr), y(kr)=imag(y(kr)); end
%Make imag part zero, keep real part
ki=abs(imag(y))<lim; if any(ki), y(ki)=real(y(ki)); end


return




function lnprod=lnprodterms(v,vnegint,n)
%Given vector A, and term n compute
%Log[ (v_1)(n) ... (v_t)(n) ] where
%(v)(n)= (v)*(v+1)*(v+2)*...*(v+n-1) [Pochhammer symbol]
%Where vnegint is logical index of negative integers within V

%Step 1: Compute non-negative integers (valid input for lngamma)
lnprod=sum(lngamma(v(~vnegint)+n)-lngamma(v(~vnegint)));

%Step 2: Add the remaining terms that use negative integers
if (any(vnegint))&&(n>0)  %Exclude n==0 term
    A=v(vnegint); A=A(:).'; B=0:n-1; %Correct orientation for bsxfun
    sumv=bsxfun(@plus,A,B);
    %Trivial Example: A=[4 7 16].', B=0:3, sumv=[A A A A]+[B;B;B]
    lnprod=lnprod+sum(log(sumv(:))); %sum of logs -> log of products
end

return






%THE CODE BELOW IS FROM
%http://www.mathworks.com/matlabcentral/fileexchange/ ...
%  ... 978-special-functions-math-library/content/gammaln.m
%
%Modifications:
%   1. Renamed to lngamma so as not to override built-in gamma or gammaln
%   2. Replaced i with 1i
% -Mike Sheppard (7/20/2011)
%

function [f] = lngamma(z)
% GAMMALOG  Natural Log of the Gamma function valid in the entire complex plane.
%           This routine uses an excellent Lanczos series approximation
%           for the complex ln(Gamma) function.
%
%usage: [f] = gammaln(z)
%             z may be complex and of any size.
%             Also  n! = prod(1:n) = exp(gammalog(n+1))
%
%tested under version 5.3.1
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%
%see also:   GAMMA GAMMALN GAMMAINC PSI
%see also:   mhelp GAMMA
%see also:   mhelp lnGAMMA

%Paul Godfrey
%pgodfrey@conexant.com
%07-13-01


siz = size(z);
z=z(:);
zz=z;

f = 0.*z; % reserve space in advance

p=find(real(z)<0);
if ~isempty(p)
    z(p)=-z(p);
end

%Lanczos approximation for the complex plane

g=607/128; % best results when 4<=g<=5

c = [  0.99999999999999709182;
    57.156235665862923517;
    -59.597960355475491248;
    14.136097974741747174;
    -0.49191381609762019978;
    .33994649984811888699e-4;
    .46523628927048575665e-4;
    -.98374475304879564677e-4;
    .15808870322491248884e-3;
    -.21026444172410488319e-3;
    .21743961811521264320e-3;
    -.16431810653676389022e-3;
    .84418223983852743293e-4;
    -.26190838401581408670e-4;
    .36899182659531622704e-5];

s=0;
for k=size(c,1):-1:2
    s=s+c(k)./(z+(k-2));
end

zg=z+g-0.5;
s2pi= 0.9189385332046727417803297;

f=(s2pi + log(c(1)+s)) - zg + (z-0.5).*log(zg);

f(z==1 | z==2) = 0.0;

if ~isempty(p)
    lpi= 1.14472988584940017414342735 + 1i*pi;
    f(p)=lpi-log(zz(p))-f(p)-log(sin(pi*zz(p)));
end

p=find(round(zz)==zz & imag(zz)==0 & real(zz)<=0);
if ~isempty(p)
    f(p)=Inf;
end

f=reshape(f,siz);

return

Contact us