Home > wafo > wstats > wggamfit.m

# wggamfit

## PURPOSE

Parameter estimates for Generalized Gamma data.

## SYNOPSIS

[phat, cov,pci]=wggamfit(data1, plotflag,chat0);

## DESCRIPTION

``` WGGAMFIT Parameter estimates for Generalized Gamma data.

CALL:  [phat, cov] = wggamfit(data, plotflag)

phat = [a,b,c] = the maximum likelihood estimates of the
parameters of the Generalized Gamma distribution
(see wggampdf) given the data.
cov  = estimated asymptotic covariance matrix of estimates
data = data vector
plotflag = 0, do not plot
> 0, plot the empiricial distribution function and the
estimated cdf (see empdistr for options)(default)

Example:
R = wggamrnd(2,2,2,1,100);
[phat, cov] = wggamfit(R,2)

## CROSS-REFERENCE INFORMATION

This function calls:
 empdistr Computes and plots the empirical CDF wgamfit Parameter estimates for Gamma data. wggambfit Is an internal routine for wggamfit wggamcdf Generalized Gamma cumulative distribution function wnorminv Inverse of the Normal distribution function deblank Remove trailing blanks. diff Difference and approximate derivative. error Display message and abort function. fmins fminsearch Multidimensional unconstrained nonlinear minimization (Nelder-Mead). fzero Scalar nonlinear zero finding. gammaln Logarithm of gamma function. linspace Linearly spaced vector. mean Average or mean value. nan Not-a-Number. optimset Create/alter OPTIM OPTIONS structure. std Standard deviation. str2num Convert string matrix to numeric array. title Graph title. version MATLAB version number.
This function is called by:
 dist2dfit Parameter estimates for DIST2D data.

## SOURCE CODE

```001 function [phat, cov,pci]=wggamfit(data1, plotflag,chat0);
002 %WGGAMFIT Parameter estimates for Generalized Gamma data.
003 %
004 % CALL:  [phat, cov] = wggamfit(data, plotflag)
005 %
006 %     phat = [a,b,c] = the maximum likelihood estimates of the
007 %            parameters of the Generalized Gamma distribution
008 %            (see wggampdf) given the data.
009 %     cov  = estimated asymptotic covariance matrix of estimates
010 %     data = data vector
011 % plotflag = 0, do not plot
012 %          > 0, plot the empiricial distribution function and the
013 %               estimated cdf (see empdistr for options)(default)
014 %
015 % Example:
016 %   R = wggamrnd(2,2,2,1,100);
017 %   [phat, cov] = wggamfit(R,2)
018 %
020
021 % Reference: Cohen & Whittle, (1988) "Parameter Estimation in Reliability
022 % and Life Span Models", p. 220 ff, Marcel Dekker.
023
024 % Tested on: Matlab 5.3
025 % History:
026 % revised pab 23.01.2004
028 % revised pab 27.01.2001
029 %  - Improved the stability of estimation by adding a linear search to
030 %    find the largest sign change with positive derivative.
031 % revised pab 13.11.2000
032 %  - added check on estimated parameters -> return NaN when any is less
033 %    than zero
034 % revised pab
035 %  -changed the order of a,b,c
036 %  -safer call to fzero
037 % rewritten ms 10.08.2000
038
039
040 error(nargchk(1,3,nargin))
041 if nargin < 2 | isempty(plotflag), plotflag=1;end
042
043 cov   = []; pci=[];
044 data  = sort(data1(:)); % make sure it is a vector
045 chatDeterministic = 0;
046 % Use Weibull start for c (assuming a=1)
047 ld   = log(data);
048 if nargin>2 & ~isempty(chat0) & isfinite(chat0)
049   chat = chat0;
050   chatDeterministic = 1;
051 else
052   start = 1./(6^(1/2)/pi*std(ld)); % approx sqrt(a)*c
053
054   if 1,
055     % Do linear search to find a sign change pab 28.01.2001
056     % -> more stable estimation
057     x  = [linspace(max(eps,start/100),start,25) ...
058       linspace(start+1,10*start+20,15)];
059     %x  = [linspace(eps,start,50) linspace(start+1,start+20,10)];
060     ll = zeros(size(x));
061     ll(1) = wggambfit(x(1),data,ld);
062     for ix=2:length(x),
063       ll(ix) = wggambfit(x(ix),data,ld);
064       if ll(ix-1)<ll(ix)&ll(ix)>0 ,x(ix+1:end)=[];ll(ix+1:end)=[];break,end
065     end
066
067     ind = find(isnan(ll));
068     ll(ind) = [];
069     x(ind)  = [];
070     dll=diff(ll)./diff(x);
071      %figure(1),  plot(x,ll,(x(1:end-1)+x(2:end))/2 ,dll,'r');axis([0 inf -0.1 0.1]),figure(2)
072     sl = sign(ll);
073     %choose the largest sign change with positive derivative
074     ind = max(find(sl(1:end-1)<sl(2:end) ) );
075     if ~isempty(ind)
076       start = x(ind:ind+1);
077     elseif any(dll>0)
078       start = x(min(find(dll>0)));
079     end
080   end
081   %start
082   mvrs = version;
083   ix   = find(mvrs=='.');
084   mvrs = str2num(mvrs(1:ix(2)-1));
085   def=1; % What is def? See 'wggambfun'
086   if (def==1)
087     if (mvrs > 5.2),
088       chat = fzero('wggambfit',start,optimset,data,ld);
089     else
090       chat = fzero('wggambfit',start,sqrt(eps),[],data,ld);
091     end
092   else
093     N = length(data);
094     F = [0.5./N:1./N:(N - 0.5)./N]';
095     size(F)
096     size(data)
097     chat = mean(start)
098     dc   = data.^chat;
099     ahat = -(chat*(mean(ld)-sum(dc.*ld)/sum(dc)))^(-1);
100     bhat = (mean(dc)/ahat).^(1/chat);
101     start = [ahat,bhat,chat]
102     if  (mvrs  >  5.2)
103       phat = fminsearch('wggambfit',start,optimset,data,F,def);
104     else
105       phat = fmins('wggambfit',start,sqrt(eps),[],data,F,def);
106     end
107     chat = phat(3);
108   end
109 end
110 if 0,
111   phatAB = wgamfit(data.^chat,0);
112   ahat = phatAB(1);
113   bhat = phatAB(2).^(1/chat);
114 else
115   dc   = data.^chat;
116   ahat = -(chat*(mean(ld)-sum(dc.*ld)/sum(dc)))^(-1);
117   bhat = (mean(dc)/ahat).^(1/chat);
118 end
119 phat0 = [ahat, bhat, chat];
120 if 0,
121   phat = fminsearch('loglike',phat0,[],data,'wggampdf');
122 else
123   phat = phat0;
124 end
125 if any(phat<=0|isnan(phat)),
126   phat(1:3)=NaN;cov=phat;pci=cov;
127   return,
128   chat = start(end)+10; %1./(6^(1/2)/pi*std(log(data)));
129   ahat =-(chat*(mean(ld)-sum(data.^chat.*ld)/sum(data.^chat)))^(-1);
130   bhat = (mean(data.^chat)/ahat).^(1/chat);
131   phat0 = max([ahat, bhat, chat],sqrt(eps))
132   if 1, % Moments fit
133     Ex = mean(data);
134     Ex2 = mean(data.^2);
135     Ex3 = mean(data.^3);
136     Ex4 = mean(data.^4);
137     phat1 = fmins('wggambfit',phat0(1:2:3),[],[],Ex2/Ex^2,Ex3/Ex2^1.5,3)
138     phat1 = fmins('wggambfit',phat0(1:2:3),[],[],Ex3/Ex2^1.5,Ex4/Ex2^2,4)
139     %phat0  = [mean(data.^start).^(1./start) start];
140   end
141   if mvrs>5.2
142     phat = fminsearch('loglike',phat0,[],data,'wggampdf');
143    % N =length(data);
144    % phat = fminsearch('wggambfit',phat0,[],data,(.5:N)'/N,2)
145   else
146     phat = fmins('loglike',phat0,[],[],data,'wggampdf');
147   end
148   ahat = phat(1);
149   bhat = phat(2);
150   chat = phat(3);
151   if any(phat<=0|isnan(phat)),
152     phat(1:3)=NaN;cov=phat;pci=cov;return,
153   end
154 end
155
156 if nargout>1,
157   h=10^(-5);
158
159   c1=(gammaln(ahat+h)-gammaln(ahat))/h;
160   c2=(gammaln(ahat+1+h)-gammaln(ahat+1))/h;
161   c3=(gammaln(ahat+2*h)-2*gammaln(ahat+h)+gammaln(ahat))/h^2;
162   c4=(gammaln(ahat+1+2*h)-2*gammaln(ahat+1+h)+gammaln(ahat+1))/h^2;
163
164   cov=[c3, -c1/chat, chat/bhat;
165     -c1/chat, (1+ahat*c4+ahat*c2^2)/chat^2, -(1+ahat*c1)/bhat;
166     chat/bhat, -(1+ahat*c1)/bhat, ahat*chat^2/bhat^2]/length(data);
167
168   if (chatDeterministic)
169     cov(3,:) = 0;
170     cov(:,3) = 0;
171   end
172   %[LL, cov] = loglike(phat,data,'wggampdf');
173 end
174 if nargout>2
175   alpha2 = ones(1,3)*0.05/2;
176   var = diag(cov).';
177   pci = wnorminv([alpha2;1-alpha2], [phat;phat],[var;var]);
178 end
179 if plotflag
180 %  sd=sort(data);
181   sd = data;
182   empdistr(sd,[sd, wggamcdf(sd,ahat,bhat,chat)],plotflag)
183   title([deblank(['Empirical and Generalized Gamma estimated cdf'])])
184 end
185
186
187
188
189
190
191```

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