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)

## CROSS-REFERENCE INFORMATION

This function calls:
 empdistr Computes and plots the empirical CDF winvgcdf Inverse Gaussian cumulative distribution function wnorminv Inverse of the Normal distribution function deblank Remove trailing blanks. error Display message and abort function. hold Hold current graph. mean Average or mean value. permute Permute array dimensions. title Graph title.
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 %
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
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.
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