Python中神经元转移熵计算结果异常的排查请求

问题描述 #

我需要计算两个神经元之间的转移熵，目的是量化：相较于仅知晓神经元2的活动，神经元1的活动能在多大程度上降低对神经元2后续活动的不确定性。但我编写的Python函数得到了反向结果——本该是神经元2预测神经元1，实际输出却是神经元1预测神经元2，请求排查函数中的计算错误。

实现代码 #

def transferEntropy_average2(iiv1,iiv2,L=100):
    """The interpretation is: How much does iiv2 predict iiv1"""
    # iiv1 and iiv2 are time vectors which denote neuronal activity as time series
    # The idea is to discretizise them, and use that to calculate probabilities
    # of different activity levels
    
    #I call neuron1 iiv1, and neuron2 iiv2. We're looking for how much neuron2 predics neuron1
    
    #Set some parameters
    N_quantiles = 100
    PY = np.zeros(N_quantiles) #Initialize PY (not used)
    quantiles = np.linspace(0, 100, N_quantiles + 1)[:-1] #Initialize stuff 
    q1 = np.percentile(iiv1, quantiles) # What are the percentiles for discretization 
    q2 = np.percentile(iiv2, quantiles) # And for the other 

    #Allocate memory
    #Joint probability of neuron 2 exhibiting one activity level, and exhibiting it
    # at t-1 before this activity level.
    PY_joint = np.zeros((N_quantiles,N_quantiles,L))
    #PXY is a 2D joint probability matrix just like above, but with the activity
    # of neuron 1 as well, so basically p(Y_t,Y_t-i,X_t-i) doing this for L i's
    PXY = np.zeros((N_quantiles,N_quantiles*N_quantiles,L))
    
    for i in range(L+1,len(iiv1)): # I do this for multiple time shifts (L)
        quantile1 = np.searchsorted(q1, iiv1[i], side='right') - 1 #How active was neuron1 (quantized)
        PY[quantile1] += 1
        for j in range(L):
            quantile2 = np.searchsorted(q1, iiv1[i-(j+1)], side='right') - 1 #How active was neuron1 at t-i
            quantile3 = np.searchsorted(q2, iiv2[i-(j+1)], side='right') - 1 #How active was neuron2 at t-i
            PY_joint[quantile1,quantile2,j] = PY_joint[quantile1,quantile2,j] + 1
            #This mess below is just because I wanted to avoid 3D matrices. It's like a "striped" 2D matrix. 
            #Every column (axis=1) is like "at t-1, neuron1 showed quantile2 amount of activity, while neuron2 showed quantile3 amount of activity"
            PXY[quantile1,quantile2+(quantile3*N_quantiles),j] = PXY[quantile1,quantile2+(quantile3*N_quantiles),j] + 1
    
    
    #then normalize
    PY = PY / np.sum(PY) #Again, not used, but for debugging purposes
    
    transferSum = 0 #Initialize the summation of transfer entropy
    #Transfer entropy for multiple time-shifts can be simplified as the mean of transfer entropy for all shifts
    
    transferTrace = np.zeros(L) # I use this if I want to check at which shift the transfer entropy is the highest
    
    for i in range(L): # For every shift
        PY_joint[:,:,i] = PY_joint[:,:,i] / np.sum(PY_joint[:,:,i]) # normalize joint probability matrix because sum of all p(Y,X) should add up to one
        PXY[:,:,i] = PXY[:,:,i] / np.sum(PXY[:,:,i]) # Same here, normalize the joint probability matrix
        
        #This is the "marginal X" but is actually Y_t-L,X_t-L, basically the probability of all neuronal activities in one past shift, irrespective of what neuron2 at t showed
        pX = np.sum(PXY[:,:,i],axis=0)
        
        #Calculate transfer entropy
        HYIYt = 0 #Conditional entropy of neuron2 H(activity(neuron2)|activity(neuron2 at t-i))
        HYIYt_Xt = 0
        # Marginal probability of shifted Y
        pYt = np.sum(PY_joint[:,:,i], axis=(0)) 
        
        #Now we're all set to calculate conditional entropies.
        #Here, I calculate the conditional entropy of neuron2 conditioned on the activity of neuron2 at every earlier time
        # and the conditional entropy of neuron2 conditioned on the activity of neuron2 and neuron1 at every earlier time
        for n1 in range(N_quantiles):
            for n2 in range(N_quantiles):
                
                if pYt[n2] == 0 or PY_joint[n1,n2,i] == 0: # If p=0, then H=0
                    result = 0
                else:
                    result = PY_joint[n1,n2,i] * np.log2(PY_joint[n1,n2,i] / pYt[n2])
                HYIYt = HYIYt - result
                for n3 in range(N_quantiles):
                    
                    if pX[n2+(n3*N_quantiles)] == 0 or PXY[n1,n2+(n3*N_quantiles),i] == 0:
                        result = 0
                    else:
                        result = (PXY[n1,n2+(n3*N_quantiles),i] * np.log2(PXY[n1,n2+(n3*N_quantiles),i] / pX[n2+(n3*N_quantiles)]))
                    HYIYt_Xt = HYIYt_Xt - result
        
        #Now we have calculated both conditional entropies, and the transfer entropy can be calculated
        # I add these up in "transfer sum" so that we can calculate the mean transfer entropy later
        transferSum = transferSum + (HYIYt - HYIYt_Xt)
        transferTrace[i] = (HYIYt - HYIYt_Xt) # Save the transfer entopy at this time shift
    
    transfer_entropy = transferSum / L
    return transfer_entropy

