Probability Density Function (PDF) of Gamma distribution
FORMAT f = spm_Gpdf(g,h,l)
x - Gamma-variate (Gamma has range [0,Inf) )
h - Shape parameter (h>0)
l - Scale parameter (l>0)
f - PDF of Gamma-distribution with shape & scale parameters h & l
_______________________________________________________________________
spm_Gpdf implements the Probability Density Function of the Gamma
distribution.
Definition:
-----------------------------------------------------------------------
The PDF of the Gamma distribution with shape parameter h and scale l
is defined for h>0 & l>0 and for x in [0,Inf) by: (See Evans et al.,
Ch18, but note that this reference uses the alternative
parameterisation of the Gamma with scale parameter c=1/l)
l^h * x^(h-1) exp(-lx)
f(x) = ---------------------
gamma(h)
Variate relationships: (Evans et al., Ch18 & Ch8)
-----------------------------------------------------------------------
For natural (strictly +ve integer) shape h this is an Erlang distribution.
The Standard Gamma distribution has a single parameter, the shape h.
The scale taken as l=1.
The Chi-squared distribution with v degrees of freedom is equivalent
to the Gamma distribution with scale parameter 1/2 and shape parameter v/2.
Algorithm:
-----------------------------------------------------------------------
Direct computation using logs to avoid roundoff errors.
References:
-----------------------------------------------------------------------
Evans M, Hastings N, Peacock B (1993)
"Statistical Distributions"
2nd Ed. Wiley, New York
Abramowitz M, Stegun IA, (1964)
"Handbook of Mathematical Functions"
US Government Printing Office
Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
"Numerical Recipes in C"
Cambridge
_______________________________________________________________________
@(#)spm_Gpdf.m 2.2 Andrew Holmes 99/04/26
0001 function f = spm_Gpdf(x,h,l) 0002 % Probability Density Function (PDF) of Gamma distribution 0003 % FORMAT f = spm_Gpdf(g,h,l) 0004 % 0005 % x - Gamma-variate (Gamma has range [0,Inf) ) 0006 % h - Shape parameter (h>0) 0007 % l - Scale parameter (l>0) 0008 % f - PDF of Gamma-distribution with shape & scale parameters h & l 0009 %_______________________________________________________________________ 0010 % 0011 % spm_Gpdf implements the Probability Density Function of the Gamma 0012 % distribution. 0013 % 0014 % Definition: 0015 %----------------------------------------------------------------------- 0016 % The PDF of the Gamma distribution with shape parameter h and scale l 0017 % is defined for h>0 & l>0 and for x in [0,Inf) by: (See Evans et al., 0018 % Ch18, but note that this reference uses the alternative 0019 % parameterisation of the Gamma with scale parameter c=1/l) 0020 % 0021 % l^h * x^(h-1) exp(-lx) 0022 % f(x) = --------------------- 0023 % gamma(h) 0024 % 0025 % Variate relationships: (Evans et al., Ch18 & Ch8) 0026 %----------------------------------------------------------------------- 0027 % For natural (strictly +ve integer) shape h this is an Erlang distribution. 0028 % 0029 % The Standard Gamma distribution has a single parameter, the shape h. 0030 % The scale taken as l=1. 0031 % 0032 % The Chi-squared distribution with v degrees of freedom is equivalent 0033 % to the Gamma distribution with scale parameter 1/2 and shape parameter v/2. 0034 % 0035 % Algorithm: 0036 %----------------------------------------------------------------------- 0037 % Direct computation using logs to avoid roundoff errors. 0038 % 0039 % References: 0040 %----------------------------------------------------------------------- 0041 % Evans M, Hastings N, Peacock B (1993) 0042 % "Statistical Distributions" 0043 % 2nd Ed. Wiley, New York 0044 % 0045 % Abramowitz M, Stegun IA, (1964) 0046 % "Handbook of Mathematical Functions" 0047 % US Government Printing Office 0048 % 0049 % Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992) 0050 % "Numerical Recipes in C" 0051 % Cambridge 0052 %_______________________________________________________________________ 0053 % @(#)spm_Gpdf.m 2.2 Andrew Holmes 99/04/26 0054 0055 %-Format arguments, note & check sizes 0056 %----------------------------------------------------------------------- 0057 if nargin<3, error('Insufficient arguments'), end 0058 0059 ad = [ndims(x);ndims(h);ndims(l)]; 0060 rd = max(ad); 0061 as = [ [size(x),ones(1,rd-ad(1))];... 0062 [size(h),ones(1,rd-ad(2))];... 0063 [size(l),ones(1,rd-ad(3))] ]; 0064 rs = max(as); 0065 xa = prod(as,2)>1; 0066 if sum(xa)>1 & any(any(diff(as(xa,:)),1)) 0067 error('non-scalar args must match in size'), end 0068 0069 %-Computation 0070 %----------------------------------------------------------------------- 0071 %-Initialise result to zeros 0072 f = zeros(rs); 0073 0074 %-Only defined for strictly positive h & l. Return NaN if undefined. 0075 md = ( ones(size(x)) & h>0 & l>0 ); 0076 if any(~md(:)), f(~md) = NaN; 0077 warning('Returning NaN for out of range arguments'), end 0078 0079 %-Degenerate cases at x==0: h<1 => f=Inf; h==1 => f=l; h>1 => f=0 0080 ml = ( md & x==0 & h<1 ); 0081 f(ml) = Inf; 0082 ml = ( md & x==0 & h==1 ); if xa(3), mll=ml; else mll=1; end 0083 f(ml) = l(mll); 0084 0085 %-Compute where defined and x>0 0086 Q = find( md & x>0 ); 0087 if isempty(Q), return, end 0088 if xa(1), Qx=Q; else Qx=1; end 0089 if xa(2), Qh=Q; else Qh=1; end 0090 if xa(3), Ql=Q; else Ql=1; end 0091 0092 %-Compute 0093 f(Q) = exp( (h(Qh)-1).*log(x(Qx)) +h(Qh).*log(l(Ql)) - l(Ql).*x(Qx)... 0094 -gammaln(h(Qh)) );