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)
Check if all input arguments are either scalar or of common size. | |
Gamma cumulative distribution function | |
Gamma probability density function | |
Inverse of the Normal distribution function | |
Display message and abort function. | |
Not-a-Number. | |
Write formatted data to string. |
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
Comments or corrections to the WAFO group