Home > wafo > wstats > mdist2dstat.m

mdist2dstat

PURPOSE ^

Mean and variance for the MDIST2D distribution.

SYNOPSIS ^

[M ,V ,Tm ,cvar, tolm, tolv] = mdist2dstat(phat,condon,cvar,tol)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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