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')

## CROSS-REFERENCE INFORMATION

This function calls:
 gaussq Numerically evaluates a integral using a Gauss quadrature. mdist2dpdf Joint 2D PDF due to Plackett given as f{x1}*f{x2}*G(x1,x2;Psi). wgamstat Mean and variance for the Gamma distribution. wgumbstat Mean and variance for the Gumbel distribution. wlognstat Mean and variance for the Lognormal distribution. wraylstat Mean and variance for the Rayleigh distribution. wweibstat Mean and variance for the Weibull distribution. error Display message and abort function. fmin fminbnd Scalar bounded nonlinear function minimization. linspace Linearly spaced vector. lower Convert string to lowercase. str2num Convert string matrix to numeric array. version MATLAB version number.
This function is called by:
 mdist2dcinv Inverse of the conditional cdf of X2 given X1. mdist2dstatplot Computes and plots the conditional mean and standard deviation.

## 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 %
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