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)

CROSS-REFERENCE INFORMATION

This function calls:
 error Display message and abort function. inv Matrix inverse.
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 %
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