はじパタ4章 with Python

環境構築するだけしてnumpyすらまともに扱えないのは恥ずかしいので、はじパタのRによる実行例をPythonに移植しながらPythonの科学計算ライブラリに慣れつつ、機械学習にも慣れようという魂胆。
とりあえず第4章から。

irisのデータを用意

scikit-learnには様々な学習データが用意されていて、お決まりのようにその中にアヤメの学習データも含まれている。

from sklearn import datasets
iris = datasets.load_iris()
print iris.data

iris.dataが学習データの本体。
iris.feature_namesがそれぞれの列のラベル。

In [28]: iris.feature_names
Out[28]:
['sepal length (cm)',
 'sepal width (cm)',
 'petal length (cm)',
 'petal width (cm)']

ちなみに、各行の識別クラスはiris.targetで得られるし、iris.targetのクラス値と実際のラベルの対応関係はiris.target_namesで得られる。

In [33]: iris.target
Out[33]:
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

In [34]: iris.target_names
Out[34]:
array(['setosa', 'versicolor', 'virginica'],
      dtype='|S10')

この辺はTab補完で提供されている機能を確認できるIPythonが最高に便利。

データの標準化（P.39 実行例4.1）

ここでは花弁の長さと幅のデータだけを使っている。
花弁は英語でpetalなので、iris.dataの3,4行目のデータを使う。

In [43]: data = iris.data[:,2:4]

