matlab实现传递熵计算（程序压缩优化）

之前写过一篇博客，参考了一个计算TE的matlab代码。见 matlab中实现1阶传递熵计算（利用最简单的概率密度估计方法）
最近有时间把代码简单优化了一下，减少了循环，提升计算效率。因为TE本身的计算复杂度就很高，因此，通过代码的优化来降低一些计算复杂度。同时，提升了代码的可读性，更容易理解。
下面是优化后的matlab程序：

%X，Y是时间序列，一维向量
%pieces是将数据离散化多少份，一般在数据量足够大的情况下分的越多，能够得到的信息越多
%j是用于计算的时间序列长度（大概选择2000，3000等即可）
function [txy, tyx]=transfer_entropy(X,Y,pieces,j)
%% parpare some data that use in reconstruction phase space
d_x=[];   select=1:j;
d_x(:,1) = X(select+1,1); d_x(:,2) = X(select,1); d_x(:,3)=Y(select+1,1); d_x(:,4) = Y(select,1);
X_max = max(X(:,1)); X_min = min(X(:,1)); Y_max = max(Y(:,1)); Y_min = min(Y(:,1)); 
delta1 = (X_max-X_min)/(2*pieces); delta2 = (Y_max-Y_min)/(2*pieces);
%% 计算概率 p(x(t)),p(x(t+1)),P(x(t+1),x(t)),p(x(t+1),x(t),y(t)),p(x(t),y(t)),p(y(t)),p(y(t+1)),p(y(t+1))
L1 = linspace(X_min+delta1,X_max-delta1,pieces);  L2 = linspace(Y_min+delta2,Y_max-delta2,pieces);

dist1 = zeros(pieces(1,1),2);
for q1 = 1:pieces
    k1 = L1(q1);k2=L2(q1);
    dist1(q1,1) = sum((d_x(:,2)>=(k1-delta1) )& (d_x(:,2)<=(k1+delta1)))
    dist1(q1,2) = sum((d_x(:,4)>=(k2-delta2) )& (d_x(:,4)<=(k2+delta2)))
end
dist1(:,1)=dist1(:,1)/sum(dist1(:,1)); dist1(:,2)=dist1(:,2)/sum(dist1(:,2));

dist2 = zeros(pieces(1,1),pieces(1,1),3);
for q1 = 1:pieces
    for q2 = 1:pieces
        k1 = L1(q1);k2=L1(q2);
        k3 = L2(q1);k4=L2(q2);
        dist2(q1,q2,1) = sum((d_x(:,1)>=(k1-delta1)) & (d_x(:,1)<=(k1+delta1)) & (d_x(:,2)>=(k2-delta1)) & (d_x(:,2)<=(k2+delta1)))
        dist2(q1,q2,2) = sum((d_x(:,3)>=(k3-delta2)) & (d_x(:,3)<=(k3+delta2)) & (d_x(:,4)>=(k4-delta2)) & (d_x(:,4)<=(k4+delta2)))
        dist2(q1,q2,3) = sum((d_x(:,2)>=(k1-delta1)) & (d_x(:,2)<=(k1+delta1)) & (d_x(:,4)>=(k4-delta2)) & (d_x(:,4)<=(k4+delta2)))
    end
end
dist2(:,:,1)=dist2(:,:,1)/sum(sum(dist2(:,:,1)));dist2(:,:,2)=dist2(:,:,2)/sum(sum(dist2(:,:,2)));dist2(:,:,3)=dist2(:,:,3)/sum(sum(dist2(:,:,3)));

dist3=zeros(pieces(1,1),pieces(1,1),pieces(1,1),2);
for q1=1:pieces
    for q2=1:pieces
        for q3=1:pieces
            k1=L1(q1);k2=L1(q2);k3=L1(q3);
            k4=L2(q1);k5=L2(q2);k6=L2(q3);
            dist3(q1,q2,q3,1) = sum((d_x(:,1)>=(k1-delta1)) & (d_x(:,1)<=(k1+delta1)) & (d_x(:,2)>=(k2-delta1)) & (d_x(:,2)<=(k2+delta1)) & (d_x(:,4)>=(k6-delta2)) & (d_x(:,4)<=(k6+delta2)))  
            dist3(q1,q2,q3,2) = sum((d_x(:,3)>=(k4-delta2)) & (d_x(:,3)<=(k4+delta2)) & (d_x(:,4)>=(k5-delta2)) & (d_x(:,4)<=(k5+delta2)) & (d_x(:,2)>=(k3-delta1)) & (d_x(:,2)<=(k3+delta1)))
        end
    end
end
dist3(:,:,:,1)=dist3(:,:,:,1)/sum(sum(sum(dist3(:,:,:,1)))); dist3(:,:,:,2)=dist3(:,:,:,2)/sum(sum(sum(dist3(:,:,:,2))));

%% 计算TE
sum_f_1=0; sum_f_2=0;
for k1 = 1:pieces(1,1)
    for k2 = 1:pieces(1,1)
        if dist2(k1,k2,2)~=0 && dist1(k2,2)~=0  
           sum_f_1=sum_f_1-dist2(k1,k2,2)*log2(dist2(k1,k2,2)/dist1(k2,2));
        end
        if dist2(k1,k2,1)~=0 && dist1(k2,1)~=0
           sum_f_2=sum_f_2-dist2(k1,k2,1)*log2(dist2(k1,k2,1)/dist1(k2,1));
        end
    end
end
disp(sum_f_1);disp(sum_f_2)
sum_s_1 = 0;sum_s_2 = 0;
for k1 = 1:pieces(1,1)
    for k2 = 1:pieces(1,1)
        for k3 = 1:pieces(1,1)
            if dist3(k1,k2,k3,2)~=0 && dist2(k3,k2,3)~=0
               sum_s_1=sum_s_1-dist3(k1,k2,k3,2)*log2(dist3(k1,k2,k3,2)/dist2(k3,k2,3));
            end
            if dist3(k1,k2,k3,1)~=0 && dist2(k2,k3,3)~=0
               sum_s_2=sum_s_2-dist3(k1,k2,k3,1)*log2(dist3(k1,k2,k3,1)/dist2(k2,k3,3));
            end
        end
    end
end
disp(sum_s_1); disp(sum_s_2)
txy = sum_f_1-sum_s_1;
tyx = sum_f_2-sum_s_2;
end

参考论文：
[1] Measuring information transfer.
[2] Finding the Direction of Disturbance Propagation in a Chemical Process Using Transfer Entropy.
[3] Detection of Cause-Effect Relations Based on Information Granulation and Transfer Entropy.