001 function y = weib2dpdf(x1,x2,a1,b1,a2,b2,c12,condon)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047 error(nargchk(3,8,nargin))
048 if nargin < 3,
049 error('Requires 3 input arguments.');
050 elseif ((nargin < 7) &(prod(size((a1)))~=5)|isempty(a1)),
051 error('Requires either 7 input arguments or that input argument 3 must have 5 entries .');
052 elseif (nargin < 5) &prod(size((a1)))==5,
053 if nargin<4|isempty(b1), condon=0;else condon=b1; end
054 b1=a1(2);
055 a2=a1(3); b2=a1(4);
056 c12=a1(5); a1=a1(1);
057 elseif (nargin < 8)|isempty(condon),
058 condon=0;
059 end
060
061
062
063 [errorcode x1, x2, a1,b1,a2, b2, c12] = comnsize(x1,x2,a1,b1,a2,b2,c12);
064 if errorcode > 0
065 error('Requires non-scalar arguments to match in size.');
066 end
067 xs=size(x1);
068
069
070
071
072
073 y = zeros(xs);
074
075 ok = ((a1 > 0) .*(b1 > 0).*(a2 > 0).*(b2 > 0).*(abs(c12)<1));
076
077 k1 = find(~ok);
078 if any(k1)
079 tmp = NaN;
080 y(k1) = tmp(ones(size(k1)));
081 end
082
083 k = find((x2 > 0).*(x1 > 0) & ok);
084 if any(k),
085
086 [y(k) ierr] =besseli(0,2*c12(k).*((x1(k)./a1(k)).^(b1(k)/2).*....
087 (x2(k)./a2(k)).^(b2(k)/2))./(1-c12(k).^2),1 );
088
089 switch ierr(1),
090 case 0,
091 case 1, error('Illegal arguments.')
092 case 2, disp('Overflow. Return Inf.')
093 case 3, disp('Some loss of accuracy in argument reduction.')
094 case 4, error('Complete loss of accuracy, z or nu too large.')
095 case 5, error('No convergence. Return NaN.')
096 end
097 y(k)=log(y(k));
098
099 switch condon
100 case 0,
101 y(k)=exp(- ((x1(k)./a1(k)).^b1(k) +(x2(k)./a2(k)).^b2(k) ...
102 -abs(2*c12(k).*((x1(k)./a1(k)).^(b1(k)/2).*...
103 (x2(k)./a2(k)).^(b2(k)/2))))./(1-c12(k).^2) +y(k) ) ...
104 .*b1(k).*(x1(k)./a1(k)).^(b1(k)-1)./a1(k).*b2(k).*...
105 (x2(k)./a2(k)).^(b2(k) - 1) ./a2(k)./(1-c12(k).^2);
106 case 1,
107 y(k) =exp(- (c12(k).*(x1(k)./a1(k)).^(b1(k)/2) - ...
108 (x2(k)./a2(k)).^(b2(k)/2)).^2./(1-c12(k).^2 ) +y(k) )...
109 .*(b2(k)./(a2(k).*(1-c12(k).^2))).*(x2(k)./a2(k)).^(b2(k) - 1);
110 case 2,
111 y(k)=exp(- ((x1(k)./a1(k)).^(b1(k)/2) -c12(k).*...
112 (x2(k)./a2(k)).^(b2(k)/2)).^2./(1-c12(k).^2) +y(k) )...
113 .*(b1(k)./(a1(k).*(1-c12(k).^2))).*(x1(k)./a1(k)).^(b1(k)-1);
114 case 3,
115 y(k)=x1(k).*exp(- ((x1(k)./a1(k)).^(b1(k)/2) -c12(k).*...
116 (x2(k)./a2(k)).^(b2(k)/2)).^2./(1-c12(k).^2) +y(k) )...
117 .*(b1(k)./(a1(k).*(1-c12(k).^2))).*(x1(k)./a1(k)).^(b1(k)-1);
118 case 4,
119 y(k)=x1(k).^2.*exp(- ((x1(k)./a1(k)).^(b1(k)/2) -c12(k).*...
120 (x2(k)./a2(k)).^(b2(k)/2)).^2./(1-c12(k).^2) +y(k) )...
121 .*(b1(k)./(a1(k).*(1-c12(k).^2))).*(x1(k)./a1(k)).^(b1(k)-1);
122 otherwise , error('Illegal value for CONDON')
123 end
124
125 kc=find(isinf(y(k)));
126 if any(kc),
127 disp('Computational problem occured: returning NaNs') ;
128 y(k(kc))=NaN;
129 end
130 end
131
132
133 switch condon
134 case 0,k1 = find( (x1 == 0 & b1 < 1& x2>0)|( x2 == 0 & b2 < 1& x1>0 )| ...
135 ( x1==0 & x2 == 0 & b1 < 1& b2+b1-2<0) );
136 case 1, k1 = find( x2 == 0 & b2 < 1 );
137 case {2}, k1 = find( x1 == 0 & b1 < 1 );
138 case {3,4}, k1=[];
139 end
140 if any(k1)
141 tmp = Inf;
142 y(k1) = tmp(ones(size(k1)));
143 end
144
145
146 switch condon,
147 case 0,
148 k2 = find((x1 == 0 & b1 == 1&x2>0) );
149 if any(k2),
150 y(k2) = b2(k2) .* (x2(k2)./a2(k2)).^ (b2(k2) - ...
151 1)./a2(k2)./(1-c12(k2).^2)./a1(k2).*....
152 exp(- ( (x2(k2)./a2(k2)) .^ b2(k2) )./(1-c12(k2).^2) );
153 end
154 k3 = find((x2 == 0 & b2 == 1&x1>0) );
155 if any(k3),
156 y(k3) = b1(k3) .* (x1(k3)./a1(k3)).^ (b1(k3) - ...
157 1)./a1(k3)./a2(k3)./(1-c12(k3).^2).*...
158 exp(- ((x1(k3)./a1(k3)) .^ b1(k3) )./(1-c12(k3).^2) );
159 end
160 k4 = find( x1==0 & x2==0 & b1 == 1& b2==1) ;
161 if any(k4),
162 y(k4) = 1./(a1(k4).*a2(k4).*(1-c12(k4).^2));
163 end
164
165 case 1,
166 k3 = find(x2 == 0 & b2 == 1 );
167 if any(k3),
168 y(k3) = exp(- (c12(k3).^2.*(x1(k3)./a1(k3)) .^ b1(k3) ...
169 )./(1-c12(k3).^2) )./a2(k3)./(1-c12(k3).^2) ;
170 end
171
172 case {2},
173 k3 = find(x1 == 0 & b1 == 1 );
174 if any(k3),
175 y(k3) = exp(- (c12(k3).^2.*(x2(k3)./a2(k3)) .^ b2(k3) ...
176 )./(1-c12(k3).^2) )./a1(k3)./(1-c12(k3).^2) ;
177 end
178 case {3,4},
179 end
180
181
182
183
184
185
186