Home > wafo > wstats > mdist2dlike.m

mdist2dlike

PURPOSE ^

MDIST log-likelihood function.

SYNOPSIS ^

[LL, cov] = mdist2dlike(params,data1,data2,dist)

DESCRIPTION ^

 MDIST2DLIKE MDIST log-likelihood function.
 
  CALL:  [L, cov] = mdist2dlike(params,x1,x2,dist)
 
       L = the MDIST log-likelihood 
    cov  = Asymptotic covariance matrix of phat (if phat is estimated by
                   a maximum likelihood method).
  params = [phat.x{:}] is the distribution parameters 
   x,x2  = data vector or matrices with common size.
   dist  = list of marginal distributions of x1 and x2, respectively 
           Options are: 'tgumbel', 'gumbel', 
           'lognormal','rayleigh','weibull','gamma'.
 
    MDIST2DLIKE is a utility function for maximum likelihood estimation. 
 
  Example: 
   [L, C] = mdist2dlike([1 2 2 10],x1,x2,{'weibull','rayleigh'})
 
  See also   mdist2dfit, mdist2dpdf, mdist2dcdf, mdist2drnd

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [LL, cov] = mdist2dlike(params,data1,data2,dist)
002 %MDIST2DLIKE MDIST log-likelihood function.
003 %
004 % CALL:  [L, cov] = mdist2dlike(params,x1,x2,dist)
005 %
006 %      L = the MDIST log-likelihood 
007 %   cov  = Asymptotic covariance matrix of phat (if phat is estimated by
008 %                  a maximum likelihood method).
009 % params = [phat.x{:}] is the distribution parameters 
010 %  x,x2  = data vector or matrices with common size.
011 %  dist  = list of marginal distributions of x1 and x2, respectively 
012 %          Options are: 'tgumbel', 'gumbel', 
013 %          'lognormal','rayleigh','weibull','gamma'.
014 %
015 %   MDIST2DLIKE is a utility function for maximum likelihood estimation. 
016 %
017 % Example: 
018 %  [L, C] = mdist2dlike([1 2 2 10],x1,x2,{'weibull','rayleigh'})
019 %
020 % See also   mdist2dfit, mdist2dpdf, mdist2dcdf, mdist2drnd
021 %
022 
023 
024 %  tested on: matlab 5.2
025 % history
026 % revised  pab 03.11.2000
027 % - improved the calculation of cov
028 % revised pab 8.11.1999
029 %  - updated header info
030 %  by Per A. Brodtkorb 01.02.99 
031 
032 
033 if nargin < 3, 
034     error('Requires at least FOUR input arguments'); 
035 end
036 
037 if (nargin< 4)|isempty(dist), 
038   error('Too few inputs')
039 else
040   HDIST=lower(dist{2});
041   VDIST=lower(dist{1});
042 end
043 
044 
045 
046 switch VDIST(1:2),
047   case 'ra', nv=1;
048  otherwise, nv=2;
049 end
050 switch HDIST(1:2),
051   case 'ra', nh=1;
052  otherwise, nh=2;
053 end
054 
055 if nv+nh+1~=length(params)
056  error('param is not the right size')
057 end
058  
059 data1=data1(:);
060 data2=data2(:);
061 [n, m] = size(data1);
062 [n2, m2] = size(data2);
063 if n~=n2
064   error('data1 and data2  must have equal size')
065 end
066 
067 if nargout == 2 & max(m,n) == 1
068   error('To compute the 2nd output, the 2nd input must have at least two elements.');
069 end
070 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
071 phat.dist=dist;
072 
073 x = mdist2dpdf(data1,data2,phat,0)+eps;
074 LL = -sum(log(x)); % log likelihood function
075 
076 if nargout > 1
077   np=nv+nh+1; % # of parameters we estimate
078   sparam=params;
079   delta = eps^.4;
080   delta2=delta^2;
081 
082   dist0 ='mdist2dpdf';
083   % Approximate 1/(nE( (d L(x|theta)/dtheta)^2)) with
084   %             1/(d^2 L(theta|x)/dtheta^2) 
085   %  using central differences
086     
087   H = zeros(np);             % Hessian matrix
088   for ix=1:np,
089     sparam = params;
090     sparam(ix)= params(ix)+delta;
091     phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
092     x  = feval(dist0,data1,data2,phat)+eps; 
093     fp = sum(log(x));
094     sparam(ix) = params(ix)-delta;
095     phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
096     x  = feval(dist0,data1,data2,phat)+eps; 
097     fm = sum(log(x));
098     H(ix,ix) = (fp+2*LL+fm)/delta2;
099     for iy=ix+1:np,
100       sparam(ix) = params(ix)+delta;
101       sparam(iy) = params(iy)+delta;
102       phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
103       
104       x   = feval(dist0,data1,data2,phat)+eps; 
105       fpp = sum(log(x));
106       sparam(iy) = params(iy)-delta;
107       phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
108       x   = feval(dist0,data1,data2,phat)+eps; 
109       fpm = sum(log(x));
110       sparam(ix) = params(ix)-delta;
111       phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
112       x   = feval(dist0,data1,data2,phat)+eps; 
113       fmm = sum(log(x));
114       sparam(iy) = params(iy)+delta;
115       phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
116       x   = feval(dist0,data1,data2,phat)+eps; 
117       fmp = sum(log(x));
118       H(ix,iy) = (fpp-fmp-fpm+fmm)/(4*delta2);
119       H(iy,ix) = H(ix,iy);
120     end
121   end
122   % invert the Hessian matrix (i.e. invert the observed information number)
123   cov = -H\eye(np); 
124   
125 end
126

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