这是重现Moon等人的尝试。1995年，作者的副本可以通过这里获得。

作为基准，我们估计一个简单正弦信号的时滞互信息$\sin(0.02\pi t)$使用高斯核：

$$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}},$$ 其中选择为 其中是多元核的维数，是数据点的数量。$h$ $$h=\left\{\dfrac{4}{(d+2)}\right\}^{1/d+4}n^{-1/(d+4)},$$ $d$$n$

我写了一个 Matlab 函数（代码在文末）。

假设是原始时间序列有 400 个数据点，并且是。作者将估计的互信息 tau的函数：$x_t$$\sin(0.02\pi t)$$t=1,2,\cdots,400$$x_{t-\tau}$$\sin[0.02\pi (t-\tau)]$$I_{x_t,x_{t-\tau}}$$\tau$

但是，我从我的代码中得到的是

虽然定性特征（如滞后 0、50 和 100 处的发散互信息），但量级相差甚远，整体形状不正确。

我的第一个怀疑是我对联合 pdf 和边际 pdf 使用了不同的窗口宽度，它们可能不一致。所以我的问题是：如何在互信息估计的背景下一致地构建联合和边际 pdf 估计？

第二个问题是：除了可能不一致的 pdf KDE 之外，代码中还有什么问题吗？

编辑1：

我之所以选择 KDE 而不是其他方法，是因为我认为 KDE 在概念上简单直观，并且可能在实现上也很简单。除了 MI，我还想估计传输熵和可能的其他信息论度量。我觉得估计联合概率足以适用于各种信息论测量。如果您对我上面的论点有意见，请告诉我。

我没有使用任何只能在 Matlab 中找到的东西。欢迎使用 R、Python 等编写代码。

编辑 2：Matlab 代码

Matlab代码：

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