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
Computes and plots the empirical CDF | |
Inverse Gaussian cumulative distribution function | |
Inverse of the Normal distribution function | |
Remove trailing blanks. | |
Display message and abort function. | |
Hold current graph. | |
Average or mean value. | |
Permute array dimensions. | |
Graph title. |
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
Comments or corrections to the WAFO group