In [44]: data
Out[44]:
array([[ 1.4,  0.2],
       [ 1.4,  0.2],
       [ 1.3,  0.2],
       [ 1.5,  0.2],
       [ 1.4,  0.2],
       [ 1.7,  0.4],
       [ 1.4,  0.3],
       [ 1.5,  0.2],
       [ 1.4,  0.2],
       ...(略)

このように列を取り出すときはnumpy.ndarrayに対してコロンをうまく使う。

まずは平均ベクトルを出してみる。

In [61]: np.mean(data[:,0]), np.mean(data[:,1])
Out[61]: (3.7586666666666693, 1.1986666666666672)

Rだと行列のまま処理してくれるんだけどな…Pythonでもベクトルに切り出さずに処理できるうまい方法はないかな。
…と思ったらあった。

In [79]: np.mean(data, axis=0)
Out[79]: array([ 3.75866667,  1.19866667])

sum, var, sort, any, meanのように多くの関数がaxisという引数をとるらしい。
これは何かというと、

axis=0 : 縦方向
axis=1 : 横方向

を一つのデータとして、行列を扱うらしい。
今回の平均ベクトルを得る場合は、行列を縦に見るとpetal lengthとpetal widthで、それぞれの平均を得たいので、axis=0で縦方向と見なす。

次。共分散行列。

In [96]: np.cov(data, rowvar=0)
Out[96]:
array([[ 3.11317942,  1.29638747],
       [ 1.29638747,  0.58241432]])

rowvarオプションはデフォルトではnon-zeroになっていて、ようするにrowvar=trueということだから各行を変数として見ることになる。
だが、今回は例によって縦方向の各列を変数として見たいので、rowvar=0とする。
まあ行列を転置すれば同じことが実現できるが、不必要に転置するのはあまり好ましくないように思われるのでここではやらないことにする。
行列の転置はdata.Tである。

続いてメインのデータの標準化。
できるだけライブラリ関数でパッと出す方法がないかどうか調べていたところ、やっぱりあるらしい。
日本語では基準化変量とかいうのでstandardizedなんちゃらかと思っていたら、まさかのz-scoreだった。そりゃ見つかりませんわ。

In [110]: scipy.stats.zscore(data, axis=0)
Out[110]:
array([[ -1.34127240e+00,  -1.31297673e+00],
       [ -1.34127240e+00,  -1.31297673e+00],
       [ -1.39813811e+00,  -1.31297673e+00],
       [ -1.28440670e+00,  -1.31297673e+00],
       [ -1.34127240e+00,  -1.31297673e+00],
       [ -1.17067529e+00,  -1.05003079e+00],
       [ -1.34127240e+00,  -1.18150376e+00],
       [ -1.28440670e+00,  -1.31297673e+00],
       [ -1.34127240e+00,  -1.31297673e+00],
       [ -1.28440670e+00,  -1.44444970e+00],
       [ -1.28440670e+00,  -1.31297673e+00],
       [ -1.22754100e+00,  -1.31297673e+00],
       ...(略)

最後に標準化前と標準化後のデータをプロットしよう。
グラフのプロットはmatplotlibを使えばよく、散布図用のscatterという関数が用意されているのでそれを使う。
ここまでの流れを一旦スクリプトにまとめておく。

from sklearn import datasets
import scipy
import matplotlib.pyplot as pl

iris = datasets.load_iris()
data = iris.data[:, 2:4]
zdata = scipy.stats.zscore(data)

pl.scatter(zdata[:,0], zdata[:,1])
pl.show()

結果は以下の通り。

標準化前
f:id:levelfour:20140723115802p:plain

標準化後
f:id:levelfour:20140723121346p:plain

おお〜ちょっと感動。

無相関化（P.42 実行例4.2）

まずはさっき求めた共分散行列の固有値と固有ベクトルを求める。

In [13]: cov
Out[13]:
array([[ 3.11317942,  1.29638747],
       [ 1.29638747,  0.58241432]])

In [14]: v, S = np.linalg.eig(cov)

In [15]: v
Out[15]: array([ 3.65937449,  0.03621925])

In [16]: S
Out[16]:
array([[ 0.92154695, -0.38826694],
       [ 0.38826694,  0.92154695]])

covはさっき求めた共分散行列。numpy.linalg.eigは引数でもらった行列の固有値と固有ベクトルを計算して、tupleで返す関数。
なので、ここではvは固有値（のベクトル）、Sは固有ベクトル（を成分とする行列）になる。
なお、はじパタでは
$\Sigma=\begin{pmatrix}0.92&0.39\\0.39&-0.92\end{pmatrix}$
となっていてnumpyで得られたSの2列目と符号が違うけれど、固有ベクトルはそもそも一意に定まらないので問題はないだろう。

ここで得られたSを使って、dataを無相関化変換する。

In [47]: sdata = np.dot(S.T,data.T).T

In [48]: sdata
Out[48]:
array([[ 1.36781912, -0.35926433],
       [ 1.36781912, -0.35926433],
       [ 1.27566442, -0.32043764],
       [ 1.45997381, -0.39809103],
       [ 1.36781912, -0.35926433],
       [ 1.72193659, -0.29143502],
       [ 1.40664581, -0.26710964],
       [ 1.45997381, -0.39809103],
       [ 1.36781912, -0.35926433],
       [ 1.42114712, -0.49024572],
       [ 1.45997381, -0.39809103],
       ...(略)

まずは、numpyにおいて行列積はnumpy.dotを使うことに留意する。
単なる*(アスタリスク)は別の意味でオーバーロードされていて、単に行列の要素同士の積をそのままとる。
この場合はSとdataはそもそもshapeが違うのでエラーになる。
積の取り方がわかれば話は早く、
${\bf y}={\bf S}^T{\bf x}$
で線形変換すればよいのだから、上のようなコードになる。dataも転置をとっているのは、変換前の観測データは縦ベクトルにしないといけないからである。

あとは変換したsdataを先程と同じようにプロットするだけである。無相関化前のグラフは先程のグラフと同じだが、比較の意味を込めて再掲する。

無相関化前
f:id:levelfour:20140723115802p:plain

無相関化後
f:id:levelfour:20140723125748p:plain

コードはこのようになった。

from sklearn import datasets
import numpy as np
import scipy
import matplotlib.pyplot as pl

iris = datasets.load_iris()
data = iris.data[:, 2:4]
cov = np.cov(data, rowvar=0)
v, S = np.linalg.eig(cov)
sdata = np.dot(S.T, data.T).T

pl.scatter(sdata[:,0], sdata[:,1])
pl.ylim([-2.0, 3.0])
pl.show()

無相関化後した後の分布がx軸に平行になっていることがよくわかるように、y軸のrangeをある程度広めにとってみた。
デフォルトだとmatplotlibがうまく拡大してしまうので、ただの散布図になってよくわからない。
上のようにylimで調整してやればOK。

白色化（P.44 実行例4.3）

from sklearn import datasets
import numpy as np
import scipy
import matplotlib.pyplot as pl

iris = datasets.load_iris()
data = iris.data[:, 2:4]
mean = np.mean(data)			# 平均ベクトル
cov = np.cov(data, rowvar=0)	# 共分散行列
v, S = np.linalg.eig(cov)		# 固有値、固有ベクトル
# 対角化された行列（はじパタではΛ）
R = np.dot(np.dot(np.linalg.inv(S), cov), S)

# R ** (-1/2)
RSQ = np.linalg.inv(R ** .5)
# 白色化されたデータ
wdata = np.dot(np.dot(RSQ, S.T), (data - mean).T).T

# 白色化されたデータの共分散行列
print np.cov(wdata, rowvar=0)

pl.scatter(wdata[:,0], wdata[:,1])
pl.xlim([-3, 3])
pl.show()

詳しい証明ははじパタに譲るとして、やらなきゃいけないのは
${\bf u}={\bf \Lambda}^{-\frac{1}{2}}{\bf S}^T({\bf x}-{\bf \mu})$
という変換を行うこと。
コードでは若干わかりにくくなってしまっているが、上の数式をそのまま落とし込んでいるだけである。
実行例4.1、実行例4.2を実装・理解できたなら、これも問題ないはずである。