pearson_score.m

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

function [phi, phi2]=pearson_score(x,b,a,C);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PEARSON_SCORE - Calculates the score function and its derivative
% for a distribution from the Pearson system.
%
% Input:
%   x     = the realization matrix (sample values)
%   b     = the distribution parameters vector [b0 b1 b2]
%   a     = the distribution parameters vector [a0 a1]    
%   C     = the saturation limit (scalar); 
%           After saturation -C < phi(x) < C for all x
% Output:
%   phi   = the score function matrix corresponding to x  
%   phi2  = the derivative of the score function matrix corresponding to x  
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Juha Karvanen, Jan Eriksson, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.
%
% Last modified: 11.9.2001
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The classification is given in
%   Karvanen J. and Koivunen V.:
%   "Blind separation methods based on Pearson system and its extensions",
%   Signal Processing, 2002, to appear
%
% The estimates are given by the equations (4) and (5) in
%   Karvanen, J.,Eriksson, J., and Koivunen, V.:
%   "Pearson System Based Method for Blind Separation", 
%   Proceedings of Second International Workshop on
%   Independent Component Analysis and Blind Signal Separation,
%   Helsinki 2000, pp. 585--590
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if b(3)~=0 & b(3)~=-1,
delta=b(2)^2/(4*b(1)*b(3));
if a(2)==0,
ddelta=1;
else
ddelta=1+(delta-1)/(1+(4*b(3)/a(2)*(1+b(3)/a(2)))*delta);
end;
%ddelta is used to determine the type of distribution. 

if ddelta>=1, %unbound domain
C_1=(-a(2)-C*b(2)-sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
C_2=(-a(2)-C*b(2)+sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
C_3=(a(2)-C*b(2)-sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));
C_4=(a(2)-C*b(2)+sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));

if b(2)<0, %left bound domain
minscorelimit=max([C_1 C_2 C_3 C_4]);
if isreal(minscorelimit),
x=max(x,minscorelimit); 
end
else %right bound domain
maxscorelimit=min([C_1 C_2 C_3 C_4]);
if isreal(maxscorelimit),
x=min(x,maxscorelimit); 
end
end;

%elseif  (delta<=0 & delta>-1/(4*a(2)*b(3)*(1+a(2)*b(3)))),  %beta of first kind
elseif ddelta<0 | (ddelta==0 & b(3)>0),
%disp('beta')
C_1=(-a(2)-C*b(2)-sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
C_2=(-a(2)-C*b(2)+sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
C_3=(a(2)-C*b(2)-sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));
C_4=(a(2)-C*b(2)+sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));
CC=sort([C_1 C_2 C_3 C_4]);

minscorelimit=CC(2);
maxscorelimit=CC(3);

x=min(max(x,minscorelimit),maxscorelimit);

end;
elseif b(3)==0 & b(2)~=0,  %special case: Gamma distribution
C_1=(b(1)*C+a(1))/(a(2)-C*b(2))
C_2=(-b(1)*C+a(1))/(a(2)+C*b(2))

if b(2)<0, %left bound Gamma distribution
minscorelimit=max([C_1 C_2]);
if isreal(minscorelimit),
x=max(x,minscorelimit); 
end
else %right bound Gamma distribution
maxscorelimit=min([C_1 C_2]);
if isreal(maxscorelimit),
x=min(x,maxscorelimit); 
end
end;

end;

%the actual score function is calculated here:
denominator=b(1)+b(2)*x+b(3)*x.^2;
phi=-(a(2)*x-a(1))./denominator;
phi2=-(b(1)*a(2)+b(2)*a(1)+2*a(1)*b(3)*x-a(2)*b(3)*x.^2)./(denominator.^2);

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