Home > wafo > wstats > weib2dfit.m

weib2dfit

PURPOSE ^

Parameter estimates for 2D Weibull data.

SYNOPSIS ^

[phat, cov,pci]=weib2dfit(data1,data2,method,given,gparam,alpha)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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