Home > wafo > wstats > mdist2dfit.m

mdist2dfit

PURPOSE ^

Parameter estimates for MDIST2D data.

SYNOPSIS ^

phato=mdist2dfit(V,H,dist,alpha,method)

DESCRIPTION ^

  MDIST2DFIT Parameter estimates for MDIST2D data. 
  
   CALL: Phat = mdist2dfit(x1,x2,dist,alpha,method) 
  
    phat = structure array containing 
           x      = cellarray of distribution parameters 
                    x{1} = Phat1 marginal parameters of x1 
                    x{2} = Phat2 marginal parameters of x2 
                    x{3} = Psi the interaction parameter between x1 and x2. 
           CI     = Confidence interval of distribution parameters 
           dist   = as explained below 
           alpha  = ------||---------- 
           method = ------||---------- 
           note   = Memorandum string 
           date   = date and time of creation  
  
  x1,x2  = input data 
  dist   = list of distributions used in the fitting of X1 and X2,  
           respectively. Options are 'tgumbel', 'gumbel',  
           'lognormal','rayleigh','weibull','gamma'. 
  alpha  = confidence level constant (If not given CI will not be computed) 
  method = fitting method used 'SML' or 'MML' (Simultainous or Marginal  
           Maximum Likelihood method) (default 'SML') 
  
  Example: 
   phat.x={1 2  2 }; 
   phat.dist={'rayl','rayl'}; 
   [R1,R2] = mdist2drnd(1000,phat); 
   Phat2 = mdist2dfit(R1,R2,{'weibull','rayleigh'},.05,'SML') 
   
  See also  mdist2drnd,  mdist2dpdf, mdist2dcdf, mdist2dprb

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function  phato=mdist2dfit(V,H,dist,alpha,method) 
002 % MDIST2DFIT Parameter estimates for MDIST2D data. 
003 % 
004 %  CALL: Phat = mdist2dfit(x1,x2,dist,alpha,method) 
005 % 
006 %   phat = structure array containing 
007 %          x      = cellarray of distribution parameters 
008 %                   x{1} = Phat1 marginal parameters of x1 
009 %                   x{2} = Phat2 marginal parameters of x2 
010 %                   x{3} = Psi the interaction parameter between x1 and x2. 
011 %          CI     = Confidence interval of distribution parameters 
012 %          dist   = as explained below 
013 %          alpha  = ------||---------- 
014 %          method = ------||---------- 
015 %          note   = Memorandum string 
016 %          date   = date and time of creation  
017 % 
018 % x1,x2  = input data 
019 % dist   = list of distributions used in the fitting of X1 and X2,  
020 %          respectively. Options are 'tgumbel', 'gumbel',  
021 %          'lognormal','rayleigh','weibull','gamma'. 
022 % alpha  = confidence level constant (If not given CI will not be computed) 
023 % method = fitting method used 'SML' or 'MML' (Simultainous or Marginal  
024 %          Maximum Likelihood method) (default 'SML') 
025 % 
026 % Example: 
027 %  phat.x={1 2  2 }; 
028 %  phat.dist={'rayl','rayl'}; 
029 %  [R1,R2] = mdist2drnd(1000,phat); 
030 %  Phat2 = mdist2dfit(R1,R2,{'weibull','rayleigh'},.05,'SML') 
031 %  
032 % See also  mdist2drnd,  mdist2dpdf, mdist2dcdf, mdist2dprb 
033  
034 % References: 
035 % Plackett, R. L. (1965) "A class of bivariate distributions." 
036 %                                J. Am. Stat. Assoc. 60. 516-22 
037 %      [1]  Michel K. Ochi, 
038 %       OCEAN TECHNOLOGY series 6 
039 %      "OCEAN WAVES, The stochastic approach", Cambridge 
040 %      1998 p. 133-134. 
041  
042  
043 %  tested on: matlab 5.2 
044 % history 
045 % Revised pab nov2004 
046 %  -replaced fmins with fminsearch   
047 % revised pab Dec2003 
048 % revised pab 8.11.1999 
049 %  - updated header info 
050 %  - changed phat from vector to structure 
051 % Per A. Brodtkorb 29.01.1999 
052  
053  
054 if (nargin< 3) 
055   error('To few input arguments') 
056 end 
057   CDIST=lower(dist{1}); 
058   UDIST=lower(dist{2}); 
059  
060 if (nargin<5)|isempty(method) 
061   method='SML'; % options MML, MOM or SML (simultaneous maximum likelihood) 
062  ;% marginal maximum likelihood  
063 end 
064 if (nargin < 4)|isempty(alpha) 
065     alpha = []; 
066 else 
067   p_int = [alpha/2; 1-alpha/2]; 
068 end 
069 V = V(:); 
070 H = H(:); 
071 [n, m] = size(V); 
072 [n2, m2] = size(H); 
073 if n~=n2 
074   error('V and H  must have equal size') 
075 end 
076  
077 [ ph, phci]=distfit(H,UDIST,method,alpha); 
078 [ pv,  pvci]=distfit(V,CDIST,method,alpha); 
079  
080 rho= corrcoef(V,H); 
081 rho=rho(2,1); 
082 psi=findPsi(rho);   
083 phat=[pv ph psi ]; 
084 phatCI=[pvci phci]; 
085 switch lower(method(1:3)), 
086   case 'sml', %simultaneous MLE 
087    pinit=phat; 
088    mvrs=version;ix=find(mvrs=='.'); 
089    if str2num(mvrs(1:ix(2)-1))>5.2, 
090      phat = fminsearch('mdist2dlike',pinit,[],V,H,dist); 
091    else 
092      phat = fmins('mdist2dlike',pinit,[],[],V,H,dist); 
093    end 
094   case 'mml',% marginal MLE is already found 
095       
096 end 
097  
098 if ~isempty(alpha) %  
099    [logL,cov]=mdist2dlike(phat,V,H,dist); 
100    sa = diag(cov)'; 
101    phatCI = wnorminv(repmat(p_int,1,length(phat)),[phat; phat],[sa;sa]); 
102  end 
103  phato.x=cell(2,1); 
104 phato.CI=cell(1,1); 
105 n1=length(pv); 
106 [phato.x{1}]=phat(1:n1); 
107 [phato.x{2}]=phat(n1+1:end-1); 
108 [phato.x{3}]=phat(end); 
109  
110 [phato.CI{1}]=phatCI; 
111 phato.dist=dist; 
112 phato.alpha=alpha; 
113 phato.method=method; 
114 phato.note=[]; 
115 phato.date=datestr(now); 
116  
117 function k=findPsi(rho) 
118   % finds psi by a simple iteration scheme  
119   rtol=1e-6; %relativ tolerance 
120   switch rho 
121   case 0,  k=1; return 
122   case 1, k=inf;return 
123   case -1, k=0; return 
124   otherwise,  
125     if rho<0, 
126       kold=0.5; 
127     else 
128      kold=4; 
129     end 
130   end 
131   kold2=kold-0.001; 
132   rho2=getrho(kold2); 
133   for ix=1:500, 
134     rho1=getrho(kold); 
135     k=kold-(rho1-rho).*(kold-kold2)./( rho1-rho2); 
136    %    disp(['k=' num2str(k) ]), pause 
137     if abs(kold-k)<abs(rtol.*k),break 
138     elseif k<=0, k=linspace(0,20,1000); 
139       [tmp I]=min(abs(getrho(k)-rho)); %linear search 
140       k=k(I); 
141       break; 
142     end 
143     kold=k; 
144     if ix>=500, disp('could not find k'),end 
145   end 
146   return 
147  
148 function r=getrho(psi) 
149 r=zeros(size(psi)); 
150 k1=find(psi==0); 
151 if any(k1), 
152  r(k)=-1; 
153 end 
154 k3=find((psi<1.*psi>0)|psi>1); 
155 if any(k3), 
156  r(k3)=(psi(k3).^2-1-2.*psi(k3).*log(psi(k3)))./((psi(k3)-1).^2) ; 
157 end 
158 return 
159  
160 function [pvhat,pvci]=distfit(tmp,dist,method,alpha) 
161    
162   if strcmp(lower(method(2:3)),'ml'), 
163     switch lower(dist(1:2)), 
164       case 'ra', [pvhat tmp2] = wraylfit(tmp,0) ; 
165       case 'tg', [pvhat tmp2] = wgumbfit(tmp,0); 
166       case 'gu', [pvhat tmp2] = wgumbfit(tmp,0); 
167       case 'lo', [pvhat tmp2] = wnormfit(log(tmp),0); 
168       case 'ga', [pvhat tmp2] = wgamfit(tmp,0); 
169       case 'we', [pvhat tmp2] = wweibfit(tmp,0); 
170       otherwise, error('Unknown distribution') 
171     end 
172     if prod(size(tmp2))~=length(tmp2); 
173       sa=diag(tmp2)'; 
174     else 
175       sa = tmp2(:)'; 
176     end 
177     p_int=[alpha/2; 1-alpha/2]; 
178     pvci = wnorminv(repmat(p_int,1,length(pvhat)),[pvhat; pvhat],[sa;sa]); 
179   else  % MOM fit 
180     pvci=NaN; 
181     error('MOM is not implemented for any distribution')     
182     switch dist(1:2) 
183       case {'tg','gu'} ,   
184       case 'we',  
185       case 'lo',  
186       case 'ga', error('MOM is not implemented for Gamma distribution')      
187       case 'ra', error('MOM is not implemented for Rayleigh distribution')      
188       otherwise , error('Unknown distribution') 
189     end 
190   end 
191   return 
192  
193

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