WGAMINV Inverse of the Gamma distribution function CALL: x = wgaminv(F,a,b) x = inverse cdf for the Gamma distribution evaluated at F a = parameter, a>0 b = parameter, b>0 (default b=1) Example: F = linspace(0,1,100); x = wgaminv(F,5); plot(F,x) See also wgamcdf
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. |
Inverse of the Chi squared distribution function | |
Parameter estimates for Exponential data. |
001 function x = wgaminv(F,a,b) 002 %WGAMINV Inverse of the Gamma distribution function 003 % 004 % CALL: x = wgaminv(F,a,b) 005 % 006 % x = inverse cdf for the Gamma distribution evaluated at F 007 % a = parameter, a>0 008 % b = parameter, b>0 (default b=1) 009 % 010 % Example: 011 % F = linspace(0,1,100); 012 % x = wgaminv(F,5); 013 % plot(F,x) 014 % 015 % See also wgamcdf 016 017 % Reference: Johnson, Kotz and Balakrishnan (1994) 018 % "Continuous Univariate Distributions, vol. 1", p. 494 ff 019 % Wiley 020 021 % Tested on; Matlab 5.3 022 % History: 023 % adapted from stixbox ms 26.06.2000 024 % Revised by jr 01-August-2000 025 % - Added approximation for higher degrees of freedom (see References) 026 % added b parameter ms 23.08.2000 027 % revised pab 23.10.2000 028 % - added comnsize, nargchk 029 030 error(nargchk(2,3,nargin)) 031 if nargin<3|isempty(b), b=1;end 032 033 [errorcode F a b] = comnsize(F,a,b); 034 if errorcode > 0 035 error('F, a and b must be of common size or scalar.'); 036 end 037 x=zeros(size(F)); 038 039 ok = (0<=F& F<=1 & a>0 & b>0); 040 041 042 k = find(a>130 & ok); 043 if any(k), % This approximation is from Johnson et al, p. 348, Eq. (17.33). 044 za = wnorminv(F(k),0,1); 045 x(k) = a(k) + sqrt(a(k)).*( sqrt(a(k)).*( (1-(1/9).*a(k).^(-1)+... 046 (1/3).*a(k).^(-0.5).*za(k)).^3 -1) ); 047 end 048 049 k1 = find(0<F& F<1& a<=130 & ok); 050 if any(k1) 051 a1 = a(k1); 052 x1 = max(a1-1,0.1); % initial guess 053 F1=F(k1); 054 dx = ones(size(x1)); 055 iy=0; 056 max_count=100; 057 ix =find(x1); 058 while (any(ix) & iy<max_count) 059 xi=x1(ix);ai=a1(ix); 060 dx(ix) = (wgamcdf(xi,ai) - F1(ix)) ./wgampdf(xi,ai); 061 xi = xi -dx(ix); 062 % Make sure that the current guess is larger than zero. 063 x1(ix) = xi + 0.5*(dx(ix) - xi).* (xi<=0); 064 iy=iy+1; 065 ix=find((abs(dx) > sqrt(eps)*abs(x1)) & abs(dx) > sqrt(eps)); 066 %disp(['Iteration ',num2str(iy),... 067 %' Number of points left: ' num2str(length(ix)) ]), 068 end 069 x(k1)=x1; 070 if iy == max_count, 071 disp('Warning: wgaminv did not converge.'); 072 str = 'The last step was: '; 073 outstr = sprintf([str,'%13.8f'],dx(ix)); 074 fprintf(outstr); 075 end 076 end 077 078 079 080 k2=find(F==1&ok); 081 if any(k2) 082 tmp=inf; 083 x(k2)=tmp(ones(size(k2))); 084 end 085 086 x=x.*b; 087 088 k3=find(~ok); 089 if any(k3) 090 tmp=NaN; 091 x(k3)=tmp(ones(size(k3))); 092 end 093 094 095 096 097 098 099 100 101 102 103 104 105 106 107 108 109 110
Comments or corrections to the WAFO group