pythonでKL距離(KLダイバージェンス)

データを生成

N=10000個だけ正規分布、パレート分布（自由度10）、べき分布からサンプルを生成。

import matplotlib.pyplot as plt
import numpy as np

# サンプル数
N=10000

# 各分布からサンプルをN個生成
x = np.random.normal(size=N)
x2 = np.random.normal(size=N)

y = np.random.pareto(10, size=N)
y2 = np.random.pareto(10, size=N)

z = np.random.power(5, size=N)
z2 = np.random.power(5, size=N)

# 各分布から生成した点のヒストグラムをプロット
plt.figure(figsize=(10, 5))
hist, _ = np.histogram(x)
plt.plot(hist, label="normal")
hist, _ = np.histogram(x2)
plt.plot(hist, label="normal2")
hist, _ = np.histogram(y)
plt.plot(hist, label="pareto")
hist, _ = np.histogram(y2)
plt.plot(hist, label="pareto2")
hist, _ = np.histogram(z)
plt.plot(hist, label="powar")
hist, _ = np.histogram(z2)
plt.plot(hist, label="powar2")
plt.legend(title="distribution")
plt.grid()

KL-divergence

本来ならば aの分布（a_hist）に０が含まれると0にしないといけないが、 簡易的に分布全体に小さい値 epsilon=.00001 を足して計算できるようにする。

def KLD(a, b, bins=10, epsilon=.00001):
    # サンプルをヒストグラムに, 共に同じ数のビンで区切る
    a_hist, _ = np.histogram(a, bins=bins) 
    b_hist, _ = np.histogram(b, bins=bins)
    
    # 合計を1にするために全合計で割る
    a_hist = (a_hist+epsilon)/np.sum(a_hist)
    b_hist = (b_hist+epsilon)/np.sum(b_hist)
    
    # 本来なら a の分布に0が含まれているなら0, bの分布に0が含まれているなら inf にする
    return np.sum([ai * np.log(ai / bi) for ai, bi in zip(a_hist, b_hist)])

実行結果

x~x2, y~y2, z~z2間は同じ分布同士なので小さくなる。

print(KLD(x, y), KLD(y, z), KLD(z, x))
print(KLD(x, x2), KLD(y, y2), KLD(z, z2))

4.853129620306818 7.017965941238099 2.9077428243608163
0.003968976036639689 0.01426772591715473 0.012341556339116544

参考文献

Marsh, Charles. "Introduction to continuous entropy." Department of Computer Science, Princeton University (2013).

https://www.crmarsh.com/static/pdf/Charles_Marsh_Continuous_Entropy.pdf