Home > wafo > wstats > wgaminv.m

wgaminv

PURPOSE ^

Inverse of the Gamma distribution function

SYNOPSIS ^

x = wgaminv(F,a,b)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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