Home > wafo > wstats > winvgfit.m

winvgfit

PURPOSE ^

Parameter estimates for Inverse Gaussian data.

SYNOPSIS ^

[phat, var,ciL,ciU] = winvgfit(data,plotflag)

DESCRIPTION ^

 WINVGFIT Parameter estimates for Inverse Gaussian data.
 
  CALL:  [phat var] = winvgfit(data, plotflag)
 
    phat  = [m, l] = maximum likelihood estimate of the parameters of
            the distribution (see winvgpdf)
    var   = estimated asymptotic variance of phat (cov(m,l)=0)
    data  = data matrix
  plotflag = 0, do not plot
           > 0, plot the empiricial distribution function and the
                estimated cdf (see empdistr for options)(default)
 
  Example:
    R=winvgrnd(2,2,100,2);
    [phat, var]=winvgfit(R)
 
  See also  winvgcdf, empdistr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [phat, var,ciL,ciU] = winvgfit(data,plotflag)
002 %WINVGFIT Parameter estimates for Inverse Gaussian data.
003 %
004 % CALL:  [phat var] = winvgfit(data, plotflag)
005 %
006 %   phat  = [m, l] = maximum likelihood estimate of the parameters of
007 %           the distribution (see winvgpdf)
008 %   var   = estimated asymptotic variance of phat (cov(m,l)=0)
009 %   data  = data matrix
010 % plotflag = 0, do not plot
011 %          > 0, plot the empiricial distribution function and the
012 %               estimated cdf (see empdistr for options)(default)
013 %
014 % Example:
015 %   R=winvgrnd(2,2,100,2);
016 %   [phat, var]=winvgfit(R)
017 %
018 % See also  winvgcdf, empdistr
019 
020 % Reference: Chhikara & Folks, "The Inverse Gaussian Distribution", p. 53
021 
022 
023 %tested on: matlab 5.x
024 % History:
025 % revised pab 24.10.2000
026 % - added  nargchk
027 % - fixed some bugs when data is a matrix 
028 % - cov changed to var = variance since cov(m,l)=0
029 % added ciU, and ciL which is the 95% confidence interval  upper & lower
030 %     limits, respectively, for phat.
031 % added ms 14.08.2000
032 
033 error(nargchk(1,2,nargin))
034 if nargin<2|isempty(plotflag),  plotflag=1; end
035 sz = size(data);
036 Nsz=length(sz);
037 dim = min(find(sz~=1));  %1st non-singleton dimension
038 % make sure dim=1 is the first non-singleton dimension
039 if isempty(dim) | dim ~= 1, 
040   order = [dim 1:dim-1 dim+1:Nsz];
041   data  = permute(data,order);
042   sz    = size(data);
043 end
044 m = prod(sz(2:end));
045 n0 =sz(1);
046 
047 mhat=mean(data);
048 lhat=1./(mean(data(:,:).^(-1)-mhat(ones(n0,1),:).^(-1)));
049 phat=[mhat(:),lhat(:)];
050 
051 
052 if nargout>1
053   %nl/lhat is chi2(n-1)
054   n=max(6,n0);
055   var=[mhat(:).^3./lhat(:), lhat(:).^2*(n^3/(n-3)/(n-5)-n^3/(n-3)^2)]/n; 
056 end
057 
058 if nargout>2, % nl/lhat ~ chi2(n-1)
059   alpha2=ones(1,2)*0.05/2;
060   ciL = wnorminv(alpha2(ones(m,1),:),phat,var);
061   ciU = wnorminv(1-alpha2(ones(m,1),:),phat,var);
062 end
063 
064 if plotflag 
065   sd=sort(data);
066   empdistr(sd(:,1),[sd(:,1) winvgcdf(sd(:,1),mhat(1),lhat(1))],plotflag), hold on
067   for ix=2:m, empdistr(sd(:,ix),[sd(:,ix) winvgcdf(sd(:,ix),mhat(ix),lhat(ix))],plotflag),end 
068   hold off
069   title([deblank(['Empirical and Inverse Gaussian estimated cdf'])])
070 end
071 
072 
073 
074 
075 
076 
077 
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