2つのヒストグラムからKL情報量を計算する

2つの確率密度関数からKL情報量を計算できるなら，（経験的）ヒストグラム密度からでもKL情報量を計算できるじゃんというメモ．

離散分布のKL情報量

ヒストグラムを使うので離散でのKL情報量の計算になる．（カーネル密度推定すれば連続のようにも扱えるが，結局数学的に積分できない・無限の範囲を扱えないのでこっちになるだろう）

$$ D_{KL}(p||q) = \sum_{i} p(i) \log{(\frac{p(i)}{q(i)})} $$

$i$: ヒスグラムの各ビンのindex
$p, q$: 各ヒストグラムの密度

histogramの密度化

最初，足して1にすればいいと誤解していたので残しておく．実際は幅も考慮して，幅×カウントの総和が1になるようにすることが密度(density)化である．（ビン幅=1なら前者でも問題ない）

ビンの数を$n$として，ビン幅のベクトルを定義する．

$$ \vec{b} = ( bin_{1} - bin_{0}, \cdots, bin_{n} - bin_{ n-1 } ) = ( b_{1} , \cdots, b_{n}) $$

合計して1により以下の式を満たすことが密度化になる．

$$ \sum_{i=1}^{n} b_{i} p( i ) =1 $$

ビンのi番目の密度は，

${ p( b_{ i } ) = \frac{\text{counts}( b_{ i } ) }{ \sum{ \text{counts}( b_{ i } ) \cdot b_{ i } } } = \frac{ \text{counts}( b_{ i } ) }{ \sum{ \text{counts}( b_{ i } ) \cdot ( \text{bin}_{ i } - \text{bin}_{ i - 1 } ) } } }$

numpy.histogramの場合，density=Trueを引数に入れるだけで変換してくれる．

実装（コード）

gist: histogramからKL · GitHub

histogram作ってKL情報量を計算すればいいだけなので，簡単．問題は対象データがヒストグラムを作れるほど手に入るかどうか．

histogramの作成：今回はnumpy(scipy)のものを使用する．ビンを揃えるためにbinsの指定を必ず行う．デフォだと20分割．

また，0があると対数尤度が計算できなくなるので，小さい値を入れておく．この例では1e-9を密度に足してからさらに密度化する．

標準正規分布からhistogramを作ると以下になる．

import scipy as sp
import matplotlib.pyplot as plt
import scipy.stats as stats

sample1 = stats.norm(loc=0, scale=1).rvs(100)
sample2 = stats.norm(loc=0, scale=1).rvs(100)

def create_hists(data1, data2, bins, eps=1e-9):
    hist1, bins = sp.histogram(sample1, bins, density=True)
    hist2, bins = sp.histogram(sample2, bins, density=True)
    hist1 = (hist1 + eps)
    hist1 = hist1 / (sp.diff(bins) * hist1.sum())
    hist2 = (hist2 + eps)
    hist2 = hist2 / (sp.diff(bins) * hist2.sum())
    return hist1, hist2

hist1, hist2 = create_hists(sample1, sample2, bins)

KL情報の計算は，numpy(scipy)で書くと簡単，1行で書ける．密度$p$と対数尤度の差をベクトルに見立て，内積を取れば総和になる．

def D_KL(p, q):
    return p.dot(sp.log(p/q))

例

正規分布や二項分布を使って計算してみる．合わせてヒストグラムも描いてみる．

標準正規分布の値をヒストグラムにしてKLを計算

# sample 1000
import scipy as sp
import matplotlib.pyplot as plt
import scipy.stats as stats

N = 1000
sample1 = stats.norm(loc=0, scale=1).rvs(N)
sample2 = stats.norm(loc=0, scale=1).rvs(N)

hist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)

plt.bar(bins[:-1], hist1, label="sample 1")
plt.bar(bins[:-1], hist2, label="sample 2")
plt.legend()
plt.grid(True)
plt.show()

D_KL(hist1, hist2)

> 0.05993882587548316

標本の大きさ1000程度だと0.2-0.06ほどの値をとるのでこれをヒストグラム化すればKL情報量の分布も作れる．

$N=10,000$に標本を増やすと：

N = 10000
sample1 = stats.norm(loc=0, scale=1).rvs(N)
sample2 = stats.norm(loc=0, scale=1).rvs(N)

hist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)

plt.bar(bins[:-1], hist1, label="sample 1")
plt.bar(bins[:-1], hist2, label="sample 2")
plt.legend()
plt.grid(True)
plt.savefig("hist_10000.png")
plt.show()

D_KL(hist1, hist2)

回数を10倍増やすとKL情報量の精度も上がりおよそ1/10になる．

> 0.008646225250418648

二項分布と近似した正規分布

二項分布$Bin(n, p)$を正規分布に近似にすると $N(np, \sqrt{np(1-p)}^2)$．ただし試行回数パラメータnが十分大きいときに有効．

Bin(n=10, p=1/2)，N(np=5, np(1-p)=10/4)のKLを求める：

# sample 10000
import scipy as sp
import matplotlib.pyplot as plt
import scipy.stats as stats

N = 10000
n = 10
p = 1/2
sample1 = stats.binom(n, p).rvs(N)
sample2 = stats.norm(loc=n*p, scale=sp.sqrt(n*p*(1-p))).rvs(N)

bins = sp.linspace(5-8,5+8, 20+1)
hist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)

plt.bar(bins[:-1], hist1, label="bin(10, 1/2)")
plt.bar(bins[:-1], hist2, label="$N(np, np(1-p))$")
plt.legend()
plt.grid(True)
plt.savefig("Bin(10,0.5)_hist.png")
plt.show()

D_KL(hist1, hist2)

> 0.31268925739647135

正規分布への近似がよくないので0に近づかない．

二項分布Bin(n=10, p=1/2)と近似正規分布間のKL(N=10,000)

試行回数$n=10,000$回にすると

N = 10000
n = 10000
p = 1/2
sample1 = stats.binom(n, p).rvs(N)
sample2 = stats.norm(loc=n*p, scale=sp.sqrt(n*p*(1-p))).rvs(N)

bins = sp.linspace(n*p-200,n*p+200, 200+1)
hist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)

plt.bar(bins[:-1], hist1, label="$Bin(10000, 1/2)$")
plt.bar(bins[:-1], hist2, label="$N(np, np(1-p))$")
plt.legend()
plt.grid(True)
plt.savefig("Bin(10000,0.5)_hist.png")
plt.show()

D_KL(hist1, hist2)