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: This function is called by:

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