Home > wafo > wstats > weib2dlike.m

weib2dlike

PURPOSE ^

2D Weibull log-likelihood function.

SYNOPSIS ^

[LL, C] = weib2dlike(param1,data1,data2,given,gparam)

DESCRIPTION ^

  WEIB2DLIKE 2D Weibull log-likelihood function.
 
  CALL:  [L, cov] = weib2dlike(phat,data1,data2) 
 
    L           = log-likelihood of the parameters given the data
    cov         = Asymptotic covariance matrix of phat (if phat is estimated by
                  a maximum likelihood method).
    phat        = [A1 B1 A2 B2 C12] vector of distribution parameters
    data1,data2 = data vectors
 
    WEIB2DLIKE is a utility function for maximum likelihood estimation. 
    The PDF is defined by:
 
  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)/N) 
   where 
     N=1-C12^2, xn1=X1/A1,  xn2=X2/A2 and 
     I0 is the modified bessel function of zeroth order.
 
    See also  weib2dpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [LL, C] = weib2dlike(param1,data1,data2,given,gparam)
002 % WEIB2DLIKE 2D Weibull log-likelihood function.
003 %
004 % CALL:  [L, cov] = weib2dlike(phat,data1,data2) 
005 %
006 %   L           = log-likelihood of the parameters given the data
007 %   cov         = Asymptotic covariance matrix of phat (if phat is estimated by
008 %                 a maximum likelihood method).
009 %   phat        = [A1 B1 A2 B2 C12] vector of distribution parameters
010 %   data1,data2 = data vectors
011 %
012 %   WEIB2DLIKE is a utility function for maximum likelihood estimation. 
013 %   The PDF is defined by:
014 %
015 % f(X1,X2) = B1*B2*xn1^(B1-1)*xn2^(B2-1)/A1/B1/N*...
016 %            exp{-[xn1^B1 +xn2^B2 ]/N }*I0(2*C12*xn1^(B1/2)/N) 
017 %  where 
018 %    N=1-C12^2, xn1=X1/A1,  xn2=X2/A2 and 
019 %    I0 is the modified bessel function of zeroth order.
020 %
021 %   See also  weib2dpdf
022 
023 %tested on: matlab 5.1
024 % history:
025 % revised pab 1.11.2000
026 % - improoved the calculation of cov.
027 %  by Per A. Brodtkorb 14.11.98 
028 
029 % Secret options:
030 %   given       =  a vector with  Number to the given parameter: [ 1 3] means
031 %                  parameter 1 and 3 are fixed
032 %   gparam      = values of the given parameters which we consider fixed
033 
034 
035 error(nargchk(3,5,nargin))
036 
037 data1 = data1(:);
038 data2 = data2(:);
039 n  = length(data1);
040 n2 = length(data2);
041 if n~=n2
042   error('data1 and data2  must have equal size')
043 end
044 
045 sparams=zeros(1,5);
046 sparams(given)=gparam;
047 iz=1:5;iz(given)=[];
048 sparams(iz)=param1;
049 
050 x = weib2dpdf(data1,data2,sparams)+eps;
051 LL = -sum(log(x)); % log likelihood function
052 
053 if nargout > 1
054   params = sparams;
055   delta = eps^.4;
056   delta2=delta^2;
057   np=length(param1);
058   dist ='weib2dpdf';
059   % Approximate 1/(nE( (d L(x|theta)/dtheta)^2)) with
060   %             1/(d^2 L(theta|x)/dtheta^2) 
061   %  using central differences
062     
063   H = zeros(np);             % Hessian matrix
064   for ix=1:np,
065     sparam = params;
066     iw = iz(ix);
067     sparam(iw)= params(iw)+delta;
068     x  = feval(dist,data1,data2,sparam)+eps; 
069     fp = sum(log(x));
070     sparam(iw) = params(iw)-delta;
071     x  = feval(dist,data1,data2,sparam)+eps; 
072     fm = sum(log(x));
073     H(ix,ix) = (fp+2*LL+fm)/delta2;
074     for iy=ix+1:np,
075       iu = iz(iy);
076       sparam(iw) = params(iw)+delta;
077       sparam(iu) = params(iu)+delta;
078       x   = feval(dist,data1,data2,sparam)+eps; 
079       fpp = sum(log(x));
080       sparam(iu) = params(iu)-delta;
081       x   = feval(dist,data1,data2,sparam)+eps; 
082       fpm = sum(log(x));
083       sparam(iw) = params(iw)-delta;
084       x   = feval(dist,data1,data2,sparam)+eps; 
085       fmm = sum(log(x));
086       sparam(iu) = params(iu)+delta;
087       x   = feval(dist,data1,data2,sparam)+eps; 
088       fmp = sum(log(x));
089       H(ix,iy) = (fpp-fmp-fpm+fmm)/(4*delta2);
090       H(iy,ix) = H(ix,iy);
091     end
092   end
093   % invert the Hessian matrix (i.e. invert the observed information number)
094   C = -H\eye(np); 
095 end
096

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