gbd_momentfit.m

PearsonICA程序

function beta=gbd_momentfit(alpha3,alpha4,samplemin,samplemax,samplen);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GBD_MOMENTFIT - Estimates the parameters of the Generalized Beta 
% Distribution (GBD) using the first four sample moments
% (the method of moments). Alternatively, the sample minimum and maximum 
% can be used instead of the sample mean and variance.
%
% Input
%   alpha3    = sample 3rd moment
%   alpha4    = sample 4th moment
%   samplemin = sample mean     or sample minimum
%   samplemax = sample variance or sample maximum
%   samplelen =                    sample length
%  
% Output
%   beta = the distribution parameters vector [beta(1) beta(2) beta(3) beta(4)]
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Jan Eriksson, Juha Karvanen, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.
%
% Last modified: 5.9.2000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<4,  disp(['Usege: gbd_momentfit(alpha3,alpha4,samplemin,', ...
		     'samplemax,samplen) or']);
              disp(['gbd_momentfit(alpha3,', ...
		     'alpha4,alpha1,alpha2)']); 
              error('Not enough input arguments.');
end

if or(alpha4<=1+alpha3^2,alpha4>=3+2*alpha3^2),
              warning(['The moment values are not covered by the '...
		       'GBD distribution.']);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Moment estimates for beta(3) and beta(4). These are obtained from
% the equations (15) and (16) in
%   Eriksson, J., Karvanen, J., and Koivunen, V.:
%   "Source Distribution Adaptive Maximum Likelihood Estimation of
%   ICA Model", Proceedings of Second International Workshop on
%   Independent Component Analysis and Blind Signal Separation,
%   Helsinki 2000, pp. 227--232
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s=alpha3^2;
k=alpha4;
a=-((3*(-1 + k - s)*(k^2*(-32 + s) + ...
      k*(96 + 78*s - sqrt(s)*sqrt(k^2*(-32 + s) - ...
           9*s*(7 + 4*s) + 6*k*(16 + 13*s))) - ...
      3*(21*s + 12*s^2 + sqrt(s)*sqrt(k^2*(-32 + s) - ...
           9*s*(7 + 4*s) + 6*k*(16 + 13*s)))))/ ...
    ((2*k - 3*(2 + s))*(k^2*(-32 + s) - 9*s*(7 + 4*s) + ...
    6*k*(16 + 13*s))));
 
 b=-(3*(-1 + k - s)*(k^2*(-32 + s) - 63*s - 36*s^2 + ...
     3*sqrt(s)*sqrt(k^2*(-32 + s) - 9*s*(7 + 4*s) + ...
        6*k*(16 + 13*s)) + ...
     k*(96 + 78*s + sqrt(s)*sqrt(k^2*(-32 + s) - ...
          9*s*(7 + 4*s) + 6*k*(16 + 13*s)))))/ ...
  (2*k^3*(-32 + s) + k^2*(384 + 246*s - 3*s^2) + ...
   27*s*(14 + 15*s + 4*s^2) - 18*k*(32 + 49*s + 17*s^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: A=beta(3)+1, B=beta(4)+1 correspond to the Statistics Toolbox
% values used for the ordinary Beta Distribution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if alpha3>0,
   beta(3)=a-1; 
   beta(4)=b-1; 
else
   beta(3)=b-1; 
   beta(4)=a-1; 
end;

if nargin==5
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Estimates for beta(1) and beta(2) based on sample minimum and
 % maximum. The PDF of any GBD distribution is nonzero on interval [beta(1),
 % beta(1)+beta(2)], therefore this "ad-hoc" estimation guarantees
 % that the corresponding pdf is non-zero for all realizations 
 % (this is not the case with the moment estimation, which may
 % lead to problems when calculating the score function...).
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 beta(1)=(samplen+1)*samplemin/samplen;
 beta(2)=(samplen+1)*samplemax/samplen-beta(1);
else
 alpha1=samplemin;
 alpha2=samplemax;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Moment estimates for beta(1) and beta(2). These are given by
 % equations (17) and (18) in
 %   Eriksson, J., Karvanen, J., and Koivunen, V.: 
 %   "Source Distribution Adaptive Maximum Likelihood Estimation of
 %   ICA Model", Proceedings of Second International Workshop on
 %   Independent Component Analysis and Blind Signal Separation,
 %   Helsinki 2000, pp. 227--232
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 beta(2)=(beta(3)+beta(4)+2)*...
   sqrt(alpha2*(beta(3)+beta(4)+3)/((beta(3)+1)*(beta(4)+1)));
 beta(1)=alpha1-beta(2)*(beta(3)+1)/(beta(3)+beta(4)+2);
end

% The end of the function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%