Estimating mutual information of sine time series using 2D kernel density estimators

Question

This is an attempt to reproduce Moon et al. 1995, and the author's copy can be obtained through here.

As a benchmark, we estimate the time-lagged mutual information of a simple sine signal $\sin(0.02\pi t)$ using Gaussian kernels: $$K(y)=\dfrac{\exp\left(-\dfrac{(y-y_i)^TS^{-1}(y-y_i)}{2h^2}\right)}{(2\pi)^{d/2}h^d\det(S)^{1/2}},$$ where $h$ is chosen to be $$h=\left\{\dfrac{4}{(d+2)}\right\}^{1/d+4}n^{-1/(d+4)},$$ where $d$ is the dimension of the multivariate kernel, and $n$ is the number of data points.

I wrote a Matlab function (the code is at the end of the post).

Suppose $x_t$ is the original time series $\sin(0.02\pi t)$ with 400 data points $t=1,2,\cdots,400$, and $x_{t-\tau}$ is $\sin[0.02\pi (t-\tau)]$. The author plotted the estimated mutual information $I_{x_t,x_{t-\tau}}$ as a function of $\tau$:

However, what I got from my code is

While the qualitative features (like the diverging mutual information at lag 0, 50 and 100), but the magnitude is far off and the overall shape is not right.

My 1st suspicion is the fact that I used different window widths for joint pdf and marginal pdf, and they may not be consistent. So my question is: how to consistently construct joint and marginal pdf estimate in the context of mutual information estimation?

Second question is: is there anything wrong in the code besides the possibly inconsistent pdf KDEs?

Edit 1:

The reason why I choose KDE over other methods is that I think KDE is simple and intuitive in concept and possibly also simple in implementation. Besides MI, I also want to estimate transfer entropy and possibly other information-theoretic measures as well. I feel that estimating joint probability is general enough to apply on various information-theoretic measures. Please let me know if you have comment on my argument above.

I didn't use anything that can be found only in Matlab. Code in R, Python, etc. is welcome.

Edit 2: the Matlab code

The Matlab code:

function [MI] = get_MI(xt, xt_lag)
% xt is original time series, xt_lag is the lagged one, both are column
% vectors
x_pair = [xt' xt_lag']; n=length(xt);

d=1; h_d1=(4/(d+2))^(1/(d+4)) * n^(-1/(d+4));

d=2; h_d2=(4/(d+2))^(1/(d+4)) * n^(-1/(d+4));

MI = 0.;

Sxy_pair = cov(x_pair); invS_pair = Sxy_pair\eye(2); detS_pair = det(Sxy_pair);

Sxy_xt = cov(xt'); invS_xt = Sxy_xt\eye(1); detS_xt = det(Sxy_xt);

Sxy_lag = cov(xt_lag'); invS_lag = Sxy_lag\eye(1); detS_lag = det(Sxy_lag);

for ix=1:n
    Psq = p_mkde(x_pair(ix,:)', x_pair', h_d2, invS_pair, detS_pair);
    Ps = p_mkde(x_pair(ix,1), xt, h_d1, invS_xt, detS_xt);
    Pq = p_mkde(x_pair(ix,2), xt_lag, h_d1, invS_lag, detS_lag);

    MI = MI + (log2(Psq) - log2(Ps) - log2(Pq));
end

MI = MI/n;

end

function [pxy]=p_mkde(x,X,h,invS,detS)

s1=size(X);
d=s1(1);
N=s1(2);

pxy_sum=0;
for ix=1:N
    p2=(x-X(:,ix))'*invS*(x-X(:,ix));
    pxy_sum=pxy_sum+1/sqrt((2*pi)^d*detS)*exp(-p2/(2*h^2));
end
pxy=1/(N*h^d)*pxy_sum;
end

Answer 1

Let me start by restating that this is a terrible MI estimation procedure. Please don't use this in any serious business, unless you are absolutely sure you know what you are doing.

You are making two critical mistakes in your code:

1. Your scale doesn't match because you didn't do the integration properly. Mutual information is defined as in the paper is for probability, not probability density. But your kernel density estimator gives you back probability density which integrates to one, but does not sum to one on the points on a grid. If you want to use the sum formula as in the paper, you need to convert the density back to probability.

2. You are adding the probabilities only at the points where your data live, and counting multiple times on those regions. You should be integrating over the support of the entire distribution! An easy way to do this would be to make a uniform grid to integrate over. I think you should carefully review the definition of mutual information. (You need to add the probabilities, not the point wise MI estimates over the data.)

This is not a code-review site, but a couple of small things on your code.

• never use reserved matlab function names as variables: sum is one of them
• You don't need inv, so don't use it. Use \ and / instead
• try not to use i as a looping variable. i is for imaginary numbers in MATLAB (so is j)
• hint: MATLAB has a KDE for 1D which you can verify your code against: ksdensity

EDIT: I misinterpreted the code. I was wrong. OP is replacing expectation with sum over samples, so there's no problem using density nor having a grid.