Home > wafo > wstats > wgamfit.m

wgamfit

PURPOSE ^

Parameter estimates for Gamma data.

SYNOPSIS ^

[phat, cov,pci]=wgamfit(data1, plotflag);

DESCRIPTION ^

 WGAMFIT Parameter estimates for Gamma data. 
  
  CALL: [phat, cov] = wgamfit(data, plotflag) 
  
      phat = [a,b] = the maximum likelihood estimates of the   
             parameters of the Gamma distribution 
             (see wgamcdf) given the data. 
      cov  = asymptotic covariance matrix of estimates  
      data = data vector 
  plotflag = 0, do not plot 
           > 0, plot the empiricial distribution function and the 
                estimated cdf (see empdistr for options)(default) 
            
  Example: 
    R = wgamrnd(5,1,1,100); 
    phat = wgamfit(R,2)  
  
  See also  wgamcdf, empdistr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [phat, cov,pci]=wgamfit(data1, plotflag); 
002 %WGAMFIT Parameter estimates for Gamma data. 
003 % 
004 % CALL: [phat, cov] = wgamfit(data, plotflag) 
005 % 
006 %     phat = [a,b] = the maximum likelihood estimates of the   
007 %            parameters of the Gamma distribution 
008 %            (see wgamcdf) given the data. 
009 %     cov  = asymptotic covariance matrix of estimates  
010 %     data = data vector 
011 % plotflag = 0, do not plot 
012 %          > 0, plot the empiricial distribution function and the 
013 %               estimated cdf (see empdistr for options)(default) 
014 %           
015 % Example: 
016 %   R = wgamrnd(5,1,1,100); 
017 %   phat = wgamfit(R,2)  
018 % 
019 % See also  wgamcdf, empdistr 
020  
021 % Reference: 
022 % Bowman & Shenton, "Properties of estimators for the Gamma distribution" 
023 % Marcel Dekker 
024  
025 % tested on: matlab 5.3 
026 % History: 
027 % revised pab 21.01.2004 
028 % revised pab 24.10.2000 
029 % - added check on fzero in order to run on matlab 5.2 
030 % - added pci 
031 % added ms 23.08.2000 
032  
033  
034 error(nargchk(1,2,nargin)) 
035 if any(data1<=0) 
036   error('data must be strictly positive!') 
037 end 
038 if nargin<2|isempty(plotflag),plotflag=1;end 
039  
040 data=data1(:); % make sure it is a vector 
041  
042 meanData    = mean(data); 
043 logData     = log(data); 
044 meanLogData = mean(logData); 
045 ahat0 = -(meanLogData-mean(data.*logData)/meanData)^(-1); 
046 if 1, 
047   G     = exp(meanLogData); 
048   start = 1./(2*log(meanData/G)); 
049   %start = ahat0 
050   mvrs=version; 
051   ix=find(mvrs=='.'); 
052   if str2num(mvrs(1:ix(2)-1))>5.2, 
053     ahat = fzero('wgamafit',start,optimset,meanData,G); 
054   else 
055     ahat = fzero('wgamafit',start,sqrt(eps),[],meanData,G);  
056   end     
057 else 
058   ahat = ahat0;   
059 end 
060 bhat = meanData/ahat; 
061  
062 phat=[ahat, bhat]; 
063  
064 if nargout>1, 
065   [LL, cov] = loglike(phat,data,'wgampdf'); 
066 end 
067 if nargout>2 
068   alpha2 = ones(1,2)*0.05/2; 
069   var = diag(cov).'; 
070   pci = wnorminv([alpha2;1-alpha2], [phat;phat],[var;var]); 
071 end 
072  
073 if plotflag  
074   sd = sort(data);  
075   empdistr(sd,[sd, wgamcdf(sd,ahat,bhat)],plotflag) 
076   title([deblank(['Empirical and Gamma estimated cdf'])]) 
077 end 
078

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