错误分析与修正建议 #

1. 目标与预测变量完全颠倒 #

函数注释明确要计算iiv2预测iiv1的程度，但代码全程将iiv1作为目标变量（未来时刻的神经元活动），iiv2作为辅助变量，完全搞反了因果方向：

• 循环中用iiv1[i]作为待预测的Y_t，实际应该用iiv2[i]作为目标（因为要量化的是神经元1对神经元2后续活动的预测能力）。
• PY_joint统计的是iiv1[t]与iiv1[t-j-1]的联合概率，正确的应该是iiv2[t]与iiv2[t-j-1]的联合概率（对应条件熵H(Y_t | Y_{t-j})）。

2. 条件熵物理含义混淆 #

代码注释称HYIYt是神经元2的条件熵，但实际计算的是神经元1的条件熵。需将所有目标变量替换为iiv2，过去的自变量替换为iiv2的历史，辅助变量替换为iiv1的历史：

• 正确的条件熵定义：
  • HYIYt = H(Y_t | Y_{t-j})，其中Y是iiv2
  • HYIYt_Xt = H(Y_t | Y_{t-j}, X_{t-j})，其中X是iiv1

3. 概率矩阵索引与维度错误 #

• PY_joint应存储(Y_t, Y_{t-j})的联合概率，即Y_t为iiv2当前值，Y_{t-j}为iiv2历史值。
• PXY的索引quantile2+(quantile3*N_quantiles)中，quantile2应为iiv2的历史值，quantile3为iiv1的历史值，对应联合状态(Y_{t-j}, X_{t-j})。

修正核心代码片段 #

# 修正循环中的变量赋值逻辑
for i in range(L+1, len(iiv2)):  # 以iiv2的长度为循环上限，目标是预测iiv2的后续活动
    # Y_t是iiv2的当前值（待预测的未来值）
    quantile1 = np.searchsorted(q2, iiv2[i], side='right') - 1  
    PY[quantile1] += 1
    for j in range(L):
        # Y_{t-j}是iiv2的历史值
        quantile2 = np.searchsorted(q2, iiv2[i-(j+1)], side='right') - 1  
        # X_{t-j}是iiv1的历史值
        quantile3 = np.searchsorted(q1, iiv1[i-(j+1)], side='right') - 1  
        # 统计(Y_t, Y_{t-j})的联合概率
        PY_joint[quantile1, quantile2, j] += 1  
        # 统计(Y_t, Y_{t-j}, X_{t-j})的联合概率
        PXY[quantile1, quantile2 + (quantile3 * N_quantiles), j] += 1

4. 归一化与条件熵计算一致性 #

修正目标变量后，需确保条件熵计算时的概率矩阵对应正确的联合分布：

• pYt应为PY_joint对Y_t维度求和得到的Y_{t-j}边缘概率
• pX应为PXY对Y_t维度求和得到的(Y_{t-j}, X_{t-j})联合边缘概率

内容的提问来源于stack exchange，提问作者FilterFeeder