[Python] Possible error with Transfer Entropy Calculation in Python

I want to calculate the transfer entropy between two neurons. I have two time-series vectors, representing the activity of the two neurons over time. I am using transfer entropy to answer the question "how much to does knowledge of neuron 1's activity reduce uncertainty of what neuron 2's activity is going to be at a later time point, compared to just knowing the activity of neuron 2".

I've written a function in python that seems to do the job, however I am getting results which suggest there may be some mistake which makes it seem neuron1 is predicting neuron2, where it should be the other way around.

I'd be super happy if someone could look at my function and see if they spot any errors in my calculation.

Here's a github link if you prefer reading it there: https://github.com/JonesBareBones/transferEntropy/tree/main

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, 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 collumn (axis=1) is like "at t-1, neuron1 showed quantile2 ammount of activity, while neuron2 showed quantile3 ammount 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 = (HYIYt - HYIYt_Xt) # Save the transfer entopy at this time shift

transfer_entropy = transferSum / L
return transfer_entropy

Continue reading...