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) 
  
  See also   wggamcdf, empdistr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 % 
019 % See also   wggamcdf, empdistr 
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 
027 %  - added chat0  
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