gbd_momentfit.m

I developed an algorithm for using local ICA in denoising multidimensional data. It uses delay embed

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;

tmp0=(k^2*(-32+s)-9*s*(7+4*s)+6*k*(16+13*s));
tmp1=sqrt(s)*sqrt(tmp0);
tmp2=(2*k-3*(2+s))*tmp0;

a=3 * (1-k+s) * (tmp0-(k+3)*tmp1) / tmp2;
 
b=3 * (1-k+s) * (tmp0+(k+3)*tmp1) / tmp2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%