function [I] = topk(digit_x,digit_y,k);

[n,m] = size(digit_x);

I1 = find(digit_y==1);
I0 = find(digit_y==0);

Py1 = length(I1)/n;                   % P(y = 1) 
P1  = sum(digit_x(:,:))/n;            % P(x_i = 1)
% Should be normalized by # of y=1/0 examples, not by total length of the data
P11 = sum(digit_x(I1,:))/length(I1);  % P(x_i = 1|y = 1)
P10 = sum(digit_x(I0,:))/length(I0);  % P(x_i = 1|y = 0)

MI =    Py1*( P11.*log2(P11 ./P1)+(1-P11).*log2((1-P11)./(1-P1)) ) + ...
    (1-Py1)*( P10.*log2(P10 ./P1)+(1-P10).*log2((1-P10)./(1-P1)) );

% rank the components based on mutual information

[MItmp,I] = sort(-MI);

% select k most informative components

I = I(1:k);
