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: This function is called by:

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