python 如何计算转移熵？

这个工具箱里有：

https://github.com/nmtimme/Neuroscience-Information-Theory-Toolbox

上述的工具箱是matlab撰写的， 里边很多函数，没有整理；

下边的这个工具箱，相对完整，

Time Series Measures — PyInform 0.2.0 documentation

这个工具箱考虑了：

source  X， target Y

以及背景噪声：background，W

import numpy as np
from pyinform.transferentropy import *

k = 2
#生成两个随机序列
source = [0,1,1,1,1,0,0,0,0]
target = [0,0,1,1,1,1,0,0,0]
print("the tranfer entropy of X -> Y is {:.2f}".format(transfer_entropy(source, target, k=k)))

the tranfer entropy of X -> Y is 0.68

计算转移熵的时候，只能在0 1 序列上进行嘛？如果在普通连续信号（比如神经元的钙信号）上计算呢？

----------------------

还有一个博客，介绍了如何计算传递熵，写的还算详细：

机器学习——建立因果连系（传递熵）_Alphoseven的博客-程序员宅基地 - 程序员宅基地

最核心的一块是，计算联合条件概率和边缘概率分布。以及条件概率分布：

-----------------------

还有一个计算转移熵的仓库：

https://github.com/majianthu/pycopent/blob/master/copent/copent.py

#---
from scipy.special import digamma
from scipy.stats import rankdata as rank
from math import gamma, log, pi
from numpy import array, ndarray, abs, max, sum, sqrt, square, vstack, zeros
from numpy.random import normal as rnorm

##### calculating distance matrix
def dist(x, dtype = 2):
    (N,d) = x.shape
    xd = ndarray(shape = (N,N), dtype = float)

    for i in range(0,N):
        for j in range(i,N):
            if dtype == 1:
                xd[i,j] = sqrt(sum(square(x[i,:]-x[j,:])))
                xd[j,i] = xd[i,j]
            else: ## dtype = 2
                xd[i,j] = max(abs(x[i,:]-x[j,:]))
                xd[j,i] = xd[i,j]

    return xd
##### calculating the distance between two samples
def dist_ij(xi,xj, dtype = 2):
    if dtype == 1:
        return sqrt(sum(square(xi-xj)))
    else: ## dtype = 2
        return max(abs(xi-xj))
def construct_empirical_copula(x):
    (N,d) = x.shape
    xc = np.zeros([N,d])
    for i in range(0,d):
        xc[:,i] = rank(x[:,i]) / N

    return xc
##### constructing empirical copula density [1]
def construct_empirical_copula(x):
    (N,d) = x.shape
    xc = zeros([N,d])
    for i in range(0,d):
        xc[:,i] = rank(x[:,i]) / N

    return xc
##### Estimating entropy with kNN method [2]
def entknn(x, k = 3, dtype = 2, mode = 1):
    (N,d) = x.shape

    g1 = digamma(N) - digamma(k)

    if dtype == 1: # euciledean distance
        cd = pi**(d/2) / 2**d / gamma(1+d/2)
    else:  # maximum distance
        cd = 1;

    logd = 0
    if mode == 1: # for speed / small data
        distx = dist(x, dtype)
        for i in range(0,N):
            distx[i,:].sort()
            logd = logd + log( 2 * distx[i,k] ) * d / N
    else: # 2, for space / large data
        for i in range(0,N):
            disti = []
            for j in range(0,N):
                disti.append( dist_ij(x[i,:],x[j,:],dtype) )
            disti.sort()
            logd = logd + log( 2 * disti[k] ) * d / N

    return (g1 + log(cd) + logd)
##### 2-step Nonparametric estimation of copula entropy [1]
def copent(x, k = 3, dtype = 2, mode = 1, log0 = False):
    xarray = array(x)

    if log0:
        (N,d) = xarray.shape
        max1 = max(abs(xarray), axis = 0)
        for i in range(0,d):
            if max1[i] == 0:
                xarray[:,i] = rnorm(0,1,N)
            else:
                xarray[:,i] = xarray[:,i] + rnorm(0,1,N) * max1[i] * 0.000005

    xc = construct_empirical_copula(xarray)

    try:
        return -entknn(xc, k, dtype, mode)
    except ValueError: # log0 error
        return copent(x, k, dtype, mode, log0 = True)
