Home > wafo > wstats > wggaminv.m

# wggaminv

## PURPOSE

Inverse of the Generalized Gamma distribution function

## SYNOPSIS

x = wggaminv(F,a,b,c)

## DESCRIPTION

``` WGGAMINV Inverse of the Generalized Gamma distribution function

CALL:  x = wggaminv(F,a,b,c)

x = inverse cdf for the Generalized Gamma distribution evaluated at F
a,b,c  = parameters (see wggampdf)

Example:
x = linspace(0,1,200);
p1 = wggaminv(x,0.5,1,1); p2 = wggaminv(x,1,2,1);
p3 = wggaminv(x,2,1,2); p4 = wggaminv(x,2,2,2);
plot(x,p1,x,p2,x,p3,x,p4)```

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. wgamcdf Gamma cumulative distribution function wgampdf Gamma probability density function wnorminv Inverse of the Normal distribution function error Display message and abort function. nan Not-a-Number. sprintf Write formatted data to string.
This function is called by:

## SOURCE CODE

```001 function x = wggaminv(F,a,b,c)
002 %WGGAMINV Inverse of the Generalized Gamma distribution function
003 %
004 % CALL:  x = wggaminv(F,a,b,c)
005 %
006 %        x = inverse cdf for the Generalized Gamma distribution evaluated at F
007 %   a,b,c  = parameters (see wggampdf)
008 %
009 %
010 % Example:
011 %   x = linspace(0,1,200);
012 %   p1 = wggaminv(x,0.5,1,1); p2 = wggaminv(x,1,2,1);
013 %   p3 = wggaminv(x,2,1,2); p4 = wggaminv(x,2,2,2);
014 %   plot(x,p1,x,p2,x,p3,x,p4)
015
016 % Reference: Cohen & Whittle, (1988) "Parameter Estimation in Reliability
017 % and Life Span Models", p. 220 ff, Marcel Dekker.
018
019 % Tested on; Matlab 5.3
020 % History:
021 % adapted from stixbox ms 10.08.2000
022 % revised pab 23.10.2000
023 %  - added comnsize, nargchk
024
025
026
027 error(nargchk(2,4,nargin))
028 if nargin<3|isempty(b), b=1;end
029 if nargin<3|isempty(c), c=1;end
030
031 [errorcode F a b c] = comnsize(F,a,b,c);
032 if errorcode > 0
033   error('F, a b and c must be of common size or scalar.');
034 end
035
036 x=zeros(size(F));
037
038 ok = (0<=F& F<=1 & a>0 & b>0 & c>0);
039
040
041 k = find(a>130 & ok);
042 if any(k),  % This approximation is from Johnson et al, p. 348, Eq. (17.33).
043   za = wnorminv(F(k),0,1);
044   x(k) = a(k) + sqrt(a(k)).*( sqrt(a(k)).*( (1-(1/9).*a(k).^(-1)+...
045       (1/3).*a(k).^(-0.5).*za(k)).^3 -1) );
046 end
047
048 k1 = find(0<F& F<1& a<=130 & ok);
049 if any(k1)
050   a1 = a(k1);
051   x1 = max(a(k1)-1,0.1); % initial guess
052   F1=F(k1);
053   dx = ones(size(x1));
054   iy=0;
055   max_count=100;
056   ix =find(x1);
057   while (any(ix) & iy<max_count)
058     xi=x1(ix);ai=a1(ix);
059     dx(ix) = (wgamcdf(xi,ai) - F1(ix)) ./wgampdf(xi,ai);
060     xi = xi -dx(ix);
061     % Make sure that the current guess is larger than zero.
062     x1(ix) = xi + 0.5*(dx(ix) - xi).* (xi<=0);
063     iy=iy+1;
064     ix=find((abs(dx) > sqrt(eps)*abs(x1))  &  abs(dx) > sqrt(eps));
065     %disp(['Iteration ',num2str(iy),...
066     %'  Number of points left:  ' num2str(length(ix)) ]),
067   end
068   x(k1)=x1;
069   if iy == max_count,
070     disp('Warning: wgaminv did not converge.');
071     str = 'The last step was:  ';
072     outstr = sprintf([str,'%13.8f'],dx(ix));
073     fprintf(outstr);
074   end
075 end
076
077
078
079 k2=find(F==1&ok);
080 if any(k2)
081   tmp=inf;
082   x(k2)=tmp(ones(size(k2)));
083 end
084
085 x=x.^(1./c).*b;
086
087 k3=find(~ok);
088 if any(k3)
089   tmp=NaN;
090   x(k3)=tmp(ones(size(k3)));
091 end
092 return
093
094
095
096```

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