001 function [phat, cov,pci]=wggamfit(data1, plotflag,chat0);
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
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(:));
045 chatDeterministic = 0;
046
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));
053
054 if 1,
055
056
057 x = [linspace(max(eps,start/100),start,25) ...
058 linspace(start+1,10*start+20,15)];
059
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
072 sl = sign(ll);
073
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
082 mvrs = version;
083 ix = find(mvrs=='.');
084 mvrs = str2num(mvrs(1:ix(2)-1));
085 def=1;
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;
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,
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
140 end
141 if mvrs>5.2
142 phat = fminsearch('loglike',phat0,[],data,'wggampdf');
143
144
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
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
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