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
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
Comments or corrections to the WAFO group