Home > wafo > wstats > dist2dfit.m

dist2dfit

PURPOSE ^

Parameter estimates for DIST2D data.

SYNOPSIS ^

[phat] =dist2dfit(V,H,dist,res,method,monitor,chat0)

DESCRIPTION ^

 DIST2DFIT Parameter estimates for DIST2D data. 
  
  CALL:  phat = dist2dfit(x1,x2,dist,[res,noverlap C],method,monitor); 
  
         phat = structure array containing 
                x      = cellarray of distribution parameters 
                CI     = Confidence interval of distribution parameters 
                dist   = as explained below 
                method = ------||---------- 
                note   = string 
                date   = date and time of creation  
        x1,x2 = input data 
         dist = list of distributions used in the fitting of X1 given X2  
                and X2, respectively. Options are 'tgumbel', 'gumbel',  
                'lognormal','rayleigh','trayleigh','weibull','tweibull', 
                'gamma','ggamma'. 
          res = resolution in the conditional fitting  (default range(x2)/12) 
     noverlap = 0,1,2,3..., nr. of groups above and below to overlap in 
                the conditional fitting (default 0) 
            C = defines the location parameter by the following: 
                Cv = min({X1|X2=x}-C,0) 
                if no value is given then Cv=0.  
       method = list containing method used for fitting of  
                X1 given X2 and X2, respectively. 
                Options 'MLE' or 'MOM'  (default {'MLE',MLE'}) 
      monitor = 1 if monitor the fitting  (default=0)  
  
   DIST2DFIT fits the following distribution: F{x1|X2=x2}*F{x2}.  
   The parameter(s) of the unconditional distribution of X2, 
   F{x2}, is  returned in phat.x{2}. The parameters of the conditional 
   distribution of X1 given X2 are returned in phat.x{1}. The first column 
   in PHAT.X{1} contains the X2 values the parameters in column 2 and 3 are 
   conditioned on.  PHAT.CI{1} and PHAT.CI{2} gives 95% confidence 
   interval if MLE are selected for METHOD{1} and METHOD{2}, respectively. 
  
  Example: 
    R = wraylrnd(2,1000,2); x1 = linspace(0,10)'; 
    phat0=struct('x',[]); phat0.x={[x1,2*ones(size(x1))] 2 }; 
    phat0.dist={'rayl','rayl'}; 
    phat = dist2dfit(R(:,1),R(:,2),{'rayl','rayl'});  
    sphat = dist2dsmfun2(phat,x1,0); % smooth the parameters 
    dist2dparamplot(phat,sphat);  
    figure(2) % compare fitted model with theoretical 
    f = dist2dpdf2(x1,x1,phat);  
    fs = dist2dpdf2(x1,x1,sphat);  
    fe = dist2dpdf2(x1,x1,phat0);  
    pdfplot(fs); hold on, 
    pdfplot(fe,'k--') 
    plot(R(:,1),R(:,2),'.'),hold off 
  
   See also    dist2drnd,  dist2dpdf dist2dcdf dist2dprb

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function  [phat] =dist2dfit(V,H,dist,res,method,monitor,chat0) 
0002 %DIST2DFIT Parameter estimates for DIST2D data. 
0003 % 
0004 % CALL:  phat = dist2dfit(x1,x2,dist,[res,noverlap C],method,monitor); 
0005 % 
0006 %        phat = structure array containing 
0007 %               x      = cellarray of distribution parameters 
0008 %               CI     = Confidence interval of distribution parameters 
0009 %               dist   = as explained below 
0010 %               method = ------||---------- 
0011 %               note   = string 
0012 %               date   = date and time of creation  
0013 %       x1,x2 = input data 
0014 %        dist = list of distributions used in the fitting of X1 given X2  
0015 %               and X2, respectively. Options are 'tgumbel', 'gumbel',  
0016 %               'lognormal','rayleigh','trayleigh','weibull','tweibull', 
0017 %               'gamma','ggamma'. 
0018 %         res = resolution in the conditional fitting  (default range(x2)/12) 
0019 %    noverlap = 0,1,2,3..., nr. of groups above and below to overlap in 
0020 %               the conditional fitting (default 0) 
0021 %           C = defines the location parameter by the following: 
0022 %               Cv = min({X1|X2=x}-C,0) 
0023 %               if no value is given then Cv=0.  
0024 %      method = list containing method used for fitting of  
0025 %               X1 given X2 and X2, respectively. 
0026 %               Options 'MLE' or 'MOM'  (default {'MLE',MLE'}) 
0027 %     monitor = 1 if monitor the fitting  (default=0)  
0028 % 
0029 %  DIST2DFIT fits the following distribution: F{x1|X2=x2}*F{x2}.  
0030 %  The parameter(s) of the unconditional distribution of X2, 
0031 %  F{x2}, is  returned in phat.x{2}. The parameters of the conditional 
0032 %  distribution of X1 given X2 are returned in phat.x{1}. The first column 
0033 %  in PHAT.X{1} contains the X2 values the parameters in column 2 and 3 are 
0034 %  conditioned on.  PHAT.CI{1} and PHAT.CI{2} gives 95% confidence 
0035 %  interval if MLE are selected for METHOD{1} and METHOD{2}, respectively. 
0036 % 
0037 % Example: 
0038 %   R = wraylrnd(2,1000,2); x1 = linspace(0,10)'; 
0039 %   phat0=struct('x',[]); phat0.x={[x1,2*ones(size(x1))] 2 }; 
0040 %   phat0.dist={'rayl','rayl'}; 
0041 %   phat = dist2dfit(R(:,1),R(:,2),{'rayl','rayl'});  
0042 %   sphat = dist2dsmfun2(phat,x1,0); % smooth the parameters 
0043 %   dist2dparamplot(phat,sphat);  
0044 %   figure(2) % compare fitted model with theoretical 
0045 %   f = dist2dpdf2(x1,x1,phat);  
0046 %   fs = dist2dpdf2(x1,x1,sphat);  
0047 %   fe = dist2dpdf2(x1,x1,phat0);  
0048 %   pdfplot(fs); hold on, 
0049 %   pdfplot(fe,'k--') 
0050 %   plot(R(:,1),R(:,2),'.'),hold off 
0051 % 
0052 %  See also    dist2drnd,  dist2dpdf dist2dcdf dist2dprb 
0053  
0054  
0055 % tested on: matlab 5.2 
0056 % history: 
0057 % revised pab Sept 2005 
0058 % -fixed abug: [a{1}] = 1 is not allowed anymore!! replaced with a{1} = 1; 
0059 % revised pab Jul2004 
0060 % Added secret option of C0 to wggamfit  
0061 % revised pab 03.12.2000 
0062 % -added truncated weibull and rayleigh 
0063 % revised pab 12.11.2000 
0064 %  - added wggampdf option 
0065 %  - For monitor> 0:  removed normplot ...etc and replaced it with a call 
0066 %     to empdistr instead 
0067 % Per A. Brodtkorb 20.10.1998 
0068  
0069 error(nargchk(3,7,nargin)) 
0070 ptime     = 1; %pause length if monitor=1 
0071 printflag = 0;%print if monitor=1 
0072  
0073 [HI I] = sort(H(:));%sorting H 
0074 VI     = V(I); 
0075  
0076 if (nargin< 3)|isempty(dist), 
0077   dist={'tgumbel','rayleigh'}; 
0078 else 
0079   nd=length(dist); 
0080   if nd==1, 
0081   dist={dist{1},'rayleigh'}; 
0082   end 
0083 end 
0084  
0085 if (nargin<5)|isempty(method) 
0086   method={'MLE','MLE'}; % options MLE or MOM 
0087 else 
0088   nd=length(method); 
0089   if nd==1, 
0090     method={method{1}, 'MLE'}; 
0091   end 
0092 end 
0093  
0094 if (nargin<6)|isempty(monitor) 
0095   monitor=0; %if 1 monitor the fitting  
0096 end 
0097 %monitor=1; 
0098  
0099  
0100 noverlap=0;C=[]; 
0101 if nargin<4|isempty(res) 
0102   res=range(H(:))/12; % resolution 
0103 else 
0104   if length(res)>1 
0105     noverlap=res(2); 
0106   end 
0107   if length(res)>2 
0108     C=res(3); 
0109   end 
0110   res=res(1); 
0111 end 
0112  
0113 %N=length(HI); 
0114  
0115 grp=floor(HI/res)+1; % dividing the data into groups 
0116  
0117 if monitor 
0118   [len,bin] = bincount(grp); 
0119   utp = 1/(2*max(len)); 
0120   ax = [0 max(VI),utp,1-utp]; 
0121   name1=inputname(1); 
0122   name2=inputname(2); 
0123   if isempty(name1), name1 = 'x1'; end 
0124   if isempty(name2), name2 = 'x2'; end 
0125    
0126 end 
0127 [ phhat, phci]=distfit(HI,dist{2},method{2},monitor);%unconditional fitting 
0128 if monitor 
0129   xlabel(name2) 
0130   pause(ptime) 
0131   if printflag, print -Pmhlaser ; end   
0132 end 
0133  
0134 Ngrp=grp(end); % # groups 
0135  
0136 if strcmpi(dist{1}(1:2),'ra'), 
0137   npar=1; 
0138 elseif strcmpi(dist{1}(1:2),'gg')|strcmpi(dist{1}(1:2),'tw'), 
0139   npar = 3; 
0140 else 
0141   npar=2;   
0142 end 
0143 %C = input('Specify location parameter (default none)');  
0144 if isempty(C),   
0145   pvhat=zeros(Ngrp,1+npar); 
0146 else 
0147   Cvold=0; 
0148   pvhat=zeros(Ngrp,2+npar); 
0149 end 
0150  
0151 pvci=zeros(Ngrp,2*npar); 
0152  
0153 pvhat(:,1)=linspace(res/2, (Ngrp-0.5)*res, Ngrp)'; 
0154 if nargin<7 
0155   chat0 = []; 
0156 end 
0157  
0158 Nmin = min(max(6,Ngrp),25); % Minimum number of data in groups 
0159 m = zeros(Ngrp,1); 
0160 v = m; 
0161 Ni = m; 
0162 for ix=1:Ngrp, 
0163   ind  = (grp==ix);   
0164   tmp  = VI(ind); 
0165   Ni(ix) = length(tmp); 
0166   if length(tmp)>max(4,0),% if less than 4 observations in the group  
0167     m(ix)=mean(tmp); % mean of data in group ix 
0168     v(ix)=std(tmp).^2;  % variance of data in group ix 
0169   else  
0170     m(ix)=NaN; 
0171     v(ix)=NaN; 
0172   end 
0173   %if 1 | ((ix-1)*res>2) 
0174   %  chat0 = 1.5; 
0175   %end 
0176   switch class(chat0) 
0177    case 'inline' 
0178     chat00 = chat0(pvhat(ix,1)); 
0179    case 'struct' 
0180     chat00 = ppval(chat0,pvhat(ix,1)); 
0181    otherwise 
0182     chat00 = chat0;  
0183   end 
0184    
0185   if ix>=Ngrp-noverlap+1 | ix<=noverlap 
0186     for iy=1:min(min(noverlap,ix-1),Ngrp-ix) 
0187       kup = (grp==ix+iy); 
0188       if any(kup) 
0189     ind = (ind | (grp==ix-iy) | kup); 
0190       end 
0191     end  
0192   end              
0193    
0194   tmp=VI(ind);%find data in group number ix 
0195    
0196   if length(tmp)<Nmin,% if less than Nmin observations in the group  
0197     %grp(grp==ix)=ix-1; %merge the data with group below 
0198     %if (ix>3)&~isempty(tmp), grp(grp==ix-2)=ix-1;end % also merge the 
0199                                                      % data in grp ix-2 with grp ix-1 
0200     pvhat(ix,:)=NaN; 
0201     pvci(ix,:)=NaN; 
0202   else 
0203      
0204     if ~isempty(C) 
0205       %disp('Cv') 
0206       Cv=max(min(tmp)-C,0); 
0207        
0208       Cv=min(Cv,Cvold+C/30); 
0209       tmp=tmp-Cv; 
0210       Cvold=Cv; 
0211       pvhat(ix,end)=Cv; 
0212     end 
0213     % conditional fitting  
0214     [ pvhat(ix, 2:npar+1),... 
0215       pvci(ix,1:2*npar)]=distfit(tmp(:),dist{1},method{1},monitor,chat00); 
0216     if monitor 
0217       axis(ax) 
0218       xlabel([name1 ' |  '  name2 '=' num2str(pvhat(ix,1)) ]) 
0219       pause(ptime) 
0220       if printflag, print -Pmhlaser ; end   %print -Pmhlaser 
0221     end 
0222   end 
0223 end 
0224  
0225 %phat=struct('x',[],'CI',[],'dist',dist,'method',method); 
0226 phat.x=cell(2,1); 
0227 phat.CI=cell(2,1); 
0228 phat.x{1}=pvhat; 
0229 phat.x{2}=phhat; 
0230 phat.CI{1}=pvci; 
0231 phat.CI{2}=phci; 
0232 phat.dist=dist; 
0233 phat.method=method; 
0234 phat.note=[]; 
0235 phat.date=datestr(now); 
0236 phat.res=res; 
0237 phat.noverlap=noverlap; 
0238 phat.C=C; 
0239 phat.visual=[]; 
0240 phat.csm=[]; 
0241 phat.lin=[]; 
0242 phat.stats1{1} = pvhat(:,1); % x2 
0243 phat.stats1{2} = m; %conditional mean given x2 
0244 phat.stats1{3} = v; % conditional variance given x2 
0245 phat.stats1{4} = Ni;% conditional number of samples given x2 
0246 %pvhat(end,:)=[]; 
0247 %pvhat(1,:)=[]; 
0248 %pvci(end,:)=[]; 
0249 %pvci(1,:)=[]; 
0250 function [pvhat,pvci]=distfit(tmp,dist,method,monitor,chat) 
0251   if nargin>4 & ~isempty(chat) & isfinite(chat) 
0252     chat0 = chat; 
0253   else 
0254     chat0 = []; 
0255   end 
0256   if strncmpi(method,'mle',3), 
0257     switch lower(dist(1:2)), 
0258       case 'tr', [pvhat, tmp2] = wtraylfit(tmp,monitor) ; 
0259       case 'ra', [pvhat, tmp2] = wraylfit(tmp,monitor) ; 
0260       case 'tg', [pvhat, tmp2] = wgumbfit(tmp,monitor); 
0261       case 'gu', [pvhat, tmp2] = wgumbfit(tmp,monitor); 
0262       case 'lo', [pvhat, tmp2] = wnormfit(log(tmp),monitor); 
0263       case 'ga', [pvhat, tmp2] = wgamfit(tmp,monitor); 
0264      case 'gg',  [pvhat, tmp2] = wggamfit(tmp,monitor,chat0); 
0265     %if any(pvhat<0),pvhat(1:3) =NaN; end 
0266       case 'we', [pvhat, tmp2] = wweibfit(tmp,monitor); 
0267       case 'tw', [pvhat, tmp2] = wtweibfit(tmp,monitor); 
0268     %[tmp tmp2]=gumbfit(log(tmp)); 
0269     %pvhat=[exp(tmp(2)) 1/tmp(1) ] 
0270       otherwise, error('Unknown distribution') 
0271     end 
0272     if any(isnan(pvhat)) 
0273       pvci=[pvhat pvhat]; 
0274     else 
0275       tz=size(tmp2); 
0276       alpha=0.05; % 95% confidense interval 
0277       pint = [alpha/2; 1-alpha/2]; 
0278       if tz(1)==tz(2), 
0279     sa = diag(tmp2)'; 
0280       else 
0281     sa  = tmp2(:)'; 
0282       end 
0283       nt=length(sa); 
0284       %pvhat 
0285       pvci = wnorminv(repmat(pint,[1 nt]),[pvhat;pvhat],[sa;sa]); 
0286       pvci=pvci(:).'; 
0287     end 
0288     if 0 %monitor 
0289       switch lower(dist(1:2)), 
0290     case 'ra',wraylplot(tmp); 
0291     case 'gu',wgumbplot(tmp); 
0292     case 'tg',wgumbplot(tmp); 
0293     case 'lo',wnormplot(log(tmp)); 
0294     case 'ga',distplot(tmp,'wgampdf'); 
0295     case 'gg',distplot(tmp,'wggampdf'); 
0296     case 'we',wweibplot(tmp); 
0297     otherwise, error('Unknown distribution') 
0298       end     
0299     end 
0300      
0301   else  % MOM fit 
0302     pvci=NaN; 
0303     switch lower(dist(1:2)) 
0304      case {'tg','gu'} ,  pvhat =wgumbplot(tmp); 
0305        
0306      case 'we', pvhat = wweibplot(tmp); 
0307        
0308      case 'lo', pvhat = wnormplot(log(tmp)); 
0309        
0310      case 'ga', ma=mean(tmp);sa=var(tmp); 
0311        pvhat=[ma^2/sa sa/ma]; 
0312        
0313      case 'ra', error('MOM is not implemented for Rayleigh distribution')      
0314      otherwise , error('Unknown distribution') 
0315     end 
0316     pvci = repmat(nan,1,2*length(pvhat)); 
0317      
0318     if monitor 
0319       switch lower(dist(1:2)), 
0320     case {'gu' 'tg'}, wgumbplot(tmp); 
0321     case 'lo',  wnormplot(log(tmp)); 
0322     case 'ga', distplot(tmp,'wgampdf'); 
0323     case 'we',  wweibplot(tmp); 
0324     otherwise, error('Unknown distribution') 
0325       end 
0326     end 
0327   end 
0328   return 
0329  
0330  
0331  
0332  
0333  
0334

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