function I=quad_tsallis_fast(X,Y,xcov,ycov,NS) % function I=quad_tsallis_fast(X,Y,xcov,ycov) % % Quadratic Tsallis mutual information between X and Y. % Defines mutual information as the sum of the individual % entropies, minus the joint entropy. % % This can be written as I(X;Y) = H(X)+H(Y)-H(X)H(Y) - H(X,Y), % which simplifies to \int p(x,y)^2 - \int p(x)^2 \int p(y)^2. % % Uses a Nystrom approximation to speed up the computation. % Specify the number of points to use in the approximation % with the NS parameter. % % Assumes a Parzen density estimate for both variables, % with spherical, homoscedastic Gaussian kernels. % % Currently, this requires that you have the same % number of samples for X and Y, because it assumes that aligned % samples are samples from the joint of p(x,y) % % D. Wingate 10/25/2006 % [dimx, nx] = size(X); [dimy, ny] = size(Y); if ( nx ~= ny ) error( 'Number of samples must be equal!' ); return; end; xcov = 2*xcov; ycov = 2*ycov; norm_x = 1/sqrt((2*pi*xcov)^dimx); norm_y = 1/sqrt((2*pi*ycov)^dimy); nsq = 1/(nx*nx); % the joint distribution of X and Y. % we premultiply by the gaussian width here because % it may be different for X and Y, but our rbf4nn % function can't cope with that... XY = [ 1/sqrt(xcov) * X; 1/sqrt(ycov) * Y ]; % get up to NS linearly independent points %iX = get_lin_ind_pts( X, xcov, 2*xcov, NS ); %iY = get_lin_ind_pts( Y, ycov, 2*ycov, NS ); [iX,bX,bXinv] = get_lin_ind_pts_K( X, xcov, 0.001, NS ); [iY,bY,bYinv] = get_lin_ind_pts_K( Y, ycov, 0.001, NS ); iXY = unique( [iX iY] ); % the union of X and Y % the subsampled points sX = X( :, iX ); sY = Y( :, iY ); sXY = XY( :, iXY ); % entropy of X vX = sum( rbf4nn( sX, X, xcov ), 2 ); SX = norm_x * nsq * ( vX' * bXinv * vX ); % entropy of Y vY = sum( rbf4nn( sY, Y, ycov ), 2 ); SY = norm_y * nsq * ( vY' * bYinv * vY ); % joint entropy bXY = rbf4nn( sXY, sXY, 1 ); vXY = sum( rbf4nn( sXY, XY, 1 ), 2 ); SXY = (norm_x*norm_y) * nsq * (vXY' * pinv( bXY ) * vXY ); % Equivalent! %I = HX + HY + HX*HY - HXY; I = SXY - SX*SY; %fprintf( 'HX=%.4f HY=%.4f HX*HY=%.4f HXY=%.4f ->\t %.4f\n', bobX, bobY, ... % bobX*bobY, bobXY, I ); return; % % --------------------------------------- % function [v,Kp,Kpinv] = get_lin_ind_pts_K( pts, sigsq, ald_thresh, maxsize ) cnt = size(pts,2); v = randind( maxsize, cnt ); K = rbf4nn( pts(:,v), pts(:,v), sigsq ); dp = kdict_init( K, ald_thresh ); for i=2:maxsize dp = kdict( dp, i ); end; v = v( dp.Dict )'; Kp = dp.K; Kpinv = dp.Kinv; return; % % --------------------------------------- % function [v] = get_lin_ind_pts( pts, sigsq, ald_thresh, maxsize ) dp = dict_init( @rbf4nn, sigsq, ald_thresh, pts( :, 1) ); v(1) = 1; cnt = size( pts, 2 ); for i = 2:cnt dp = dict( dp, pts( :, i) ); if ( dp.addedFlag ) v( length(v)+1 ) = i; end; if ( size(dp.Dict,2) >= maxsize ) break; end; end; return;