Home > wafo > wstats > wggambfit.m

# wggambfit

## PURPOSE

Is an internal routine for wggamfit

## SYNOPSIS

l=wggambfit(b,data,F,def)

## DESCRIPTION

WGGAMBFIT Is an internal routine for wggamfit

## CROSS-REFERENCE INFORMATION

This function calls:
 wggamcdf Generalized Gamma cumulative distribution function drawnow Flush pending graphics events. gamma Gamma function. gammaln Logarithm of gamma function. mean Average or mean value. nan Not-a-Number. num2str Convert number to string. (Fast version) plot Linear plot.
This function is called by:
 wggamfit Parameter estimates for Generalized Gamma data.

## SOURCE CODE

001 function l=wggambfit(b,data,F,def)
002 %WGGAMBFIT Is an internal routine for wggamfit
003 %
004
005 % History
006 % revised pab 21.01.2004
007 % revised pab 27.01.2001
008 if nargin<4|isempty(def),def=1;end
009
010 monitor = logical(0);
011 switch def
012   case 1,
013     %ld  = log(data);
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)), % Avoid error with gammaln for a<0 pab 27.01.2001
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, % LS-fit to empirical CDF
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,% Moment fit: data = E(x^2)/E(x)^2,F= E(x^3)/E(x^2)^(3/2)
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,% Moment fit: data = E(x^3)/E(x^2)^(3/2),F= E(x^4)/E(x^2)^(2)
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

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