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)

## 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:
 wchi2inv Inverse of the Chi squared distribution function wexpfit Parameter estimates for Exponential data.

## 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 %
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