Home > wafo > wstats > wgevlike.m

wgevlike

PURPOSE ^

Is an internal routine for wgevfit

SYNOPSIS ^

[ll,cov] = wgevlike(p,sample)

DESCRIPTION ^

 WGEVLIKE Is an internal routine for wgevfit
         
  CALL:  [L,cov] = wgevlike(phat,sample);
 
     L    = -log(f(phat|sample)), i.e., the log-likelihood function
            with parameters phat given the data.
     cov  = Asymptotic covariance matrix of phat (if phat is estimated by
            a maximum likelihood method).
   phat   = Parameters in distribution
            [k s m] = [shape scale location].  (see wgevcdf)
   sample = the vector of data points.
 
  WGEVLIKE is a utility function for maximum likelihood estimation
  and is used by routine wgevfit.
 
  Example:
    R = wgevrnd(-0.2,3,0,1,100);                      
    phat0 = [-0.2 3 0];       % initial guess
    phat = fminsearch('wgevlike',phat0,[],R)
    [L, cov] = wgevlike(phat,R)
 
  See also  wgevfit, wgevcdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [ll,cov] = wgevlike(p,sample)  
002 %WGEVLIKE Is an internal routine for wgevfit
003 %        
004 % CALL:  [L,cov] = wgevlike(phat,sample);
005 %
006 %    L    = -log(f(phat|sample)), i.e., the log-likelihood function
007 %           with parameters phat given the data.
008 %    cov  = Asymptotic covariance matrix of phat (if phat is estimated by
009 %           a maximum likelihood method).
010 %  phat   = Parameters in distribution
011 %           [k s m] = [shape scale location].  (see wgevcdf)
012 %  sample = the vector of data points.
013 %
014 % WGEVLIKE is a utility function for maximum likelihood estimation
015 % and is used by routine wgevfit.
016 %
017 % Example:
018 %   R = wgevrnd(-0.2,3,0,1,100);                      
019 %   phat0 = [-0.2 3 0];       % initial guess
020 %   phat = fminsearch('wgevlike',phat0,[],R)
021 %   [L, cov] = wgevlike(phat,R)
022 %
023 % See also  wgevfit, wgevcdf
024 
025 % Tested  on Matlab  5.3
026 %
027 % History:
028 % Revised by PJ (Pär Johannesson) 08-Mar-2000
029 %   updated for WAFO
030 %   Needed by GEVFIT for method 'ml'.
031 % Copied from WAT Ver. 1.1
032 % revised ms 14.06.2000
033 % - updated header info
034 % - changed name to wgevlike (from gevll)
035 % revised pab 01.11.2000
036 % - added cov (from wgevfit)
037 
038 error(nargchk(2,2,nargin))
039 
040 sample = sample(:); % make sure it is a vector
041 N=length(sample);
042 y = -log(1-p(1)*(sample-p(3))/p(2))/p(1);
043 ll=N*log(p(2))+(1-p(1))*sum(y)+sum(exp(-y));
044 
045 if nargout>1,
046   % Calculate the covariance matrix by inverting the observed information. 
047     shape = p(1);
048     scale = p(2);
049     location = p(3);
050     y = -log(1-shape*(sample-location)/scale)/shape;
051     P = N-sum(exp(-y));
052     Q = sum(exp((shape-1)*y)+(shape-1)*exp(shape*y));
053     R = N-sum(y.*(1-exp(-y)));
054     dLda = -(P+Q)/scale/shape;
055     dLdb = -Q/scale;
056     dLdk = -(R-(P+Q)/shape)/shape;
057     dPdb = -sum(exp((shape-1)*y))/scale;
058     dPda = sum(exp(-y))/scale/shape+dPdb/shape;
059     dQdb = sum(exp((2*shape-1)*y));
060     dQdb = (dQdb+shape*sum(exp(2*shape*y)))*(1-shape)/scale;
061     dQda = sum(exp((shape-1)*y));
062     dQda = -(dQda+shape*sum(shape*y))*(1-shape)/scale/shape + dQdb/shape;    
063     dRda = N-sum(exp(shape*y))+sum(y.*exp(-y))-sum(exp(-y));
064     dRda = dRda+sum(exp((shape-1)*y))-sum(y.*exp((shape-1)*y));
065     dRda = -dRda/scale/shape;
066     dRdk = (sum(y)-sum(y.*exp(-y))+sum(y.*y.*exp(-y))-scale*dRda)/shape;
067     dRdb = (sum(exp(shape*y).*(1-exp(-y)+y.*exp(-y))))/scale;
068     d2Lda2 = -dLda/scale-(dPda+dQda)/scale/shape;
069     d2Ldab = -(dPdb+dQdb)/scale/shape;
070     d2Ldak = -(dRda+dLda+scale*d2Lda2)/shape;
071     d2Ldb2 = -dQdb/scale;
072     d2Ldbk = -(dRdb+scale*d2Ldab)/shape;
073     d2Ldk2 = -(dRdk+dLdk+scale*d2Ldak)/shape;
074     V = - [d2Ldk2, d2Ldak, d2Ldbk;
075            d2Ldak, d2Lda2, d2Ldab;
076            d2Ldbk, d2Ldab, d2Ldb2];
077     cov = real(inv(V));
078 end
079

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