Home > wafo > wstats > weib2dcdf.m

# weib2dcdf

## PURPOSE

Joint 2D Weibull cumulative distribution function

## SYNOPSIS

[y ,tol1] = weib2dcdf(x1,x2,param,condon,tol)

## DESCRIPTION

 WEIB2DCDF Joint 2D Weibull cumulative distribution function

CALL: [F tol]= weib2dcdf(x1,x2,param,condon,rtol)

F    = joint CDF, i.e., Prob(X1<x1 X2<x2)
tol   = estimated absolute tolerance
x1,x2 = coordinates
param = [ A1, B1,A2,B2,C12] the parameters of the distribution
condon = 0 it returns the the regular cdf of X1 and X2 (default)
1 it returns the conditional cdf given X1
2 it returns the conditional cdf given X2
rtol = specified relative tolerance (Default 1e-3)

The size of F is the common size of X1 and X2.
The CDF is defined by its PDF:

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)

Example:
x = linspace(0,6,200); [X1,X2] = meshgrid(x);
phat = [2 2  3 2.5 .8];
f = weib2dpdf(X1,X2,phat);
contour(x,x,f), hold on,
plot( [0  2  2 0 0], [0 0 1 1 0],'g-'), hold off % Mark the region
weib2dcdf(2,1,phat)  % Calculate the probability of marked region

See also  weib2dpdf, wweibpdf, quad2dg, gaussq

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. gaussq Numerically evaluates a integral using a Gauss quadrature. weib2dpdf 2D Weibull probability density function (pdf). wweibcdf Weibull cumulative distribution function wweibstat Mean and variance for the Weibull distribution. error Display message and abort function. num2str Convert number to string. (Fast version) quad2dg usage: [int tol] = quad2dg('Fun',xlow,xhigh,ylow,yhigh)
This function is called by:
 weib2dcdfplot Plot conditional empirical CDF of X1 given X2=x2 weib2dcinv Inverse of the conditional 2D weibull cdf of X2 given X1.

## SOURCE CODE

