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:
 comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. gamma Gamma function. gammaln Logarithm of gamma function. inf Infinity. nan Not-a-Number.
This function is called by:
 dist2dfun is an internal function to dist2dcdf dist2dprb. dist2dpdf Joint 2D PDF computed as f(x1|X2=x2)*f(x2) jhsnlpdf Joint (Scf,Hd) PDF for nonlinear waves with a JONSWAP spectra. jhspdf Joint (Scf,Hd) PDF for linear waves with JONSWAP spectra. ltwcpdf Long Term Wave Climate PDF of significant wave height and wave period ohhpdf Marginal wave height, Hd, PDF for Bimodal Ochi-Hubble spectra. ohhspdf Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. ohhsspdf Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra. ohhvpdf Joint (Vcf,Hd) PDF for lineare waves with Ochi-Hubble spectra. thpdf Marginal wave height, Hd, PDF for Torsethaugen spectra. thsspdf Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. wchi2pdf Chi squared probability density function wgampdf Gamma probability density function

## 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
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