はしくれエンジニアもどきのメモ

情報・Web系技術・Englishの勉強メモ・備忘録です。

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

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

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

離散分布のKL情報量

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

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

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情報量の分布も作れる.

標準正規分布間のKL(N=1000)

$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

標準正規分布間のKL(N=10,000)

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

二項分布$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)

ビン幅や数でKL情報量は変わるが,試行回数を増やすと正規分布に近づいていることがわかる.

0.016068499410022707

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