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)
Check if all input arguments are either scalar or of common size. | |
Inverse Gaussian cumulative distribution function | |
Inverse Gaussian probability density function | |
Mean and variance for the Inverse Gaussian distribution. | |
Inverse of the Lognormal distribution function | |
Display message and abort function. | |
Infinity. | |
Linearly spaced vector. | |
Not-a-Number. | |
Convert number to string. (Fast version) |
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
Comments or corrections to the WAFO group