Home > wafo > wstats > weib2dstat.m

# weib2dstat

## PURPOSE

Mean and variance for the 2D Weibull distribution

## SYNOPSIS

[M ,V ,Tm ,cvar, tolm, tolv] = weib2dstat(a1,b1,a2,b2,c12,condon,cvar,tol)

## DESCRIPTION

  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

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. gaussq Numerically evaluates a integral using a Gauss quadrature. hypgf Hypergeometric function F(a,b,c,x) wweibstat Mean and variance for the Weibull distribution. ellipke Complete elliptic integral. error Display message and abort function. fmin fminbnd Scalar bounded nonlinear function minimization. gamma Gamma function. linspace Linearly spaced vector. nan Not-a-Number. str2num Convert string matrix to numeric array. ver MATLAB, Simulink and toolbox version information.
This function is called by:
 weib2dcinv Inverse of the conditional 2D weibull cdf of X2 given X1. weib2dstatplot Computes and plots the conditional mean and standard deviation

## SOURCE CODE

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

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