motivation

Je veux calculer le montant d'informations mutuelles $ I (X; Y) $ des variables de probabilité continues $ X $ et $ Y $ en Python. $ 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):
    #Distribution de probabilité simultanée p(x,y)Calculs de
    p_xy, xedges, yedges = np.histogram2d(X, Y, bins=bins, density=True)
    
    # p(x)p(y)Calculs de
    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 et dy
    dx = xedges[1] - xedges[0]
    dy = yedges[1] - yedges[0]
    
    #Éléments d'intégration
    elem = p_xy * np.ma.log(p_xy / p_x_y)
    #Montant d'informations mutuelles et p(x, y), p(x)p(y)Production
    return np.sum(elem * dx * dy), p_xy, p_x_y

point

Si vous souhaitez calculer la quantité d'informations mutuelles pour le moment, vous pouvez utiliser la fonction ci-dessus. Soit dit en passant, je laisserai quelques points importants pour la mise en œuvre.

«Densité» de «np.histogram2d»

J'étais un peu impatient car np.sum (p_xy) ne devenait pas 1 quand je pensais vaguement que la probabilité serait renvoyée si densité = True était défini. Le point à noter est que «p_xy» est la ** densité de probabilité **, pas la probabilité.

Puisque $ X $ et $ Y $ sont des variables continues, l'approximation dans l'histogramme est la densité de probabilité. Par conséquent, si vous les additionnez en tenant compte de la largeur du bac, ce sera 1.

«np.histogram» et «np.histogram2d» renvoient la densité de probabilité et les intervalles (arêtes dans le code). Il est nécessaire de calculer «dx» et «dy» à partir de ce bac.

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)

#Si vous prenez la somme des densités de probabilité sans réfléchir, ce ne sera évidemment pas 1.
print(np.sum(p_x))  #Exemple de sortie: 1.580769264599771

#Si vous prenez la somme en tenant compte de la largeur du bac, elle devient 1.
dx = edges[1] - edges[0]
print(np.sum(p_x * dx))  #Exemple de sortie: 1.0000000000000002

Calcul de `p_x_y`

P_x_y dans le code essaie de calculer $ p (x) p (y) $. En fait, j'ai d'abord calculé avec le code suivant et cela n'a pas fonctionné.

p_x_y = p_x * p_y

Correctement

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

est. Dans le premier, p_x_y est le tableau principal, et dans le second, p_x_y est le tableau secondaire.

Exemple d'exécution

Exemple d'exécution 1 (deux ondes sin)

Comme ils ne sont pas indépendants, il y a une différence entre $ p (x, y) $ et $ p (x) p (y) $, et la quantité d'informations mutuelles augmente.

import matplotlib.pyplot as plt

#sin wave et cos wave
t = np.linspace(-5, 5, num=1000)
X = np.sin(2 * np.pi * t)
Y = np.cos(3 * np.pi * t)

#Calcul du montant d'informations mutuelles
mi, p_xy, p_x_y = mutual_information(X, Y, bins=30)

#Sortie de résultat
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()

Exemple d'exécution 2 (distribution normale indépendante)

Lorsque les deux variables sont indépendantes, $ p (x, y) $ et $ p (x) p (y) $ correspondent, et la quantité d'informations mutuelles devient petite.

import matplotlib.pyplot as plt
#Deux distributions normales indépendantes
N = 10000
X = np.random.normal(size=N)
Y = np.random.normal(size=N)

#Calcul du montant d'informations mutuelles
mi, p_xy, p_x_y = mutual_information(X, Y, bins=30)
Exemple d'exécution
#Sortie de résultat
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()

Calcul de la quantité d'informations mutuelles (valeur continue) avec numpy