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
Check if all input arguments are either scalar or of common size. | |
Numerically evaluates a integral using a Gauss quadrature. | |
2D Weibull probability density function (pdf). | |
Weibull cumulative distribution function | |
Mean and variance for the Weibull distribution. | |
Display message and abort function. | |
Convert number to string. (Fast version) | |
usage: [int tol] = quad2dg('Fun',xlow,xhigh,ylow,yhigh) |
Plot conditional empirical CDF of X1 given X2=x2 | |
Inverse of the conditional 2D weibull cdf of X2 given X1. |
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 % 035 % See also weib2dpdf, wweibpdf, quad2dg, gaussq 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 043 % - added example 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) 087 [y(k) tol1(k)]=quad2dg('weib2dpdf',0,x1s(k),0 , x2s(k),tol/2,param,condon); 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(:)), 115 %[y(ix) tol1(ix)]=quadg('weib2dpdf',0,x2(ix),tol,[],x1(ix),param([3:4 1:2 5]),2); 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(:)), 127 % [y(ix) tol1(ix)]=quadg('weib2dpdf',0,x1(ix),tol,[],x2(ix),param,2); 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
Comments or corrections to the WAFO group