mutual information estimation using KDE

カーネル密度推定(Kernel density estimation)はノンパラメトリックな確率密度関数の推定法である。 相互情報量(Mutual Information)の推定には、これを用いるのが最も基本的な方法である。 確率変数X, Y間の相互情報量は次式で定義される。

\displaystyle  \int_X \int_Y p(x, y) \log \frac{p(x, y)}{p(x)p(y)} dxdy \hspace{5em} (1)

ご覧のように、積分が出てきているので、まずはこれを潰さなければならない。 ここには前回、前々回で使ったサンプル平均の考え方を用いる。 さらに各確率密度関数を推定する必要があるが、今回は上述のようにカーネル密度推定により達成する。 最終的に、以下のように相互情報量を推定する。

\displaystyle  \hat{MI}(X, Y) := \frac{1}{n} \sum^n_{i=1} \log \frac{\hat{p}(x_i, y_i)}{\hat{p}(x_i)\hat{p}(y_i)}

実行結果を以下に示す。これは真の値は出せていないので、簡単な例で解釈する。 file://c:/Users/Hiroyuki/Dropbox/nkt1546789.github.io/python_ml/MI_1.png file://c:/Users/Hiroyuki/Dropbox/nkt1546789.github.io/python_ml/MI_3.png file://c:/Users/Hiroyuki/Dropbox/nkt1546789.github.io/python_ml/MI_2.png

このように、xとyになんらかの相関がある場合、相互情報量は高くなる。 逆になんの相関も無いならば、相互情報量は低くなる。 しかし、単に相関を測るなら相関係数で良いじゃないか、となる。 一般に普通の相関係数は線形な相関しか測れない。 非線形な相関を測りたいならば、相互情報量が有用だ、という話を聞いたことがあるのでその辺りを次回確認したい。 さらに、相関係数の非線形バージョンもあるみたいなので、その辺りも調べて行きたい。

実験に使ったコードを以下に示す。

import numpy as np
import pylab as pl
from scipy import stats
from sklearn.neighbors import KernelDensity
from sklearn.grid_search import GridSearchCV
from sklearn import metrics

n = 100
mean = [0, 0]
cov = [[1, 0.4], [0.4, 1]]
x = stats.multivariate_normal.rvs(mean=mean, cov=cov, size=n)

params = {"bandwidth": np.logspace(-1, 1, 100)}
kde_xy = GridSearchCV(KernelDensity(), params).fit(x).best_estimator_
kde_x = GridSearchCV(KernelDensity(), params).fit(np.c_[x[:,0]]).best_estimator_
kde_y = GridSearchCV(KernelDensity(), params).fit(np.c_[x[:,1]]).best_estimator_

MI = np.sum(kde_xy.score_samples(x) - kde_x.score_samples(np.c_[x[:,0]]) - kde_y.score_samples(np.c_[x[:,1]]))/n
pl.scatter(x[:,0], x[:,1])
pl.title("MI(X, Y) = {0}".format(MI))
pl.show()

Estimation of entropy of a random variable (2)

前回のポスト

で疑問にあげたシチュエーションを整理すると、

サンプル\{x_1, \ldots, x_n\}が与えられ、確率密度関数は未知、

となる。ここで確率密度関数を推定して、その確率密度関数から乱数を生成(これを手法1と呼ぶ)できるのか?というのが前回の疑問。 よく考えてみると、これは前回の問題とは別の問題だ。

今回のシチュエーションでは以下のような手法(手法2と呼ぶ)を適用できるだろう。 確率密度関数を推定して、\hat{p}(x)が得られたとする。 そして、すでに与えられているサンプルに対してこれを適用すれば、確率密度のサンプル\{\hat{p}(x_1), \ldots, \hat{p}(x_n)\}が得られる。 これによりエントロピーをサンプル近似することができる。

手法1が、任意の確率変数のエントロピーを求めること(目的1)を目的としているのに対して、 手法2は、i.i.d.標本の従う確率変数のエントロピーを求めること(目的2)を目的としている。 これは似ているようで全く違う。 目的1は目的2よりも圧倒的に難しい。 というわけで通常は目的2へアプローチすることの方が実用的である。

標本が正規分布から生成されたということが既知だとすると、そのパラメータを最尤推定することにより、確率密度関数を推定できる。 そして推定した確率密度関数とサンプルを使って、確率密度のサンプルを作る。 それでもってエントロピーをサンプル近似する。 そして前回見たように、これはサンプル数を増やせば、真の値に収束していく。

コードで書くと以下のようになる。

import numpy as np
import pylab as pl
from scipy import stats

# setting
n = 1000
mu = 0
sigma = 2
print np.log(sigma*np.sqrt(2*np.pi*np.e))

# given
x = stats.norm.rvs(mu, sigma, n)

# estimate
mu_hat = x.mean()
sigma_hat = x.std()
ds = stats.norm.pdf(x, mu_hat, sigma_hat)
print - np.sum(np.log(ds)) / n

estimation of entropy of Gaussian

連続型確率変数のエントロピーの計算法を調べてもなかなか良いのが出てこなかったので苦労した。 というわけで、それに関連するポスト。

