How do I Estimate Joint Entropy Using a Histogram?

Question

I am trying to estimate the entropy for two time series, defined by random variables $X$ and $Y$, each distributed according to an unknown PDF which is to be estimated empirically (using a histogram approach, and binning data equally according to Sturge’s rule)

• In one dimension, I believe that we can use the estimated PDFs to then estimate the entropy of $X$, which is defined in this case for $n$ bins by: $H(X)= -\sum \limits_{i=1}^{n} P(x_i)\log P(x_i)$

  • Am I correct in my understanding that, considering the histogram approach, we can simply take the proportion of the total samples which lie within the bin as $P(x_i)$, multiplied by the log of the same, and perform a weighted sum across all bins using the midpoint value within each bin as the weight?

• However, to calculate the joint entropy between X and Y, we have multiple dimensions:$H(X,Y) = - \sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}p(x,y)\log p(x,y)$

  • I am not sure that performing the same procedure as above, only now in the $X$ and $Y$ direction, quite achieves this. Is the approach correct? Should we perhaps only consider the bins on the diagonal? (i.e. $i=j$)

• Finally, as I am performing an autoregressive analysis, I also have additional variables representing lagged time series (e.g. $X - \Delta$, $X - 2\Delta$ etc.). I therefore need to be able to perform the calculation in the N-dimensional case.

I am using Pandas data frames in Python to store the time series and perform the calculations, so answers with code samples are welcomed!

Answer 1

Joint entropy fundamentally relies on the joint distribution and there is no way to approximate it based on marginal distributions: histograms for each $X$ and $Y$. You need a joint histogram for $(X;Y)$ that is a 2D histogram with $n\times m$ bins.

For example, if $X$ and $Y$ are independent, $H(X,Y)=H(X)+H(Y)$ but whether they are independent or not is in the joint distribution and there is strictly no information about it in marginals.

If your dataset is large enough so that each $n\times m$ bin is far from empty, it may work. I don't think using the diagonal only is likely to produce anything useful.

Maybe you can read advanced methods for continuous entropy estimation that would work in dimension $N$. KDE is often used. I found this article.

Answer 2

Having spent more time working on this, it appears that the pdf can be estimated simply by considering the volume (in the 2d histogram case) of each bin as a proportion of the total volume. Since all bins in my implementation are the same size, I think I can use the height of each bar as a proportion of the total height (which equals number of samples N). Dealing with empty bins remains an outstanding question though, given the log term.

Using np.histogramdd() generalises to the many-dimensional case. I provide my current solution in the 2D case below, with sample time-series data. I still need to refactor the code to allow for a many-dimensional analysis, based on additional timelags, as the example doesn't yet consider lagged time series $Y - \Delta$, $X-\Delta$ etc.

import numpy as np
import pandas as pd
from itertools import tee, islice, chain
from numpy import ma
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


# Function to return midpoints of bins in any dimension
def get_midpoints(hist_bins):
    l_bound, u_bound = tee(hist_bins, 2)
    u_bound = chain(islice(u_bound, 1, None), [None])
    return [(i+j)/2 for i,j in zip(l_bound,u_bound) if j is not None]

# Dummy data
N = 1000
ix = np.arange(N)
sin = np.cos(2*np.pi*ix/float(N/2)).reshape(N,1)
cos = np.sin(2*np.pi*ix/float(N/2)).reshape(N,1)
columns = ['Y','X1']
df = pd.DataFrame(np.hstack((sin,cos)),columns=columns, index=ix)

# Split data into bins according to Sturges formula
y_bins = np.log2(len(df['Y']))+1
x_bins = np.log2(len(df['X1']))+1

# Create 2d histogram
Hist, [by, bx] = np.histogramdd((df['Y'],df['X1']),bins=(y_bins,x_bins), normed=False)

# Get midpoint of bins
midpoints_y = get_midpoints(by)
midpoints_x = get_midpoints(bx)

# Get width of bins
width_y = (np.max(midpoints_y) - np.min(midpoints_y))/len(midpoints_y)
width_x = (np.max(midpoints_x) - np.min(midpoints_y))/len(midpoints_x)

# PDF array across bins = area of histogram bin divided by total area (since bins equal size we can just use height)
p_xy = [[ Hist[y][x] for x in range(len(midpoints_x))] for y in range(len(midpoints_y)) ]
pdf = np.array(p_xy)/np.sum(p_xy) 


# Single entropy for each dimension H(X) = -sum(pdf(x) * log(pdf(x)))
H_Y = -np.sum( pdf.sum(axis=0) * ma.log2(pdf.sum(axis=0)).filled(0)) # Use masking to replace log(0) with 0
H_X1 = -np.sum( pdf.sum(axis=1) * ma.log2(pdf.sum(axis=1)).filled(0)) 

# Joint entropy H(X,Y) = -sum(pdf(x,y) * log(pdf(x,y))) 
H_XY =  -np.sum(pdf * ma.log2(pdf).filled(0)) # Use masking to replace log(0) with 0


# Output
print('Joint Entropy H(Y,X1): ' + str(H_XY))

plt.plot(df)
plt.xlabel('Time')
plt.legend(['Y','X1'])
plt.title('Two time series Y, X1')
plt.show()

# Define 3d axes
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

# Create an X-Y mesh of the 2D data
x_data, y_data = np.meshgrid(np.arange(Hist.shape[1]),np.arange(Hist.shape[0]))

# Flatten out the arrays to pass to bar3d
x_data = x_data.flatten()
y_data = y_data.flatten()
z_data = Hist.flatten()
ax.bar3d( x_data, y_data, np.zeros(len(z_data)), 1, 1, z_data)

# Finally, display the histogram
plt.show()