001 function [y ,tol1] = weib2dcdf(x1,x2,param,condon,tol)
002 %WEIB2DCDF Joint 2D Weibull cumulative distribution function
003 %
004 %  CALL: [F tol]= weib2dcdf(x1,x2,param,condon,rtol)
005 %
006 %     F    = joint CDF, i.e., Prob(X1<x1 X2<x2)
007 %    tol   = estimated absolute tolerance
008 %    x1,x2 = coordinates
009 %    param = [ A1, B1,A2,B2,C12] the parameters of the distribution
010 %   condon = 0 it returns the the regular cdf of X1 and X2 (default)
011 %            1 it returns the conditional cdf given X1
012 %            2 it returns the conditional cdf given X2
013 %     rtol = specified relative tolerance (Default 1e-3)
014 %
015 % The size of F is the common size of X1 and X2.
016 % The CDF is defined by its PDF:
017 %
018 %  f(X1,X2)=B1*B2*xn1^(B1-1)*xn2^(B2-1)/A1/B1/N*...
019 %           exp{-[xn1^B1 +xn2^B2 ]/N }*I0(2*C12*xn1^(B1/2)*xn2^(B2/2)/N)
020 %   where
021 %      N=1-C12^2, xn1=X1/A1,  xn2=X2/A2 and
022 %      I0 is the modified  bessel function of zeroth order.
023 %
024 %  (The marginal distribution of X1 is weibull with parameters A1 and B1)
025 %  (The marginal distribution of X2 is weibull with parameters A2 and B2)
026 %
027 % Example:
028 %   x = linspace(0,6,200); [X1,X2] = meshgrid(x);
029 %   phat = [2 2  3 2.5 .8];
030 %   f = weib2dpdf(X1,X2,phat);
031 %   contour(x,x,f), hold on,
032 %   plot( [0  2  2 0 0], [0 0 1 1 0],'g-'), hold off % Mark the region
033 %   weib2dcdf(2,1,phat)  % Calculate the probability of marked region
034 %
036
037 % tested on: matlab5.1
038 %history:
039 % revised pab Dec2003
040 %  call weibstat -> wweibstat
041 % revised pab 29.10.2000
042 % - updated to wstats
044 %  Per A. Brodtkorb 17.11.98
045
046 % NB! weibpdf must be modified to correspond to
047 % pdf=x^(b-1)/a^b*exp(-(x/a)^b)
048
049
050 error(nargchk(3,5,nargin))
051 if (nargin <5)| isempty(tol),   tol=1e-3; end
052 if (nargin <4)| isempty(condon),  condon=0;end
053 if nargin < 3,
054   error('Requires 3 input arguments.');
055 elseif prod(size((param)))~=5|isempty(param),
056    error('input argument 3 must have 5 entries.');
057 else
058   a1=param(1);  b1=param(2);
059   a2=param(3);  b2=param(4);    c12=param(5);
060 end
061
062 [errorcode x1 x2 ] = comnsize(x1,x2);
063 if  errorcode > 0
064   error('Requires non-scalar arguments to match in size.');
065 end
066 y = zeros(size(x1));
067 tol1=y;
068
069 if 0
070   % This is a trick to get the html documentation correct.
071   k = weib2dpdf(1,1,2,3);
072 end
073
074 switch condon
075   case 0,% regular cdf
076     k = find( x1 > 0 & x2>0);
077     if 1,
078        if any(k),
079      [y(k) tol1(k)]= quad2dg('weib2dpdf',0,x1(k),0 , x2(k),  tol,a1,b1,a2,b2,c12,condon);
080        end
081
082
083      else% divide the integral in to several parts this is not correct yet
084       [m1 v1]=wweibstat(a1,b1); [m2 v2]=wweibstat(a2,b2);
085        x1s=min(m1,x1); x2s=min(m2,x2);
086        if any(k)
088        end
089        k1=find( x1(k)>m1& x2(k)<=m2);
090        if any(k1),
091      [tmp1 tmp2]=quad2dg('weib2dpdf',x1s(k(k1)),  x1(k(k1)),0,   x2s(k(k1)), tol/4,param,condon);
092      y(k(k1)) =y(k(k1))+tmp1;tol1(k(k1))=tol1(k(k1))+tmp2;
093        end
094        k1=find( x1(k)<=m1& x2(k)>m2);
095        if any(k1),
096      [tmp1 tmp2]=quad2dg('weib2dpdf',0,x1s(k(k1)), x2s(k(k1)),   x2(k(k1)),   tol/4,param,condon);
097      y(k(k1)) =y(k(k1)) +tmp1;tol1(k(k1))=tol1(k(k1))+tmp2;
098        end
099        k1=find(m1<x1(k)& m2<x2(k));
100        if any(k1),
101      [tmp1 tmp2]=quad2dg('weib2dpdf',x1s(k(k1)),x1(k(k1)), x2s(k(k1)),   x2(k(k1)),   tol/4,param,condon);
102      y(k(k)) =y(k(k)) +tmp1;tol1(k(k))=tol1(k(k))+tmp2;
103        end
104      end
105    case 1,%conditional CDF given x1
106      k = find( (x1 > 0) & (x2>0) & (c12(ones(size(x2))) ~=0 ));
107      if any(k),
108        [y(k) tol1(k)]=gaussq('weib2dpdf',0,x2(k),tol,[],x1(k),a2,b2,a1,b1,c12,2);%param([3:4 1:2 5]),2);
109      end
110      k = find( (x1==0)| (c12(ones(size(x2))) ==0));
111      if any(k),
112        y(k) =wweibcdf(x2(k),param(3),param(4));
113      end
114      %for ix=1:length(x1(:)),
116      %end
117    case 2,%conditional CDF given x2
118      k = find( (x1 > 0) & (x2>0) & (c12(ones(size(x2))) ~=0));
119      if any(k),
120        [y(k) tol1(k)]=gaussq('weib2dpdf',0,x1(k),tol,[],x2(k),a1,b1,a2,b2,c12,2);%param,2);
121      end
122      k = find( (x2==0)| (c12(ones(size(x2))) ==0));
123       if any(k),
124        y(k) =wweibcdf(x1(k),param(1),param(2));
125      end
126      %for ix=1:length(x2(:)),
128      %end
129  end
130
131 % make sure 0 <= y<=1
132 k2=find(y<0);
133 if any(k2)
134   y(k2)=zeros(size(k2));
135 end
136 k2=find(y>1);
137 if any(k2)
138   y(k2)=ones(size(k2));
139 end
140 if any(isnan(y)),
141   disp(['Warning: there are : ', num2str(sum(isnan(y))),' NaNs']);
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