function I=quad_renyi_fast(X,Y,xcov,ycov,NS) % function I=quad_renyi_fast(X,Y,xcov,ycov,NS) % % Quadratic Renyi mutual information between X and Y. % Defines mutual information as the sum of the individual % entropies, minus the joint entropy. % % 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 ); HX = norm_x * nsq * ( vX' * bXinv * vX ); % entropy of Y vY = sum( rbf4nn( sY, Y, ycov ), 2 ); HY = norm_y * nsq * ( vY' * bYinv * vY ); % joint entropy bXY = rbf4nn( sXY, sXY, 1 ); vXY = sum( rbf4nn( sXY, XY, 1 ), 2 ); HXY = (norm_x*norm_y) * nsq * (vXY' * pinv( bXY ) * vXY ); I = -log(HX) -log(HY) + log(HXY); 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;