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:
 comnsize Check if all input arguments are either scalar or of common size. weib2dcdf Joint 2D Weibull cumulative distribution function weib2dpdf 2D Weibull probability density function (pdf). weib2dstat Mean and variance for the 2D Weibull distribution wnorminv Inverse of the Normal distribution function error Display message and abort function. inf Infinity. interp1 1-D interpolation (table lookup) linspace Linearly spaced vector. nan Not-a-Number. num2str Convert number to string. (Fast version)
This function is called by:
 weib2drnd Random numbers from the 2D Weibull distribution.

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