matlab中实现1阶传递熵计算（利用最简单的概率密度估计方法）

关于时间序列的因果推理，基于传递熵方法。

传递熵（transfer entropy）是一种基于概率分布和香农熵，得出时间序列之间的因果关系的方法。这种方法首先提出是在PRL上由T.Schreiber 提出的，现在已经广泛应用于神经学、经济金融、工业等各领域中。

什么是transfer entropy ？它其实就是一个条件分布带来的探测到时间序列间的不对称性。就是一个变量 j 的历史是否对于另一个变量 i 的预测是否有影响，有影响即表示 j 到 i 之间存在信息传递，传递熵的作用就是利用信息传递来衡量变量之间的因果关系。下面看看它是怎么计算的：

通过这两部分（信息熵/香农熵）的相减得到传递熵。 这里k和l两个参数，是时间序列的阶数，也就是需要考虑多少的历史（阶数越大，考虑的历史越多）。大部分情况下都会认为是马尔科夫的，即只考虑1阶，所以 k=l=1，比较常见。但是对于复杂过程，需要考虑高阶的情况。

当k=l=1时，通过最简单的概率密度估计方法来估计概率密度函数（基于直方图/插值），进而计算传递熵。

%X是一个时间序列，是一个一维列向量%Y是一个时间序列，是一个一维列向量
%pieces是将数据离散化多少份，一般在数据量足够大的情况下分的越多，能够得到的信息越多
%j是用于计算的时间序列长度
function [en_1_2,en_2_1]=transfer_entropy(X,Y,pieces,j)
%% parpare some data that use in reconstruction phase space
d_x=[];    sit=size(X);
% templ=randperm((sit(1,1)-1));
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);
%% calculate the data disturbution 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);
count=0;
for q1=1:pieces
    k1=L1(q1);k2=L2(q1);
    count=count+1;
    count1=0;count2=0;
    for i=1:j
        if d_x(i,2)>=(k1-delta1) && d_x(i,2)<=(k1+delta1)
            count1=count1+1;
        end
        if d_x(i,4)>=(k2-delta2) && d_x(i,4)<=(k2+delta2)
            count2=count2+1;
        end
    end
    dist1(count,1)=count1;dist1(count,2)=count2;
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);
        count1=0;count2=0;count3=0;
        for i1=1:j
                if d_x(i1,1)>=(k1-delta1) && d_x(i1,1)<=(k1+delta1) && d_x(i1,2)>=(k2-delta1) && d_x(i1,2)<=(k2+delta1)
                    count1=count1+1;
                end
                if d_x(i1,3)>=(k3-delta2) && d_x(i1,3)<=(k3+delta2) && d_x(i1,4)>=(k4-delta2) && d_x(i1,4)<=(k4+delta2)
                    count2=count2+1;
                end
                if d_x(i1,2)>=(k1-delta1) && d_x(i1,2)<=(k1+delta1) && d_x(i1,4)>=(k4-delta2) && d_x(i1,4)<=(k4+delta2)
                    count3=count3+1;
                end
        end        
        dist2(q1,q2,1)=count1;dist2(q1,q2,2)=count2;dist2(q1,q2,3)=count3;        
    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);
            count1=0;count2=0;
            for i1=1:j
                if d_x(i1,1)>=(k1-delta1) && d_x(i1,1)<=(k1+delta1) && d_x(i1,2)>=(k2-delta1) && d_x(i1,2)<=(k2+delta1) && d_x(i1,4)>=(k6-delta2) && d_x(i1,4)<=(k6+delta2)  
                    count1=count1+1;
                end
                if d_x(i1,3)>=(k4-delta2) && d_x(i1,3)<=(k4+delta2) && d_x(i1,4)>=(k5-delta2) && d_x(i1,4)<=(k5+delta2) && d_x(i1,2)>=(k3-delta1) && d_x(i1,2)<=(k3+delta1)  
                    count2=count2+1;
                end
            end
            dist3(q1,q2,q3,1)=count1;dist3(q1,q2,q3,2)=count2;
        end
    end
end
dist3(:,:,:,1)=dist3(:,:,:,1)/sum(sum(sum(dist3(:,:,:,1))));dist3(:,:,:,2)=dist3(:,:,:,2)/sum(sum(sum(dist3(:,:,:,2))));
%% use the dosturbution calculate thansfer entropy
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)
en_1_2=sum_f_1-sum_s_1;
en_2_1=sum_f_2-sum_s_2;
end