001 function [x, k2]= weib2dcinv(x1,p,a1,b1,a2,b2,c12,tol);
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
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
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
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
064
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
075 cvar=linspace(0,max(x1(k)),20);
076
077 [mn v ]=weib2dstat(parm,1,cvar);
078 mn = interp1(cvar,mn,x1(k),'linear');
079 v = interp1(cvar,v,x1(k),'linear');
080
081 switch 2,
082 case 1,
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,
088 xk=abs(wnorminv(pk,mn,v));
089
090 if any(isnan(xk))
091 disp(['Warning: NaNs ' num2str(sum(isnan(xk)))])
092 end
093 case 3,
094
095
096
097 end
098
099
100 h = ones(size(pk));
101
102
103
104
105
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
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
124
125
126
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
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