Probability Density Function (PDF) of Gamma distribution
FORMAT f = pr_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
_______________________________________________________________________
Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0001 function f = pr_spm_gpdf(x,h,l) 0002 % Probability Density Function (PDF) of Gamma distribution 0003 % FORMAT f = pr_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 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience 0054 0055 % Andrew Holmes 0056 % $Id: spm_Gpdf.m 112 2005-05-04 18:20:52Z john $ 0057 0058 0059 %-Format arguments, note & check sizes 0060 %----------------------------------------------------------------------- 0061 if nargin<3, error('Insufficient arguments'), end 0062 0063 ad = [ndims(x);ndims(h);ndims(l)]; 0064 rd = max(ad); 0065 as = [ [size(x),ones(1,rd-ad(1))];... 0066 [size(h),ones(1,rd-ad(2))];... 0067 [size(l),ones(1,rd-ad(3))] ]; 0068 rs = max(as); 0069 xa = prod(as,2)>1; 0070 if sum(xa)>1 & any(any(diff(as(xa,:)),1)) 0071 error('non-scalar args must match in size'), end 0072 0073 %-Computation 0074 %----------------------------------------------------------------------- 0075 %-Initialise result to zeros 0076 f = zeros(rs); 0077 0078 %-Only defined for strictly positive h & l. Return NaN if undefined. 0079 md = ( ones(size(x)) & h>0 & l>0 ); 0080 if any(~md(:)), f(~md) = NaN; 0081 warning('Returning NaN for out of range arguments'), end 0082 0083 %-Degenerate cases at x==0: h<1 => f=Inf; h==1 => f=l; h>1 => f=0 0084 ml = ( md & x==0 & h<1 ); 0085 f(ml) = Inf; 0086 ml = ( md & x==0 & h==1 ); if xa(3), mll=ml; else mll=1; end 0087 f(ml) = l(mll); 0088 0089 %-Compute where defined and x>0 0090 Q = find( md & x>0 ); 0091 if isempty(Q), return, end 0092 if xa(1), Qx=Q; else Qx=1; end 0093 if xa(2), Qh=Q; else Qh=1; end 0094 if xa(3), Ql=Q; else Ql=1; end 0095 0096 %-Compute 0097 f(Q) = exp( (h(Qh)-1).*log(x(Qx)) +h(Qh).*log(l(Ql)) - l(Ql).*x(Qx)... 0098 -gammaln(h(Qh)) );