00001 function f = spm_Gpdf(x,h,l)
00002 % Probability Density Function (PDF) of Gamma distribution
00003 % FORMAT f = spm_Gpdf(g,h,l)
00004 %
00005 % x - Gamma-variate   (Gamma has range [0,Inf) )
00006 % h - Shape parameter (h>0)
00007 % l - Scale parameter (l>0)
00008 % f - PDF of Gamma-distribution with shape & scale parameters h & l
00009 %_______________________________________________________________________
00010 %
00011 % spm_Gpdf implements the Probability Density Function of the Gamma
00012 % distribution.
00013 %
00014 % Definition:
00015 %-----------------------------------------------------------------------
00016 % The PDF of the Gamma distribution with shape parameter h and scale l
00017 % is defined for h>0 & l>0 and for x in [0,Inf) by: (See Evans et al.,
00018 % Ch18, but note that this reference uses the alternative
00019 % parameterisation of the Gamma with scale parameter c=1/l)
00020 %
00021 %           l^h * x^(h-1) exp(-lx)
00022 %    f(x) = ---------------------
00023 %                gamma(h)
00024 %
00025 % Variate relationships: (Evans et al., Ch18 & Ch8)
00026 %-----------------------------------------------------------------------
00027 % For natural (strictly +ve integer) shape h this is an Erlang distribution.
00028 %
00029 % The Standard Gamma distribution has a single parameter, the shape h.
00030 % The scale taken as l=1.
00031 %
00032 % The Chi-squared distribution with v degrees of freedom is equivalent
00033 % to the Gamma distribution with scale parameter 1/2 and shape parameter v/2.
00034 %
00035 % Algorithm:
00036 %-----------------------------------------------------------------------
00037 % Direct computation using logs to avoid roundoff errors.
00038 %
00039 % References:
00040 %-----------------------------------------------------------------------
00041 % Evans M, Hastings N, Peacock B (1993)
00042 %       "Statistical Distributions"
00043 %        2nd Ed. Wiley, New York
00044 %
00045 % Abramowitz M, Stegun IA, (1964)
00046 %       "Handbook of Mathematical Functions"
00047 %        US Government Printing Office
00048 %
00049 % Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
00050 %       "Numerical Recipes in C"
00051 %        Cambridge
00052 %_______________________________________________________________________
00053 % @(#)spm_Gpdf.m        2.2 Andrew Holmes 99/04/26
00054 
00055 %-Format arguments, note & check sizes
00056 %-----------------------------------------------------------------------
00057 if nargin<3, error('Insufficient arguments'), end
00058 
00059 ad = [ndims(x);ndims(h);ndims(l)];
00060 rd = max(ad);
00061 as = [  [size(x),ones(1,rd-ad(1))];...
00062         [size(h),ones(1,rd-ad(2))];...
00063         [size(l),ones(1,rd-ad(3))]     ];
00064 rs = max(as);
00065 xa = prod(as,2)>1;
00066 if sum(xa)>1 & any(any(diff(as(xa,:)),1))
00067         error('non-scalar args must match in size'), end
00068 
00069 %-Computation
00070 %-----------------------------------------------------------------------
00071 %-Initialise result to zeros
00072 f = zeros(rs);
00073 
00074 %-Only defined for strictly positive h & l. Return NaN if undefined.
00075 md = ( ones(size(x))  &  h>0  &  l>0 );
00076 if any(~md(:)), f(~md) = NaN;
00077         warning('Returning NaN for out of range arguments'), end
00078 
00079 %-Degenerate cases at x==0: h<1 => f=Inf; h==1 => f=l; h>1 => f=0
00080 ml = ( md  &  x==0  &  h<1 );
00081 f(ml) = Inf;
00082 ml = ( md  &  x==0  &  h==1 ); if xa(3), mll=ml; else mll=1; end
00083 f(ml) = l(mll);
00084 
00085 %-Compute where defined and x>0
00086 Q  = find( md  &  x>0 );
00087 if isempty(Q), return, end
00088 if xa(1), Qx=Q; else Qx=1; end
00089 if xa(2), Qh=Q; else Qh=1; end
00090 if xa(3), Ql=Q; else Ql=1; end
00091 
00092 %-Compute
00093 f(Q) = exp( (h(Qh)-1).*log(x(Qx)) +h(Qh).*log(l(Ql)) - l(Ql).*x(Qx)...
00094                 -gammaln(h(Qh)) );