Home > wafo > wstats > wggampdf.m

wggampdf

PURPOSE ^

Generalized Gamma probability density function

SYNOPSIS ^

f = wggampdf(x,a,b,c);

DESCRIPTION ^

 WGGAMPDF Generalized Gamma probability density function
 
  CALL:  f = wggampdf(x,a,b,c);
 
         f = density function evaluated at x
     a,b,c = parameters   (default b=1,c=1)
 
  The generalized Gamma distribution is defined by its pdf
 
  f(x;a,b,c)=c*x^(a*c-1)/b^(a*c)*exp(-(x/b)^c)/gamma(a), x>=0, a,b,c>0
  
  Example: 
    x = linspace(0,7,200);
    p1 = wggampdf(x,1,2,1); p2 = wggampdf(x,3,1,1);
    plot(x,p1,x,p2)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function f = wggampdf(x,a,b,c);
002 %WGGAMPDF Generalized Gamma probability density function
003 %
004 % CALL:  f = wggampdf(x,a,b,c);
005 %
006 %        f = density function evaluated at x
007 %    a,b,c = parameters   (default b=1,c=1)
008 %
009 % The generalized Gamma distribution is defined by its pdf
010 %
011 % f(x;a,b,c)=c*x^(a*c-1)/b^(a*c)*exp(-(x/b)^c)/gamma(a), x>=0, a,b,c>0
012 % 
013 % Example: 
014 %   x = linspace(0,7,200);
015 %   p1 = wggampdf(x,1,2,1); p2 = wggampdf(x,3,1,1);
016 %   plot(x,p1,x,p2)
017 
018 % Reference: Cohen & Whittle, (1988) "Parameter Estimation in Reliability
019 % and Life Span Models", p. 220 ff, Marcel Dekker.
020 
021 % Tested on; Matlab 5.3
022 % History:
023 % revised pab 24.10.2000
024 %  - added comnsize, nargchk
025 % added ms 09.08.2000
026 
027 error(nargchk(2,4,nargin))
028 
029 if nargin<3|isempty(b),  b=1; end
030 if nargin<4|isempty(c),  c=1; end
031 
032 [errorcode x a b c] = comnsize(x,a,b,c);
033 
034 if errorcode > 0
035     error('x, a, b and c must be of common size or scalar.');
036 end
037 
038 % Initialize f to zero.
039 f = zeros(size(x));
040 
041 ok= ((a >0) & (b> 0) & (c>0));
042 
043 k=find(x > 0 & ok);
044 if any(k),
045   % old call
046   %  f(k)=c(k).*x(k).^(a(k).*c(k)-1)./b(k).^(a(k).*c(k))....
047   %      .*exp(-(x(k)./b(k)).^c(k))./gamma(a(k));
048   
049   % new call: more stable for large a
050   tmp = (a(k).*c(k) - 1) .* log(x(k)) - (x(k) ./ ...
051       b(k)).^c(k) - gammaln(a(k)) - a(k).*c(k) .* log(b(k));
052   f(k) = c(k).*exp(tmp);
053 end
054  
055 % Special cases for x==0: (this avoids warning messages)
056 k1 = find(x == 0 & a.*c < 1 & ok);
057 if any(k1),
058   tmp = Inf;
059   f(k1) = tmp(ones(size(k1)));
060 end
061 k2 = find(x == 0 & a.*c == 1 & ok);
062 if any(k2)
063   f(k2) = (c(k2)./b(k2)./gamma(a(k2)));
064 end
065 
066 %   Return NaN if the arguments are outside their respective limits.
067 k3 = find(~ok);     
068 if any(k3)
069   tmp = NaN;
070   f(k3) = tmp(ones(size(k3)));
071 end
072 
073 
074 
075

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003