Home > wafo > wstats > dist2dcdf.m

# dist2dcdf

## PURPOSE

Joint 2D CDF computed as int F(X1

## SYNOPSIS

[y ,eps2] = dist2dcdf(V,H,phat,condon)

## DESCRIPTION

``` DIST2DCDF Joint 2D CDF computed as int F(X1<v|X2=x2).*f(x2)dx2

CALL:  [F tol] = dist2dcdf(x1,x2,phat,condon)

F      = cdf evaluated at x1, x2
tol    = absolute tolerance of the calculated F's, i.e, abs(int-intold)
x1,x2  = evaluation points
phat   = structure array containing
x    = cellarray of distribution parameters
dist = cellarray of strings defining the distributions of
X2 and X1 given X2, respectively. Options are:
'tgumbel', 'gumbel', 'lognormal','rayleigh','weibull',
and 'gamma'.
condon = 0 regular cdf of X1 and X2, F(X1<v,X2<h), is returned (default)
1 cdf of X2 conditioned on X1, F(X2<h|X1=v), is returned
2 cdf of X1 conditioned on X2, F(X1<v|X2=h), is returned

The size of F is the common size of X1 and X2.

Example: 2D Rayleigh, ie, f(x1|X2=x2)*f(x2)
x1=linspace(0,10)';
phat.x={[x1,exp(-0.1*x1)] 2 };
phat.dist={'rayl','rayl'};
dist2dcdf(2,1,phat)
f = dist2dpdf2(x1,x1,phat);
pdfplot(f)

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. dist2dfun is an internal function to dist2dcdf dist2dprb. gaussq Numerically evaluates a integral using a Gauss quadrature. wgamcdf Gamma cumulative distribution function wgamstat Mean and variance for the Gamma distribution. wgumbcdf Gumbel cumulative distribution function. wgumbstat Mean and variance for the Gumbel distribution. wlogncdf Lognormal cumulative distribution function wlognstat Mean and variance for the Lognormal distribution. wraylcdf Rayleigh cumulative distribution function wraylstat Mean and variance for the Rayleigh distribution. wweibcdf Weibull cumulative distribution function wweibstat Mean and variance for the Weibull distribution. error Display message and abort function. lower Convert string to lowercase. num2str Convert number to string. (Fast version) strcmp Compare strings.
This function is called by:
 dist2dcdfplot Plot conditional empirical CDF of X1 given X2=x2 recfig11 Probability of exceeding H: Model (dash); data (dots) recfig12 Probability of exceeding V: Model (dash); data (dots)

## SOURCE CODE

```001 function [y ,eps2] = dist2dcdf(V,H,phat,condon)
002 %DIST2DCDF Joint 2D CDF computed as int F(X1<v|X2=x2).*f(x2)dx2
003 %
004 % CALL:  [F tol] = dist2dcdf(x1,x2,phat,condon)
005 %
006 %   F      = cdf evaluated at x1, x2
007 %   tol    = absolute tolerance of the calculated F's, i.e, abs(int-intold)
008 %   x1,x2  = evaluation points
009 %   phat   = structure array containing
010 %            x    = cellarray of distribution parameters
011 %            dist = cellarray of strings defining the distributions of
012 %                  X2 and X1 given X2, respectively. Options are:
013 %                  'tgumbel', 'gumbel', 'lognormal','rayleigh','weibull',
014 %                   and 'gamma'.
015 %   condon = 0 regular cdf of X1 and X2, F(X1<v,X2<h), is returned (default)
016 %            1 cdf of X2 conditioned on X1, F(X2<h|X1=v), is returned
017 %            2 cdf of X1 conditioned on X2, F(X1<v|X2=h), is returned
018 %
019 % The size of F is the common size of X1 and X2.
020 %
021 % Example: 2D Rayleigh, ie, f(x1|X2=x2)*f(x2)
022 %   x1=linspace(0,10)';
023 %   phat.x={[x1,exp(-0.1*x1)] 2 };
024 %   phat.dist={'rayl','rayl'};
025 %   dist2dcdf(2,1,phat)
026 %   f = dist2dpdf2(x1,x1,phat);
027 %   pdfplot(f)
028 %
030
031 % tested on: matlab 5.2
032 % history:
033 %  Per A. Brodtkorb 28.10.98
034
035
036 if (nargin < 3),
037   error('Requires three input arguments.');
038 end
039
040 if (nargin <4)| isempty(condon),
041   condon=0;
042 end
043
044 CDIST=phat.dist{1}; %conditional distribution
045 UDIST=phat.dist{2}; %unconditional distribution
046 PH=phat.x{2};
047
048 [errorcode V H ] = comnsize(V,H);
049 if  errorcode > 0
050   error('Requires non-scalar arguments to match in size.');
051 end
052
053 y = zeros(size(V));
054 eps2=y;
055
056 if strcmp('gu', lower(CDIST(1:2))),
057   if strcmp('gu',lower(UDIST(1:2))),
058     k=find(H>-inf);
059   else
060     k=find(H>=0);
061   end
062 elseif strcmp('gu',lower(UDIST(1:2))),
063   k = find(V >=0 );
064 else
065   k = find(V >= 0 & H>=0);
066 end
067
068 % NB! weibpdf must be modified to correspond to
069 % pdf=x^(b-1)/a^b*exp(-(x/a)^b)
070
071 if any(k),
072     global PHAT CONDON
073     PHAT=phat;
074     CONDON=condon;
075     eps1=1e-6;%sqrt(eps); %accuracy of the estimates
076
077     switch condon
078       case { 0,1}, % 0:no conditional CDF  1: conditional CDF given V
079     %approximate # of divisions required for resolution
080     %  H                       H
081     % int  p(h)*P(V|h)dh or % int  p(h)*p(V|h)dh
082       % 0                        0
083
084     [y(k) eps2(k) ]=gaussq('dist2dfun', 0,H(k),eps1,[],V(k));
085     if any(eps2>eps1),
086       disp(['The accuracy of the computed cdf is '  num2str(max(eps2))])
087     end
088     if condon==1,  %conditional CDF given V
089       %must normalize y(k)
090       [M1, V1]=dist1dstat(PH,UDIST);
091       hmax=M1+10*sqrt(V1); %finding a value for inf
092      %  inf
093     % int  p(h)*p(V|h)dh
094       % H
095       [tmp eps3]=gaussq('dist2dfun', H(k),hmax+H(k),eps1,[],V(k))
096       if any(eps3>eps1),
097         disp(['The accuracy of the computed cdf is '  num2str(max(eps3))])
098       end
099       y(k)=y(k)./( y(k)+tmp);
100     end
101       case 2,% conditional CDF given H
102     y(k)=dist2dfun(H(k),V(k));
103     end
104 end
105
106 function cdf1=dist1dcdffun(H,Ah,dist2 )
107    switch dist2(1:2)
108       case 'ra',  pdf1= wraylcdf(H,Ah);
109       case 'we' ,  pdf1=wweibcdf(H,Ah(1),Ah(2));
110       case 'gu' ,  pdf1=wgumbcdf(H,Ah(1),Ah(2),0);
111       case 'tg' ,  pdf1=wgumbcdf(H,Ah(1),Ah(2),1);
112       case 'ga' ,  pdf1=wgamcdf(H,Ah(1),Ah(2));
113       case 'lo' ,  pdf1=wlogncdf(H,Ah(1),Ah(2));
114       otherwise, error('unknown distribution')
115     end
116 return
117
118 function [m, v]=dist1dstat(Ah,dist2);
119 switch lower(dist2(1:2))
120       case 'ra',  [m ,v]= wraylstat(Ah);
121       case 'we' ,  [m ,v]=wweibstat(Ah(1),Ah(2));
122       case 'gu' ,  [m ,v]=wgumbstat(Ah(1),Ah(2),0);
123       case 'tg' ,  [m ,v]=wgumbstat(Ah(1),Ah(2),1);
124       case 'ga' ,  [m ,v]=wgamstat(Ah(1),Ah(2));
125       case 'lo' ,  [m ,v]=wlognstat(Ah(1),Ah(2));
126       otherwise, error('unknown distribution')
127     end
128 return
129
130```

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