確率変数Xのエントロピーは次式で定義される。

\displaystyle  h(X) = - \int^{\infty}_{-\infty} p(x)\log p(x) dx

今回の実験では、求めるエントロピーの真の値がわかっている必要があるので1変数正規分布を利用する。 1変数正規分布N(\mu, \sigma^2)に従う連続型確率変数のエントロピーは次式で求められる。

\displaystyle  ln\left(\sigma\sqrt{2\pi e}\right)

今回は標準正規分布N(0, 1)を用いる。 なおこの設定では当たり前だが、確率密度関数は既知だとする。

確率密度関数が既知なので、それに従う乱数をn個生成し、その確率密度を計算すると、n個の確率密度によるサンプル\{p(x_1), \ldots, p(x_n)\}を作ることができる。 つまりエントロピーを以下のようにサンプル近似する。

\displaystyle  \hat{h}(X) := - \frac{1}{n} \sum^n_{i=1} \log p(x)

実行結果は以下のようになった。

file://c:/Users/Hiroyuki/Dropbox/nkt1546789.github.io/python_ml/estimation_of_MI.png

ここから、サンプル数を増やしていけば、真の値に近づいていく、つまり、

\displaystyle  \hat{h}(x) \rightarrow h(x) \ (n \rightarrow \infty)

が言える。

結論として、エントロピーはサンプル近似することができるということがわかった。 ただし、ここで疑問がある。 通常は、サンプル\{x_1, \ldots, x_n\}のみが与えられていて、確率密度関数は未知である。 そこで、確率密度関数を推定して、\hat{p}(x)が得られたとする。 次は確率密度のサンプルを生成したいわけだが、これは可能なのだろうか。 すなわち任意の関数に従う乱数を生成することは可能なのだろうか。 これが疑問で仕方ないので、また調べる。

なお、使ったコードは以下である。

# coding: utf-8
import numpy as np
import pylab as pl
from scipy import stats

th = np.log(np.sqrt(2*np.pi*np.e))
print "theoritical:", th
ns = np.logspace(1, 6, 100)
hs = [np.sum(-np.log(stats.norm.pdf(stats.norm.rvs(0, 1, n))))/n for n in ns]
ths = np.repeat(th, 100)
pl.semilogx(ns, hs)
pl.semilogx(ns, ths, c="r")
pl.title("estimation of entropy of Gaussian N(0, 1)")
pl.xlabel("$n$")
pl.ylabel(r"$\hat{H}(X)$")
pl.show()

Ridge classification using python.sklearn

sklearn is a powerful module on python for machine learning. in this post, we introduce ridge regression. the solution of ridge regression is given by,

\displaystyle  \hat{\alpha} = (\Phi^T\Phi + \lambda I)^{-1}\Phi y .

the code is like this:

import numpy as np
import pylab as pl
from scipy import stats
from sklearn import linear_model

# make artificial data
n1 = 50
n2 = 50
x = np.c_[np.concatenate([stats.norm.rvs(-1, 1, n1), stats.norm.rvs(1, 1, n2)]),
          stats.norm.rvs(0, 1, 100)]
y = np.concatenate([np.repeat(-1., n1), np.repeat(1., n2)])
classifier = linear_model.RidgeCV().fit(x, y)

# for contour plot
offset = 0.1
xx, yy = np.meshgrid(np.linspace(np.min(x[:,0]-offset), np.max(x[:,0]+offset)),
                     np.linspace(np.min(x[:,1]-offset), np.max(x[:,1]+offset)))

Z = classifier.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
pl.scatter(x[:n1,0], x[:n1,1], c="b", marker="x")
pl.scatter(x[n1:,0], x[n1:,1], c="r", marker="o")
pl.contour(xx, yy, Z, levels=[0])
pl.show()
#import matplotlib.pyplot as plt
#plt.plot(x[:n1,0], x[:n1,1], "bx")
#plt.plot(x[n1:,0], x[n1:,1], "ro")
#plt.contour(xx, yy, Z, levels=[0])
#plt.show()

The results are shown below. linear classifier like this case is powerful, but it doesn’t work well if the data is not linearly separable (see graph2). The simple way to solve this problem is using kernel models. We’ll show the models and it work well even if the data is not linearly separable.

file://c:/Users/Hiroyuki/Dropbox/nkt1546789.github.io/python_ml/ridge_classification_1.png file://c:/Users/Hiroyuki/Dropbox/nkt1546789.github.io/python_ml/ridge_classification_2.png

統計的視点を捨てる

機械学習を勉強する上で統計的視点からみることは大事だと思う。

でもこれが非常に難しい。というのは統計的視点というのは解釈の連続だからだ。解釈の上にまた解釈がある。理解できない解釈があると、そこで止まってしまう。

そこで僕は統計的視点を捨てようと思う。完全に捨てるのではなく、できるだけ捨てる。

そんな僕にはイラストで学ぶ機械学習という本がベストマッチ。ほとんどの手法は最小二乗法の拡張であるという視点が僕を助けてくれる。