Home > wafo > wstats > weib2dpdf.m

# weib2dpdf

## PURPOSE

2D Weibull probability density function (pdf).

## SYNOPSIS

y = weib2dpdf(x1,x2,a1,b1,a2,b2,c12,condon)

## DESCRIPTION

``` WEIB2DPDF 2D Weibull probability density function (pdf).

CALL:  f = weib2dpdf(x1,x2,param,condon)

f    = PDF evaluated at x1 and x2.
x1,x2 = coordinates
param = [ A1, B1,A2,B2,C12] the parameters of the distribution
condon = 0 it returns the the regular pdf of X1 and X2 (default)
1 it returns the conditional pdf given X1
2 it returns the conditional pdf given X2

The size of f is the common size of the input arguments X1 and X2. A scalar input
functions as a constant matrix of the same size as the other inputs.
The PDF is defined by:

f(X1,X2)=B1*B2*xn1^(B1-1)*xn2^(B2-1)/A1/B1/N*...
exp{-[xn1^B1 +xn2^B2 ]/N }*I0(2*C12*xn1^(B1/2)*xn2^(B2/2)/N)
where
N=1-C12^2, xn1=X1/A1,  xn2=X2/A2 and
I0 is the modified  bessel function of zeroth order.

(The marginal distribution of X1 is weibull with parameters A1 and B1)
(The marginal distribution of X2 is weibull with parameters A2 and B2)
(C12 is the interaction parameter between X1 and X2)

Example:
x = linspace(0,6,200); [X1 X2]=meshgrid(x);
f = weib2dpdf(X1,X2,[2 2  3 2.5 .8]);
mesh(x,x,f),
figure(2), contour(x,x,f)

See also  weib2dcdf, wweibpdf, besseli```

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. besseli Modified Bessel function of the first kind. error Display message and abort function. inf Infinity. nan Not-a-Number.
This function is called by:
 weib2dcdf Joint 2D Weibull cumulative distribution function weib2dcinv Inverse of the conditional 2D weibull cdf of X2 given X1. weib2dlike 2D Weibull log-likelihood function. weib2dpdf2 Joint 2D Weibull probability density function

## SOURCE CODE

```001 function y = weib2dpdf(x1,x2,a1,b1,a2,b2,c12,condon)
002 %WEIB2DPDF 2D Weibull probability density function (pdf).
003 %
004 %  CALL:  f = weib2dpdf(x1,x2,param,condon)
005 %
006 %     f    = PDF evaluated at x1 and x2.
007 %    x1,x2 = coordinates
008 %    param = [ A1, B1,A2,B2,C12] the parameters of the distribution
009 %   condon = 0 it returns the the regular pdf of X1 and X2 (default)
010 %            1 it returns the conditional pdf given X1
011 %            2 it returns the conditional pdf given X2
012 %
013 %   The size of f is the common size of the input arguments X1 and X2. A scalar input
014 %   functions as a constant matrix of the same size as the other inputs.
015 %   The PDF is defined by:
016 %
017 %  f(X1,X2)=B1*B2*xn1^(B1-1)*xn2^(B2-1)/A1/B1/N*...
018 %           exp{-[xn1^B1 +xn2^B2 ]/N }*I0(2*C12*xn1^(B1/2)*xn2^(B2/2)/N)
019 %   where
020 %      N=1-C12^2, xn1=X1/A1,  xn2=X2/A2 and
021 %      I0 is the modified  bessel function of zeroth order.
022 %
023 %  (The marginal distribution of X1 is weibull with parameters A1 and B1)
024 %  (The marginal distribution of X2 is weibull with parameters A2 and B2)
025 %  (C12 is the interaction parameter between X1 and X2)
026 %
027 %Example:
028 %  x = linspace(0,6,200); [X1 X2]=meshgrid(x);
029 %  f = weib2dpdf(X1,X2,[2 2  3 2.5 .8]);
030 %  mesh(x,x,f),
031 %  figure(2), contour(x,x,f)
032 %
033 % See also  weib2dcdf, wweibpdf, besseli
034
035
036
037 %   References:
038 %     Dag Myrhaug & Håvard Rue
039 %  Journal of Ship Research, Vol 42, No3, Sept 1998, pp 199-205
040
041 %Tested on: matlab 5.2
042 % History:
043 %  by Per A. Brodtkorb 13.11.98
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 %if ndims(x1)>2, % make sure they have the right shape
069 %  a1=reshape(a1,xs);a2=reshape(a2,xs);b1=reshape(b1,xs);b2=reshape(b2,xs);
070 %  c12=reshape(c12,xs);
071 %end
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   %scales bessel with exp(-abs(x))  to avoid overflow
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, %computation OK
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,%regular pdf
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, %pdf conditioned on x1 ie. p(x2|x1)
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,%pdf conditioned on x2 ie. p(x1|x2)
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, % secret option  used by weib2dstat: returns x1*p(x1|x2)
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, % secret option  used by weib2dstat: returns x1^2*p(x1|x2)
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 % Special case for asymptotes.
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=[]; %do nothing
139 end
140 if any(k1)
141   tmp   = Inf;
142   y(k1) = tmp(ones(size(k1)));
143 end
144
145 % Special case when the marginal Weibull is the same as exponential.
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}, %do nothing
179 end
180
181
182
183
184
185
186```

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