Home > wafo > wstats > weib2dcinv.m

weib2dcinv

PURPOSE ^

Inverse of the conditional 2D weibull cdf of X2 given X1.

SYNOPSIS ^

[x, k2]= weib2dcinv(x1,p,a1,b1,a2,b2,c12,tol);

DESCRIPTION ^

 WEIB2DCINV Inverse of the conditional 2D weibull cdf of X2 given X1. 
  
  CALL: X2 =  weib2dcinv(X1 ,P,param,rtol)   
   
    X2    = inverse of the conditional weib2dcdf given X1 at the 
            probabilities in P. 
    X1    = conditional variable  
    param = [A1 B2 A2 B2 C12] distribution parameters 
  
    The size of X is the common size of the input arguments X1 and P. A scalar input   
    functions as a constant matrix of the same size as the other inputs.     
  
    WEIB2DCINV uses Newton's method to converge to the solution. 
  
  See also  weib2dcdf, weib2dstat, norminv

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [x, k2]= weib2dcinv(x1,p,a1,b1,a2,b2,c12,tol); 
002 %WEIB2DCINV Inverse of the conditional 2D weibull cdf of X2 given X1. 
003 % 
004 % CALL: X2 =  weib2dcinv(X1 ,P,param,rtol)   
005 %  
006 %   X2    = inverse of the conditional weib2dcdf given X1 at the 
007 %           probabilities in P. 
008 %   X1    = conditional variable  
009 %   param = [A1 B2 A2 B2 C12] distribution parameters 
010 % 
011 %   The size of X is the common size of the input arguments X1 and P. A scalar input   
012 %   functions as a constant matrix of the same size as the other inputs.     
013 % 
014 %   WEIB2DCINV uses Newton's method to converge to the solution. 
015 % 
016 % See also  weib2dcdf, weib2dstat, norminv 
017 % 
018  
019 % Tested on: matlab 5.2 
020 %History 
021 % revised pab 13.11.2000, minor fixes 
022 % By Per A. Brodtkorb 19.11.98 
023  
024 if nargin < 3,  
025   error('Requires 3 input arguments.');  
026 elseif (nargin < 4) &prod(size((a1)))~=5|isempty(a1), 
027   error('Requires either 5 input arguments or that input argument 1 must have 5 entries .');  
028 elseif (nargin < 5) &prod(size((a1)))==5, 
029   if (nargin <4)|isempty(b1),  tol=1e-3; else tol=b1; end; 
030    b1=a1(2); 
031   a2=a1(3);  b2=a1(4);    c12=a1(5); a1=a1(1);  
032 elseif (nargin <8)|isempty(tol),  tol=1e-3;   
033 end 
034 parm=[a1 b1 a2 b2 c12]; 
035  
036 [errorcode  x1 p a1 b1 a2 b2 c12 ] = comnsize(x1,p,a1,b1,a2, b2, c12); 
037 if errorcode > 0 
038     error('Requires non-scalar arguments to match in size.'); 
039 end 
040  
041 %   Initialize X to zero. 
042 x = zeros(size(p)); 
043  
044  ok =(0<=p & p<=1 &  a1 > 0 & b1 > 0 & a2 > 0 & b2 > 0 &  abs(c12)<1); 
045 k = find(p<0 | p>1 | ~ok); 
046 if any(k), 
047     tmp  = NaN; 
048     x(k) = tmp(ones(size(k)));  
049 end 
050  
051 % The inverse cdf of 0 is 0, and the inverse cdf of 1 is inf.   
052 k0 = find(p == 0 & ok); 
053 if any(k0), 
054     x(k0) = zeros(size(k0));  
055 end 
056  
057 k1 = find((p == 1| isinf(x1) )& ok ); 
058 if any(k1),  
059   tmp = Inf; 
060     x(k1) = tmp(ones(size(k1)));  
061 end 
062  
063 % Newton's Method 
064 % Permit no more than maxcount interations. 
065 maxcount = 100; 
066 count = 0; 
067  
068 k = find(~isinf(x1) & p > 0  &  p < 1 & ok); 
069  
070 if ~any(k),return,end 
071  
072 pk = p(k); 
073  
074 % Supply a starting guess for the iteration. 
075 cvar=linspace(0,max(x1(k)),20); % make sure that the comp. of mn and v do  
076                                 % not consume valuable time 
077 [mn v ]=weib2dstat(parm,1,cvar); %slow 
078 mn = interp1(cvar,mn,x1(k),'linear'); %fast 
079 v = interp1(cvar,v,x1(k),'linear'); %fast 
080  
081 switch 2, %method 
082   case 1,%Use a lognormal distribution.  
083   temp = log(v + mn .^ 2);  
084   mu = 2 * log(mn) - 0.5 * temp; 
085   sa = -2 * log(mn) + temp; 
086   xk = exp(wnorminv(pk,mu,sa.^2)); 
087 case 2, %   Use a normal distribution. probably the fastest choice 
088   xk=abs(wnorminv(pk,mn,v)); 
089   %xk((xk<=0))=1e-3; 
090     if any(isnan(xk)) 
091       disp(['Warning: NaNs   ' num2str(sum(isnan(xk)))]) 
092     end 
093 case 3, % use weibull distribution slowest choice 
094     %   beta=fsolve('gamma(1+1./x).^2./gamma(1+2./x)-P1.^2./(P2+P1.^2)',parm(4).*ones(size(mn)),[],[],mn,v); 
095     %   alpha=mn./gamma(1+1./beta); 
096     %   xk=wweibinv(pk,alpha,beta);  
097 end 
098 %x(k)=xk; 
099 %return 
100 h = ones(size(pk));  
101  
102 % Break out of the iteration loop for three reasons: 
103 %  1) the last update is very small (compared to x) 
104 %  2) the last update is very small (compared to sqrt(eps)) 
105 %  3) There are more than 100 iterations.  
106  
107 k2=find((abs(h) > sqrt(eps)*abs(xk))  &  abs(h) > sqrt(eps)); 
108 while(any(k2) & count < maxcount),  
109                                   
110     count = count + 1; 
111     h(k2)  = (weib2dcdf(x1(k(k2)),xk(k2),parm,1,tol) - pk(k2)) ./ weib2dpdf(x1(k(k2)),xk(k2),parm,1); 
112     
113     xnew = xk(k2) - h(k2); 
114     % Make sure that the current xnew>0. 
115     ix = find(xnew <= 0); 
116     if any(ix), 
117         xnew(ix) = xk(k2(ix)) / 10; 
118         h(k2(ix)) = xk(k2(ix))-xnew(ix); 
119     end 
120     xk(k2) = xnew; 
121      disp(['Iteration  'num2str(count) , '  Number of points left:  ' ... 
122        num2str(length(k2)) ]), 
123     %if any(isnan(xk)),  disp(['Warning: values out of range   '... 
124     %  num2str(sum(isnan(xk)))]),   end 
125      
126     % if not converged decrease the relative tolerance in the calculation of cdf 
127     if length(k2)<length(k)/4 & count>maxcount/4, tol=tol/10;,end 
128     k2=find((abs(h) > sqrt(eps)*abs(xk))  &  abs(h) > sqrt(eps));    
129 end 
130  
131  
132 % Store the converged value in the correct place 
133 x(k) = xk; 
134  
135 if count == maxcount,  
136     disp('Warning: WEIB2DCINV did not converge.'); 
137     str = ['The last steps were all less than:  ' num2str(max(abs(h(k2))))]; 
138     disp(str) 
139     ind=k(k2); 
140   else 
141     ind=[]; 
142 end 
143

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