Home > wafo > wstats > winvginv.m

# winvginv

## PURPOSE

Inverse of the Inverse Gaussian distribution function

## SYNOPSIS

x = winvginv(F,m,l)

## DESCRIPTION

``` WINVGINV Inverse of the Inverse Gaussian distribution function

CALL:  x = winvginv(F,m,l)

x = inverse cdf for the Inverse Gaussian distribution evaluated at F
m,l  = parameters (see winvgpdf)

Example:
x = linspace(0,1,200);
p1 = winvginv(x,1,1); p2 = winvginv(x,1,.25);
plot(x,p1,x,p2)```

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. winvgcdf Inverse Gaussian cumulative distribution function winvgpdf Inverse Gaussian probability density function winvgstat Mean and variance for the Inverse Gaussian distribution. wlogninv Inverse of the Lognormal distribution function error Display message and abort function. inf Infinity. linspace Linearly spaced vector. nan Not-a-Number. num2str Convert number to string. (Fast version)
This function is called by:

## SOURCE CODE

```001 function x = winvginv(F,m,l)
002 %WINVGINV Inverse of the Inverse Gaussian distribution function
003 %
004 % CALL:  x = winvginv(F,m,l)
005 %
006 %        x = inverse cdf for the Inverse Gaussian distribution evaluated at F
007 %     m,l  = parameters (see winvgpdf)
008 %
009 % Example:
010 %   x = linspace(0,1,200);
011 %   p1 = winvginv(x,1,1); p2 = winvginv(x,1,.25);
012 %   plot(x,p1,x,p2)
013
014 % Reference:
015 % Cohen & Whittle, (1988) "Parameter Estimation in Reliability
016 % and Life Span Models", p. 259 ff, Marcel Dekker.
017
018
019 % Tested on; Matlab 5.3
020 % History:
021 % adapted from stixbox ms 14.08.2000
022 % ad hoc solutions to improve convergence added ms 15.08.2000
023 % revised pab 25.10.2000
024 % - added nargchk + comnsize
025 % - changed the ad hoc solutions a bit
026
027 error(nargchk(3,3,nargin))
028 %if nargin<2|isempty(m),  m=0;  end
029 %if nargin<3|isempty(l),  l=1;  end
030
031 [errorcode, F, m, l] = comnsize (F,m, l);
032 if (errorcode > 0)
033   error ('F, m and l must be of common size or scalar');
034 end
035
036 x=zeros(size(F));
037 ok=(F>=0 & F<=1 &(m>0&(l>0)));
038
039 k=find((F>0)&(F<1) & ok);
040 if any(k)
041   mk=m(k);lk=l(k);  Fk=F(k);
042   % Need  a better starting guess here for
043   %  1) l<1 and F close to 1
044   %  2) l>5 and F close to 1 or 0
045   % Supply a starting guess
046   [mk,v]=winvgstat(mk,lk);
047   temp = log(v + mk .^ 2);
048   mu = 2 * log(mk) - 0.5 * temp;
049   sa = -2 * log(mk) + temp;
050   xk = wlogninv(Fk,mu,sa);
051   if 0, % old starting guess
052     x1=linspace(0.01*v^(1/2),4*v^(1/2),100);
053     F1=winvgcdf(x1,mk,lk);
054     for i=1:length(Fk);
055       k3=find(Fk(i)<F1);
056       if ~isempty(k3)
057     xstart(i)=x1(min(k3));
058       else
059     xstart(i)=4*v^(1/2);
060       end
061     end
062   end
063
064
065   dx = ones(size(xk));
066   iy=0;
067   ix =find(xk);
068   count=1;  max_count=100;
069   mstep=2*sqrt(v); % maximum step in each iteration
070   while (any(ix)&(count<max_count));
071     count=count+1;
072     xi=xk(ix);mi=mk(ix);li=lk(ix);
073     dx(ix) = (winvgcdf(xi,mi,li) - Fk(ix))./winvgpdf(xi,mi,li);
074     dx(ix) = min(abs(dx(ix)),mstep(ix)/count).*sign(dx(ix));
075     xi = xi - dx(ix);
076     % Make sure that the current guess is larger than zero.
077     xk(ix) = xi + 0.5*(dx(ix) - xi) .* (xi<=0);
078
079     ix=find((abs(dx) > sqrt(eps)*abs(xk))  &  abs(dx) > sqrt(eps));
080     disp(['Iteration ',num2str(count),...
081     '  Number of points left:  ' num2str(length(ix)) ]),
082   end
083
084   if (count==max_count)
085     disp('Warning: WINVGINV did not converge!')
086     disp(['The last steps were all less than: ' num2str(max(abs(dx(ix))))])
087     k(ix) % index to
088   end
089
090   x(k)=xk;
091 end
092
093 k1=find(F==1&ok);
094 if any(k1)
095   tmp=Inf;
096   x(k1)=tmp(ones(size(k1)));
097 end
098
099 k2=find(~ok);
100 if any(k1)
101   tmp=NaN;
102   x(k2)=tmp(ones(size(k2)));
103 end
104
105
106
107
108```

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