WEIB2DFIT Parameter estimates for 2D Weibull data. CALL: [phat,cov,pci] = weib2dfit(data1,data2,method) phat = maximum likelihood estimates of the parameters of the distribution cov = estimated covariance of phat pci = 95% confidense intervals for phat data = vectors of data points method = 'SMLE' : Simultanous Maximum Likelihood Estimate (Default) 'MMLE' : Marginal Maximum Likelihood Estimate Example: R = wweibrnd(1,2,10000,2); phat = weib2dfit(R(:,1),R(:,2),'smle'); See also wweibfit, weib2dlike, wnorminv, hypgf, corrcoef
Hypergeometric function F(a,b,c,x) | |
2D Weibull log-likelihood function. | |
Inverse of the Normal distribution function | |
Parameter estimates for Weibull data. | |
Correlation coefficients. | |
Complete elliptic integral. | |
Display message and abort function. | |
Multidimensional unconstrained nonlinear minimization (Nelder-Mead). | |
Gamma function. | |
Convert string to lowercase. | |
Convert string matrix to numeric array. | |
MATLAB version number. |
001 function [phat, cov,pci]=weib2dfit(data1,data2,method,given,gparam,alpha) 002 %WEIB2DFIT Parameter estimates for 2D Weibull data. 003 % 004 % CALL: [phat,cov,pci] = weib2dfit(data1,data2,method) 005 % 006 % phat = maximum likelihood estimates of the parameters of the distribution 007 % cov = estimated covariance of phat 008 % pci = 95% confidense intervals for phat 009 % data = vectors of data points 010 % method = 'SMLE' : Simultanous Maximum Likelihood Estimate (Default) 011 % 'MMLE' : Marginal Maximum Likelihood Estimate 012 % Example: 013 % R = wweibrnd(1,2,10000,2); 014 % phat = weib2dfit(R(:,1),R(:,2),'smle'); 015 % 016 % See also wweibfit, weib2dlike, wnorminv, hypgf, corrcoef 017 018 % tested on: matlab 5.2 019 %History: 020 % revised pab nov 2004 021 % fmins replaced with fminsearch 022 % revised pab 02.11.2000 023 % by Per A. Brodtkorb 14.11.98 024 025 % alpha = confidence level (default 0.05 corresponding to 95% CI) 026 % g = indices to fixed parameters not estimated (Default []) 027 % gphat = values for the fixed parameters not estimated (Default []) 028 029 % Example: 030 % R = wweibrnd(1,2,10000,2); 031 % phat0 = [2 2 0]; % set the B1 B2 and C12 respectively a priory 032 % given = [2 4 5]; % = indices to the parameters phat0 set a priory 033 % phat = weib2dfit(R(:,1),R(:,2),'smle',given,phat0); % estimate A1 and A2 034 035 036 error(nargchk(2,6,nargin)) 037 if (nargin < 3 |isempty(method)), method = 'SMLE'; end 038 if (nargin < 5 |isempty(gparam)), gparam = []; end 039 if (nargin < 4 |isempty(given)), given = []; end 040 if (nargin < 6)|isempty(alpha), alpha = 0.05; end 041 042 043 data1 = data1(:); data2 = data2(:); 044 n = length(data1); n2 = length(data2); 045 if n~=n2, error('data1 and data2 must have equal size'),end 046 047 rho = corrcoef(data1,data2); 048 rho = rho(2,1); 049 pinit = [wweibfit(data1) wweibfit(data2) sqrt(abs(rho))*sign(rho) ]; 050 %pinit(given) = gparam; 051 %pinit(5)=findk(rho,pinit); 052 053 switch lower(method(1:4)), 054 case 'smle', %simultanous MLE 055 pinit(given)=[]; 056 mvrs=version;ix=find(mvrs=='.'); 057 if str2num(mvrs(1:ix(2)-1))>5.2, 058 phat = fminsearch('weib2dlike',pinit,[],data1,data2,given,gparam); 059 else 060 phat = fmins('weib2dlike',pinit,[],[],data1,data2,given,gparam); 061 end 062 063 case 'mmle',% marginal MLE 064 phat=pinit; 065 phat(given)=gparam; 066 phat(5)=findk(rho,phat); 067 otherwise 068 phat = pinit; 069 end 070 071 if nargout > 1, 072 p_int = [alpha/2; 1-alpha/2]; 073 [logL,cov] = weib2dlike(phat,data1,data2,given,gparam); 074 sa = diag(cov).'; 075 pci = wnorminv(repmat(p_int,1,5-length(given)),[phat; phat],[sa;sa]); 076 end 077 078 function k=findk(rho,sparam) 079 % finds k by a simple iteration scheme 080 rtol=1e-6; %relativ tolerance 081 kold=sqrt(abs(rho))*sign(rho); 082 kold2=kold-0.001; 083 rho2=getrho(kold2,sparam); 084 for ix=1:500, 085 rho1=getrho(kold,sparam); 086 k=kold-(rho1-rho).*(kold-kold2)./( rho1-rho2); 087 % disp(['k=' num2str(k) ]), pause 088 if abs(kold-k)<abs(rtol.*k),break 089 elseif abs(k)>1, tmp=sqrt(abs(rho))*sign(rho):0.000001:kold; 090 [tmp2 I]=min(abs(getrho(tmp,sparam)-rho)); %linear search 091 k=tmp(I); 092 break, 093 end 094 rho2=rho1; kold2=kold; 095 kold=k; 096 if ix>=500, disp('could not find k'),end 097 end 098 return 099 100 function y=getrho(k,sparam) 101 % returns the correlationcoefficient based on: 102 a1=sparam(1);b1=sparam(2);a2=sparam(3);b2=sparam(4); 103 % and k 104 if 0, % calculating correlation correct if rayleigh 105 [K,E] = ellipke(k.^2); %complete elliptic integral of first and 106 %second kind 107 y= (E-.5*(1-k.^2).*K-pi/4)./(1-pi/4); 108 else % alternatively and possibly slower 109 qx=1./b1;qy=1./b2; 110 y= (hypgf(-qx,-qy,1,k.^2)-1).* gamma(qx).*gamma(qy) ./ .... 111 ( sqrt(2.*b1.*gamma(2.*qx) - gamma(qx).^2) .* ... 112 sqrt(2.*b2.*gamma(2.*qy) - gamma(qy).^2) ); 113 114 115 end 116 return 117
Comments or corrections to the WAFO group