MDIST2DSTAT Mean and variance for the MDIST2D distribution. CALL: [M,V] = mdist2dstat(phat,condon,cvar, tol) M,V = mean and variance of the mdist2d distribution phat = parameter structure array (see mdist2dfit) condon = 0 returns the mean, covariance and modal value of x1 and x2 (default) 1 conditional mean, variance and modal value of x2 given x1 2 conditional mean, variance and modal value of x1 given x2 cvar = conditional variable i.e., x1 or x2 depending on condon. tol = relative tolerance (default 1e-3) Example x1=linspace(0,10)'; phat.x={1 2 10}; phat.dist={'rayl','rayl'}; [M,V]=mdist2dstat(phat,2,x1); plot(x1,M,'r--',x1,sqrt(V),'k-') title(' Conditional mean and standard deviation') legend('E(x1|x2)','std(x1|x2)') xlabel('x2') See also mdist2dpdf, mdist2dfit mdist2drnd
Numerically evaluates a integral using a Gauss quadrature. | |
Joint 2D PDF due to Plackett given as f{x1}*f{x2}*G(x1,x2;Psi). | |
Mean and variance for the Gamma distribution. | |
Mean and variance for the Gumbel distribution. | |
Mean and variance for the Lognormal distribution. | |
Mean and variance for the Rayleigh distribution. | |
Mean and variance for the Weibull distribution. | |
Display message and abort function. | |
Scalar bounded nonlinear function minimization. | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Convert string matrix to numeric array. | |
MATLAB version number. |
Inverse of the conditional cdf of X2 given X1. | |
Computes and plots the conditional mean and standard deviation. |
001 function [M ,V ,Tm ,cvar, tolm, tolv] = mdist2dstat(phat,condon,cvar,tol) 002 %MDIST2DSTAT Mean and variance for the MDIST2D distribution. 003 % 004 % CALL: [M,V] = mdist2dstat(phat,condon,cvar, tol) 005 % 006 % M,V = mean and variance of the mdist2d distribution 007 % phat = parameter structure array (see mdist2dfit) 008 % condon = 0 returns the mean, covariance and modal value of x1 and x2 (default) 009 % 1 conditional mean, variance and modal value of x2 given x1 010 % 2 conditional mean, variance and modal value of x1 given x2 011 % cvar = conditional variable i.e., x1 or x2 depending on condon. 012 % tol = relative tolerance (default 1e-3) 013 % 014 % Example 015 % x1=linspace(0,10)'; 016 % phat.x={1 2 10}; 017 % phat.dist={'rayl','rayl'}; 018 % [M,V]=mdist2dstat(phat,2,x1); 019 % plot(x1,M,'r--',x1,sqrt(V),'k-') 020 % title(' Conditional mean and standard deviation') 021 % legend('E(x1|x2)','std(x1|x2)') 022 % xlabel('x2') 023 % 024 % See also mdist2dpdf, mdist2dfit mdist2drnd 025 026 % references: 027 % Plackett, R. L. (1965) "A class of bivariate distributions." 028 % J. Am. Stat. Assoc. 60. 516-22 029 % [1] Michel K. Ochi, 030 % OCEAN TECHNOLOGY series 6 031 % "OCEAN WAVES, The stochastic approach", Cambridge 032 % 1998 p. 133-134. 033 034 035 % tested on: matlab 5.2 036 % history 037 038 % revised pab jan2004 039 % - added todo comment 040 % revised pab 8.11.1999 041 % - updated header info 042 % - changed phat from vector to structure 043 % Per A. Brodtkorb 28.01.99 044 045 % TODO % modal value not implemented yet 046 047 if nargin < 1, 048 error('Requires 1 input argument.'); 049 end 050 051 if (nargin < 2)|isempty(condon), condon=0;end 052 if (nargin < 3)|isempty(cvar), cvar=[]; end 053 if (nargin <4)|isempty(tol), tol=1e-3; end; 054 VDIST=lower(phat.dist{1}); 055 HDIST=lower(phat.dist{2}); 056 057 psi=phat.x{3}; 058 PV=phat.x{1}; 059 PH=phat.x{2}; 060 061 062 [m1 v1]=dist1dstatfun(PV,VDIST);% marginal mean and covariance 063 [m2 v2]=dist1dstatfun(PH,HDIST);% marginal mean and covariance 064 065 066 % modal value not implemented yet 067 % 068 if 0 069 % This is a trick to get the html documentation correct. 070 k = mdist2dpdf(1,1,2,3); 071 end 072 073 switch condon 074 case 0, % mean and covariance 075 covar=sqrt(v1*v2)*getrho(psi); 076 077 M=[m1 m2]; 078 V=[v1 covar;covar' v2 ]; 079 %Tm=[tm1,tm2]; 080 Tm=[]; 081 case {1,2}, 082 083 if condon==1,%conditional mean and variance given V 084 txt='H'; vr=v2; mn=m2; vc=v1;mc=m1; 085 phat2=phat; 086 phat2.x=phat.x([2 1 3]); 087 phat2.dist=phat.dist([2 1]); 088 089 else%conditional mean and variance given H 090 txt='V'; vr=v1; mn=m1; vc=v2;mc=m2; 091 phat2=phat; 092 end 093 if isempty(cvar),cvar=linspace(0 ,mn+3*sqrt(vr),30)'; end 094 if 1 095 xinf=mn+mn./mc.*cvar +15*sqrt(vr); % infinite value for V or H 096 %tmp=input(['Enter an infinite value for ', txt, ': (' , num2str(xinf), ')']); 097 %if ~isempty(tmp), xinf=tmp;end 098 %disp(['Infinite value for ' ,txt,' is set to ' num2str(xinf(:)')]) 099 end 100 101 M=zeros(size(cvar));V=M;Tm=M; tolm=M;tolv=V; 102 %do the integration with a quadrature formulae 103 %E(H|V) or E(V|H) 104 [M tolm]=gaussq('mdist2dpdf',0,xinf,tol,[],cvar,phat2,3); 105 %E(H^2|V) or E(V^2|H) 106 [V tolv]=gaussq('mdist2dpdf',0,xinf,tol,[],cvar,phat2,4); 107 V=V-M.^2;%var(H|V) or var(V|H) 108 k=find(xinf<M+6*sqrt(V)); 109 if any(k), 110 xinf(k)=M(k)+6*sqrt(V(k))+xinf(k); % update the infinite value 111 disp(['Changed Infinite value for ', txt]);%, ' to ', num2str(xinf(k))]) 112 end 113 114 if nargout>2 115 Nint=length(cvar(:)); 116 mvrs=version;ix=find(mvrs=='.'); 117 if str2num(mvrs(1:ix(2)-1))>5.2, 118 119 for ix=1:Nint, 120 Tm(ix)=fminbnd('-mdist2dpdf(x,P1,P2,P3,P4)',... 121 max(0,M(ix)-sqrt(V(ix))),M(ix)+sqrt(V(ix)),[],... 122 cvar(ix),phat2,2 ); %modalvalue 123 end 124 else 125 126 for ix=1:Nint, 127 Tm(ix)=fmin('-mdist2dpdf(x,P1,P2,P3,P4)',... 128 max(0,M(ix)-sqrt(V(ix))),M(ix)+sqrt(V(ix)),[],... 129 cvar(ix),phat2,2 ); %modalvalue 130 end 131 end 132 133 end 134 135 end 136 return 137 138 139 function [me, va]=dist1dstatfun(Ah,dist2 ) 140 switch dist2(1:2) 141 case 'ra', [me va] = wraylstat(Ah(1)); 142 case 'we' , [me va]=wweibstat(Ah(1),Ah(2)); 143 case 'gu' , [me va]=wgumbstat(Ah(1),Ah(2),0); 144 case 'tg' , [me va]=wgumbstat(Ah(1),Ah(2),1); 145 case 'ga' , [me va]=wgamstat(Ah(1),Ah(2)); 146 case 'lo' , [me va]=wlognstat(Ah(1),Ah(2)); 147 otherwise, error('unknown distribution') 148 end 149 return 150 151 function r=getrho(psi) 152 r=zeros(size(psi)); 153 k1=find(psi==0); 154 if any(k1), 155 r(k)=-1; 156 end 157 k3=find((psi<1.*psi>0)|psi>1); 158 if any(k3), 159 r(k3)=(psi(k3).^2-1-2.*psi(k3).*log(psi(k3)))./((psi(k3)-1).^2) ; 160 end 161 return 162
Comments or corrections to the WAFO group