Motivation

Ich möchte den gegenseitigen Informationsbetrag $ I (X; Y) $ der kontinuierlichen Wahrscheinlichkeitsvariablen $ X $ und $ Y $ in Python berechnen. $ I(X;Y) = \int_Y \int_X p(x, y) \log \frac{p(x,y)}{p(x)p(y)} dx dy $

Code

import numpy

def mutual_information(X, Y, bins=10):
    #Gleichzeitige Wahrscheinlichkeitsverteilung p(x,y)Berechnung von
    p_xy, xedges, yedges = np.histogram2d(X, Y, bins=bins, density=True)
    
    # p(x)p(y)Berechnung von
    p_x, _ = np.histogram(X, bins=xedges, density=True)
    p_y, _ = np.histogram(Y, bins=yedges, density=True)
    p_x_y = p_x[:, np.newaxis] * p_y
    
    #dx und dy
    dx = xedges[1] - xedges[0]
    dy = yedges[1] - yedges[0]
    
    #Elemente der Integration
    elem = p_xy * np.ma.log(p_xy / p_x_y)
    #Gegenseitige Informationsmenge und p(x, y), p(x)p(y)Ausgabe
    return np.sum(elem * dx * dy), p_xy, p_x_y

Punkt

Wenn Sie die Menge der gegenseitigen Informationen vorerst berechnen möchten, können Sie die obige Funktion verwenden. Im Übrigen werde ich einige wichtige Punkte für die Umsetzung belassen.

`Dichte` von `np.histogram2d`

Ich war etwas ungeduldig, weil "np.sum (p_xy)" nicht 1 wurde, als ich vage dachte, dass die Wahrscheinlichkeit zurückgegeben werden würde, wenn "Dichte = Wahr" gesetzt würde. Der zu beachtende Punkt ist, dass "p_xy" ** Wahrscheinlichkeitsdichte ** ist, nicht Wahrscheinlichkeit.

Da $ X $ und $ Y $ kontinuierliche Variablen sind, ist die Näherung im Histogramm die Wahrscheinlichkeitsdichte. Wenn Sie sie unter Berücksichtigung der Breite des Fachs addieren, ist dies 1.

np.histogram und np.histogram2d geben die Wahrscheinlichkeitsdichte und die Bins (Kanten im Code) zurück. Es ist notwendig, "dx" und "dy" aus diesem Bin zu berechnen.

import numpy as np

N = 1000
X = np.random.normal(loc=0, scale=1, size=N)

p_x, edges = np.histogram(X, bins=10, density=True)

#Wenn Sie die Summe der Wahrscheinlichkeitsdichten ohne nachzudenken nehmen, ist dies selbstverständlich nicht 1.
print(np.sum(p_x))  #Ausgabebeispiel: 1.580769264599771

#Wenn Sie die Summe unter Berücksichtigung der Behälterbreite nehmen, wird sie zu 1.
dx = edges[1] - edges[0]
print(np.sum(p_x * dx))  #Ausgabebeispiel: 1.0000000000000002

Berechnung von `p_x_y`

P_x_y im Code versucht $ p (x) p (y) $ zu berechnen. Eigentlich habe ich zuerst mit dem folgenden Code gerechnet und es hat nicht funktioniert.

p_x_y = p_x * p_y

Korrekt

p_x_y = p_x[:, np.newaxis] * p_y

ist. Im ersten Fall ist "p_x_y" das primäre Array und im zweiten Fall ist "p_x_y" das sekundäre Array.

Ausführungsbeispiel

Ausführungsbeispiel 1 (zwei Sinuswellen)

Da sie nicht unabhängig sind, gibt es einen Unterschied zwischen $ p (x, y) $ und $ p (x) p (y) $, und die Menge der gegenseitigen Informationen nimmt zu.

import matplotlib.pyplot as plt

#Sin Wave und Cos Wave
t = np.linspace(-5, 5, num=1000)
X = np.sin(2 * np.pi * t)
Y = np.cos(3 * np.pi * t)

#Berechnung der gegenseitigen Informationsmenge
mi, p_xy, p_x_y = mutual_information(X, Y, bins=30)

#Ergebnisausgabe
plt.figure(dpi=100)
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
ax1.set_title(r'$P_{XY}(x, y)$')
ax1.imshow(p_xy)
ax2.set_title(r'$P_{X}(x) P_{Y}(y)$')
ax2.imshow(p_x_y)
plt.suptitle('MI = {}'.format(mi))
plt.show()

Ausführungsbeispiel 2 (unabhängige Normalverteilung)

Wenn die beiden Variablen unabhängig sind, stimmen $ p (x, y) $ und $ p (x) p (y) $ überein, und die Menge der gegenseitigen Informationen wird gering.

import matplotlib.pyplot as plt
#Zwei unabhängige Normalverteilungen
N = 10000
X = np.random.normal(size=N)
Y = np.random.normal(size=N)

#Berechnung der gegenseitigen Informationsmenge
mi, p_xy, p_x_y = mutual_information(X, Y, bins=30)
Ausführungsbeispiel
#Ergebnisausgabe
plt.figure(dpi=100)
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
ax1.set_title(r'$P_{XY}(x, y)$')
ax1.imshow(p_xy)
ax2.set_title(r'$P_{X}(x) P_{Y}(y)$')
ax2.imshow(p_x_y)
plt.suptitle('MI = {}'.format(mi))
plt.show()

[PYTHON] Berechnung der gegenseitigen Informationsmenge (kontinuierlicher Wert) mit numpy