WEIB2DSTAT Mean and variance for the 2D Weibull distribution CALL: [M,V,Tm] = weib2dstat( a1,b1,a2,b2,c12,condon,cvar, tol ) [M,V,Tm] = weib2dstat([ a1,b1,a2,b2,c12],condon,cvar,tol ) M , V, Tm = mean, variance and modal value, respectively ai, bi, c12 = parameters of the distribution condon = 0 Return mean, covariance and modal value of X1 and X2 (default) 1 Return the conditional values given X1 (cvar) 2 Return the conditional values given X2 (cvar) cvar = vector of conditional values (default depending on marginal mean and variance) tol = relative tolerance used in the integration. (default 1e-3) The distribution 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) Examples: [m v] = weib2dstat([2 2 3 2.5 .8]) % mean and covariance x = linspace(0,6)'; [m v] = weib2dstat([2 2 3 2.5 .8],2,x); plot(x,m,'r--', x,sqrt(v),'g-') % conditional mean and standard deviation. See also weib2dpdf, ellipke, hypgf, gamma, gaussq
Check if all input arguments are either scalar or of common size. | |
Numerically evaluates a integral using a Gauss quadrature. | |
Hypergeometric function F(a,b,c,x) | |
Mean and variance for the Weibull distribution. | |
Complete elliptic integral. | |
Display message and abort function. | |
Scalar bounded nonlinear function minimization. | |
Gamma function. | |
Linearly spaced vector. | |
Not-a-Number. | |
Convert string matrix to numeric array. | |
MATLAB, Simulink and toolbox version information. |
Inverse of the conditional 2D weibull cdf of X2 given X1. | |
Computes and plots the conditional mean and standard deviation |
001 function [M ,V ,Tm ,cvar, tolm, tolv] = weib2dstat(a1,b1,a2,b2,c12,condon,cvar,tol) 002 % WEIB2DSTAT Mean and variance for the 2D Weibull distribution 003 % 004 % CALL: [M,V,Tm] = weib2dstat( a1,b1,a2,b2,c12,condon,cvar, tol ) 005 % [M,V,Tm] = weib2dstat([ a1,b1,a2,b2,c12],condon,cvar,tol ) 006 % 007 % M , V, Tm = mean, variance and modal value, respectively 008 % ai, bi, c12 = parameters of the distribution 009 % 010 % condon = 0 Return mean, covariance and modal value of X1 and X2 (default) 011 % 1 Return the conditional values given X1 (cvar) 012 % 2 Return the conditional values given X2 (cvar) 013 % 014 % cvar = vector of conditional values 015 % (default depending on marginal mean and variance) 016 % tol = relative tolerance used in the integration. (default 1e-3) 017 % 018 % The distribution is defined by its PDF: 019 % f(X1,X2)=B1*B2*xn1^(B1-1)*xn2^(B2-1)/A1/B1/N*... 020 % exp{-[xn1^B1 +xn2^B2 ]/N }*I0(2*C12*xn1^(B1/2)*xn2^(B2/2)/N) 021 % where 022 % N=1-C12^2, xn1=X1/A1, xn2=X2/A2 and 023 % I0 is the modified bessel function of zeroth order. 024 % 025 % (The marginal distribution of X1 is weibull with parameters A1 and B1) 026 % (The marginal distribution of X2 is weibull with parameters A2 and B2) 027 % 028 % Examples: 029 % [m v] = weib2dstat([2 2 3 2.5 .8]) % mean and covariance 030 % x = linspace(0,6)'; 031 % [m v] = weib2dstat([2 2 3 2.5 .8],2,x); 032 % plot(x,m,'r--', x,sqrt(v),'g-') % conditional mean and standard deviation. 033 % 034 % See also weib2dpdf, ellipke, hypgf, gamma, gaussq 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.1 042 % history: 043 % by Per A. Brodtkorb 13.11.98 044 045 046 047 if nargin < 1, 048 error('Requires 1 input argument.'); 049 elseif (nargin < 5) &prod(size((a1)))~=5|isempty(a1), 050 error('Requires either 5 input arguments or that input argument 1 must have 5 entries .'); 051 elseif (nargin < 5) &prod(size((a1)))==5, 052 if (nargin<2)|isempty(b1), condon=0; else condon=b1; end 053 if (nargin <3)|isempty(a2), cvar=[]; else cvar=a2; end; 054 if (nargin <4)|isempty(b2), tol=1e-3; else tol=b2; end; 055 b1=a1(2); 056 a2=a1(3); b2=a1(4); c12=a1(5); a1=a1(1); 057 else 058 if (nargin < 6)|isempty(condon), condon=0;end 059 if (nargin < 7)|isempty(cvar), cvar=[]; end 060 if (nargin <8)|isempty(tol), tol=1e-3; end; 061 end 062 063 [errorcode a1 b1 a2 b2 c12] = comnsize(a1,b1,a2,b2,c12); 064 if errorcode > 0 065 error('Requires non-scalar arguments to match in size.'); 066 end 067 068 [m1,v1]=wweibstat(a1,b1); % 069 [m2,v2]=wweibstat(a2,b2); 070 tm1=zeros(size(a1));tm2=zeros(size(a2)); 071 072 073 k=find(b1>1); 074 if any(k), 075 tm1(k)=a1(k).*(1-1./b1(k)).^(1./b1(k)); 076 end 077 k2=find(b2>1); 078 if any(k2), 079 tm2(k2)=a2(k2).*(1-1./b2(k2)).^(1./b2(k2)); 080 end 081 082 switch condon 083 case 0, % mean and covariance 084 covar=zeros(size(a2)); 085 ok = a1 > 0 & b1 > 0 &a2 > 0 & b2 > 0 | abs(c12)<1; 086 k = find(~ok); 087 if any(k) 088 tmp = NaN; 089 covar(k) = tmp(ones(size(k))); 090 end 091 k = find(ok); 092 if any(k), 093 if 0, % calculating covar correct for Rayleigh 094 [K,E] = ellipke(c12(k).^2); %complete elliptic integral of first and 095 %second kind 096 covar(k) = sqrt(v1(k).*v2(k)).*(E-.5*(1-c12(k).^2).*K-pi/4)./(1-pi/4); 097 else % calculating covar correct for Weibull 098 qx=1./b1(k);qy=1./b2(k); 099 covar(k)= sqrt(v1(k).*v2(k)).*(hypgf(-qx,-qy,1,c12(k).^2)-1).* gamma(qx).*gamma(qy)./ .... 100 ( sqrt(2.*b1(k).*gamma(2*qx) - gamma(qx).^2) .* ... 101 sqrt(2.*b2(k).*gamma(2*qy) - gamma(qy).^2) ); 102 end 103 end 104 M=[m1 m2]; 105 V=[v1 covar;covar' v2 ]; 106 Tm=[tm1,tm2]; 107 case {1,2}, 108 sparam={a1 b1 a2 b2 c12}; 109 110 if condon==1,%conditional mean and variance given x1 111 txt='x2'; vr=v2; mn=m2; vc=v1;mc=m1; 112 porder=[3:4 1:2 5]; 113 else%conditional mean and variance given x2 114 txt='x1'; vr=v1; mn=m1; vc=v2;mc=m2; 115 porder=1:5; 116 end 117 if isempty(cvar),cvar=linspace(0 ,mn+3*sqrt(vr),30)'; end 118 if 1 119 120 xinf=mn+mn./mc.*cvar +15*sqrt(vr); % infinite value for x1 or x2 121 %tmp=input(['Enter an infinite value for ', txt, ': (' , num2str(xinf), ')']); 122 %if ~isempty(tmp), xinf=tmp;end 123 %disp(['Infinite value for ' ,txt,' is set to ' num2str(xinf(:)')]) 124 end 125 126 M=zeros(size(cvar));V=M;Tm=M; tolm=M;tolv=V; 127 %do the integration with a quadrature formulae 128 [M tolm]=gaussq('weib2dpdf',0,xinf,tol,[],cvar,sparam{porder},3); %E(x2|x1) or E(x1|x2) 129 [V tolv]=gaussq('weib2dpdf',0,xinf,tol,[],cvar,sparam{porder},4);%E(x2^2|x1) or E(x1^2|x2) 130 V=V-M.^2;%var(x2|x1) or var(x1|x2) 131 k=find(xinf<M+6*sqrt(V)); 132 if any(k), 133 xinf(k)=M(k)+6*sqrt(V(k))+xinf(k); % update the infinite value 134 disp(['Changed Infinite value for ', txt]);%, ' to ', num2str(xinf(k))]) 135 end 136 137 if sparam{porder(2)}>1 & nargout>2 138 %modalvalue 139 Nint=length(cvar(:)); 140 mvrs = ver; ix = find(mvrs=='.'); 141 if str2num(mvrs(1:ix(2)-1))>5.2, 142 143 for ix=1:Nint, 144 Tm(ix)=fminbnd('-weib2dpdf(x,P1,P2,P3)',max(0,M(ix)-sqrt(V(ix))),... 145 M(ix)+sqrt(V(ix)),[],cvar(ix),sparam{porder},2 ); 146 end 147 else 148 149 for ix=1:Nint, 150 Tm(ix)=fmin('-weib2dpdf(x,P1,P2,P3)',max(0,M(ix)-sqrt(V(ix))),... 151 M(ix)+sqrt(V(ix)),[],cvar(ix),sparam{porder},2 ); 152 end 153 end 154 end 155 end 156 157 158 % [M(ix) tolm(ix)]=quadg('weib2dpdf',0,xinf(ix),tol,[],cvar(ix),sparam(porder),3); %E(x2|x1) or E(x1|x2) 159 % [V(ix) tolv(ix)]=quadg('weib2dpdf',0,xinf(ix),tol,[],cvar(ix),sparam(porder),4);%E(x2^2|x1) or E(x1^2|x2) 160 %V(ix)=V(ix)-M(ix)^2; %var(x2|x1) or var(x1|x2) 161 %if Nint>ix & xinf(ix)<M(ix)+6*sqrt(V(ix)), 162 % xinf(ix+1)=xinf(ix+1)+M(ix)+10*sqrt(V(ix))-xinf(ix); % update the infinite value 163 % disp(['Changed Infinite value for ', txt, ' to ', num2str(xinf(ix))]) 164 %end
Comments or corrections to the WAFO group