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. See also weib2dpdf
Parameter estimates for 2D Weibull data. |
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 % 021 % See also weib2dpdf 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
Comments or corrections to the WAFO group