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.

## CROSS-REFERENCE INFORMATION

This function calls:
 weib2dpdf 2D Weibull probability density function (pdf). error Display message and abort function.
This function is called by:
 weib2dfit Parameter estimates for 2D Weibull data.

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