##### conditional independence test [3]
##### to test independence of (x,y) conditioned on z
def ci(x, y, z, k = 3, dtype = 2, mode = 1):
    xyz = vstack((x,y,z)).T
    yz = vstack((y,z)).T
    xz = vstack((x,z)).T
    return copent(xyz,k,dtype,mode) - copent(yz,k,dtype,mode) - copent(xz,k,dtype,mode)
##### estimating transfer entropy from y to x with lag [3]
def transent(x, y, lag = 1, k = 3, dtype = 2, mode = 1):
    xlen = len(x)
    ylen = len(y)
    if (xlen > ylen):
        l = ylen
    else:
        l = xlen
    if (l < (lag + k + 1)):
        return 0
    x1 = x[0:(l-lag)]
    x2 = x[lag:l]
    y = y[0:(l-lag)]
    return ci(x2,y,x1,k,dtype,mode)

# source = np.array([[0,1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0,0]])
source = array([0,1,1,1,1,0,0,0,0])
target = array([0,0,1,1,1,1,0,0,0])

print("te is {:.2f}".format(transent(source, target)))

te is 1.09

----------------------

奇怪不同的方法， 算出来的结果差距这么大么？

-----------------

还有一个计算转移熵的代码：

https://github.com/nmtimme/Neuroscience-Information-Theory-Toolbox

% [te_result, ci_result, all_te] = ASDFTE(asdf, j_delay, i_order, j_order, windowsize)
% Parameters:
%   asdf        - Time series in Another Spike Data Format (ASDF)
%   j_delay     - Number of bins to lag sender (j series) or a vector [default 1]
%   i_order     - Order of receiver [default 1]
%   j_order     - Order of sender [default 1]
%   windowsize  - window size used for Coincidence Index calculation (odd number only)
%
% Returns:
%   te_result - (nNeu, nNeu) NxN matrix where N(i, j) is the transfer entropy from i->j
%               If multiple j_delay is given, this is a peak value of TE over delay.
%   ci_result - (nNeu, nNeu) NxN matrix where N(i, j) is the Coincidence Index from i->j
%               Multiple delays are necessary to calculate it.
%   all_te    - (nNeu, nNeu, delays) NxNxd matrix where N(i, j, k) is the transfer entropy
%               from i->j at delay of jdelay(k). For those who need all the delays.
%
% Examples:
% >> te = ASDFTE(asdf); % gives you delay 1 TE.
% >> [maxte, ci] = ASDFTE(asdf, 1:30); % gives you TEPk and TECI based on delay of 1 to 30.
% >> [maxte, ci] = ASDFTE(asdf, 1:30, 2, 3); % gives you HOTEPk and HOTECI (k=2, l=3) based on delay of 1 to 30.
% >> [maxte, ci, te] = ASDFTE(asdf, 1:30, 1, 1, 3); % gives you TEPk and TECI with CI window of 3, also gives you TE at all delays.

%==============================================================================
% Copyright (c) 2011, The Trustees of Indiana University
% All rights reserved.
%
% Authors: Michael Hansen (mihansen@indiana.edu), Shinya Ito (itos@indiana.edu)
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
%   1. Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%
%   2. Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%
%   3. Neither the name of Indiana University nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%==============================================================================

function [te_result, ci_result, all_te] = ASDFTE(asdf, j_delay, i_order, j_order, windowsize)
% Set defaults
if nargin < 2
    j_delay = 1;
end

if nargin < 3
    i_order = 1;
end

if nargin < 4
    j_order = 1;
end

if nargin < 5
    windowsize = 5;
end

if length(j_delay) == 1
    te_result = transent(asdf, j_delay, i_order, j_order); % Single delay
    return;
else
    % Multiple delays
    num_delays = length(j_delay);
    info = asdf{end};
    num_neurons = info(1);

    % Allocate space for all matrices
    all_te = zeros(num_neurons, num_neurons, num_delays);

    % Compute TE for delay times
    for d = 1:num_delays % Change this for to parfor for parallelization.
        all_te(:, :, d) = transent(asdf, j_delay(d), i_order, j_order);
    end

    % Reduce to final matrix
    te_result = max(all_te, [], 3); % reduction in 3rd dimension

    ci_result = zeros(num_neurons);
    if nargout > 1
        for i = 1:num_neurons
            for j = 1:num_neurons
                ci_result(i, j) = CIReduce(all_te(i, j, :), windowsize);
            end
        end
    end
end % if multiple delays
 

-----------------

此外，还有个工具箱，能够计算转移熵；

JIDT

 输入如下：

 

输出如下：

所以，从结果来看，还是推荐使用pyinfor，或者JIDT