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: 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