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: 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