Home > wafo > wstats > wraylfit.m

wraylfit

PURPOSE ^

Parameter estimates for Rayleigh data.

SYNOPSIS ^

[phat, var,pci] = wraylfit(data,plotflag)

DESCRIPTION ^

 WRAYLFIT Parameter estimates for Rayleigh data.
 
  CALL:  [bhat var] = wraylfit(data, plotflag)
 
    bhat  = maximum likelihood estimate of the parameter of
            the distribution (see wraylcdf)
    var   = estimated asymptotic variance of bhat
    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=wraylrnd(2,100,2);
    [bhat var]=wraylfit(R)
 
  See also  wraylcdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [phat, var,pci] = wraylfit(data,plotflag)
002 %WRAYLFIT Parameter estimates for Rayleigh data.
003 %
004 % CALL:  [bhat var] = wraylfit(data, plotflag)
005 %
006 %   bhat  = maximum likelihood estimate of the parameter of
007 %           the distribution (see wraylcdf)
008 %   var   = estimated asymptotic variance of bhat
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=wraylrnd(2,100,2);
016 %   [bhat var]=wraylfit(R)
017 %
018 % See also  wraylcdf
019 
020 % Reference: Cohen & Whittle, (1988) "Parameter Estimation in Reliability
021 % and Life Span Models", p. 181 ff, Marcel Dekker.
022 
023 %tested on: matlab 5.x
024 % History:
025 %  by Per A. Brodtkorb 17.10.98
026 % revised ms 15.06.2000
027 % - updated header info
028 % - changed name to wraylfit (from raylfit)
029 % revised ms 11.08.2000
030 % - changed to standard *fit form
031 % -revised pab 24.10.2000
032 %  - replaced gamma with gammaln -> more robust
033 %  - added nargchk
034 
035 error(nargchk(1,2,nargin))
036 if nargin<2|isempty(plotflag),  plotflag=1; end
037 sz = size(data);
038 Nsz=length(sz);
039 dim = min(find(sz~=1));  %1st non-singleton dimension
040 % make sure dim=1 is the first non-singleton dimension
041 if isempty(dim) | dim ~= 1, 
042   order = [dim 1:dim-1 dim+1:Nsz];
043   data  = permute(data,order);
044   sz    = size(data);
045 end
046 m = prod(sz(2:end));
047 n =sz(1);
048 
049 phat=sqrt(sum(data.^2)/n/2); % MLE
050 
051 if nargout > 1 
052   if (n<200),
053     var=phat.^2.*(n-(exp(gammaln(1/2+n)-gammaln(n))).^2)/n; 
054   else    %n too large for evaluating gamma
055     var=phat.^2.*0.25/n; 
056   end
057 end
058 
059 
060 if nargout>2,
061   alpha2=0.05/2;
062   pci = [wnorminv(alpha2,phat,var);wnorminv(1-alpha2,phat,var)];
063 end
064 
065 if plotflag 
066   sd=sort(data);
067   empdistr(sd(:,1),[sd(:,1),wraylcdf(sd(:,1),phat(1))],plotflag ), hold on
068   for ix=2:m,empdistr(sd(:,ix),[sd(:,ix),wraylcdf(sd(:,ix),phat(ix))],plotflag),end
069    hold off
070   
071   title([deblank(['Empirical and Rayleigh estimated cdf'])])
072 end
073 
074 
075 
076

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