001 function l=wggambfit(b,data,F,def)
002
003
004
005
006
007
008 if nargin<4|isempty(def),def=1;end
009
010 monitor = logical(0);
011 switch def
012 case 1,
013
014 ld = F;
015 mld = mean(ld);
016 db = data.^b;
017 sdb = sum(db);
018 a = -(b*(mld-sum(db.*ld)/sdb))^(-1);
019 h = 1e-6;
020 h1 = .5*1e+6;
021 if ((a<=0) | isnan(a)),
022 l = NaN;
023 elseif (a<=h)
024 l = log(length(data)*a) - (gammaln(a+h)-gammaln(a))/h + b*mld ....
025 -log(sdb);
026 else
027 l = log(length(data)*a) - (gammaln(a+h)-gammaln(a-h))/(2*h) + b*mld ....
028 -log(sdb);
029 end
030 case 2,
031 x = data;
032 tmp = sqrt(-log(1-F));
033 tmp2 = sqrt(-log(1-wggamcdf(x,b(1),b(2),b(3))));
034 if monitor
035 plot(x,[ tmp tmp2]); drawnow
036 end
037 l = mean(abs(tmp-tmp2).^(2));
038 case 3,
039 l = ...
040 sum([(gamma(b(1))*gamma(b(1)+2/b(2))/gamma(b(1)+1/b(2))^2-data)^2+...
041 (sqrt(gamma(b(1)))*gamma(b(1)+3/b(2))/gamma(b(1)+2/b(2))^1.5-F)^2]);
042 case 4,
043 l = sum([(sqrt(gamma(b(1)))*gamma(b(1)+3/b(2))/gamma(b(1)+2/b(2))^1.5-data)^2+...
044 (gamma(b(1))*gamma(b(1)+4/b(2))/gamma(b(1)+2/b(2))^2-F)^2]);
045 end
046
047 if monitor
048 disp(['err = ' num2str(l,10) ' phat = ' num2str(b,4) ])
049 end
050