MLaPP 演習問題解答(1) Section1 EX1.1-EX1.3

MLaPPとは「Machine Learning: A Probabilistic Perspective (Adaptive Computation and Machine Learning series) (English Edition) 1st Edition, Kevin P. Murphy (著)」という本のタイトルです。2000年以降機械学習の名著と言えばPRMLでしたが2010年代に刊行されたこの本はPRMLの内容を包含しつつ新しい内容に触れています。ベイズ統計の章の後に古典統計の章を設けてそれぞれの特性とp値がうまくworkしない事例を挙げているスタイルは独特です。筆者はこの本を1年半ほど読んでいますがやっと半分くらいです。その道中出来る限りで演習問題も解いてきたのでその解答を掲載して行く予定です。

この本の図や例題、演習問題のサンプルコードはGithubで参照できますがoctabeというmatlabのクローン言語で書かれています。昨今私も含めこの業界はPythonが主流なので筆者のコードを必要に応じてPythonに書き換えています。また車輪の再発明はとにかくしないという方針ですのでscikit-learn等の既存のライブラリが使える箇所は使用しています。

また上記著者Murphy氏のGithubhttps://github.com/probml/pmtk3を参考にしています。このリモートリポジトリ丸ごとローカルにダウンロードして必要なデータを適宜参照して下さい。

Contents

1 1章概要
2 演習問題

1章概要

章の内容を要約したものが以下です。詳細は本書を当たって下さい。

1章は諸々の基本事項の説明と応用事例に触れています。

演習問題

EX1.1 MNISTデータへK近傍法を適用する

import numpy as np
from sklearn.datasets import fetch_mldata
mnist=fetch_mldata('MNIST original')
#print mnist.data.shape
print type(mnist.data)
from sklearn.neighbors import KNeighborsClassifier
neigh = KNeighborsClassifier(n_neighbors=1)

# 学習用データN個，検証用データを残りの個数に設定
N = 60000
x_train, x_test = np.split(mnist.data,      [N])
y_train, y_test = np.split(mnist.target,    [N])
N_test = y_test.size

neigh.fit(x_train, y_train)

list=[]
#print(neigh.predict(x_test[100]))
for x,y in zip(x_test,y_test):
    a=neigh.predict(x)
    b=y
    if a[0]==b:
        print a[0],b
        x=1
    else:
        x=0
    list.append(x)
#print (y_test[100])
#list.append(if a==b)
print list
print 1-(float(list.count(1))/len(x_test))

import numpy as np

from sklearn.datasets import fetch_mldata

mnist=fetch_mldata('MNIST original')

#print mnist.data.shape

print type(mnist.data)

from sklearn.neighbors import KNeighborsClassifier

neigh = KNeighborsClassifier(n_neighbors=1)

# 学習用データN個，検証用データを残りの個数に設定

N = 60000

x_train, x_test = np.split(mnist.data, [N])

y_train, y_test = np.split(mnist.target, [N])

N_test = y_test.size

neigh.fit(x_train, y_train)

list=[]

#print(neigh.predict(x_test[100]))

for x,y in zip(x_test,y_test):

a=neigh.predict(x)

b=y

if a[0]==b:

print a[0],b

x=1

else:

x=0

list.append(x)

#print (y_test[100])

#list.append(if a==b)

print list

print 1-(float(list.count(1))/len(x_test))

scikit-learnの当該ライブラリを使います。何度も言いますが車輪の再発明はしません。

EX1.2 近似K近傍法

from pyflann import *
from numpy import *
from numpy.random import *
from sklearn.datasets import fetch_mldata
mnist=fetch_mldata('MNIST original')
# 学習用データN個，検証用データを残りの個数に設定
N = 60000
x_train, x_test = np.split(mnist.data, [N])
y_train, y_test = np.split(mnist.target, [N])
N_test = y_test.size
dataset = x_train
testset = x_test
flann = FLANN()
params = flann.build_index(dataset, target_precision=0.9, log_level = "info"); print params
result, dists = flann.nn_index(testset,1, checks=params["checks"]);
list=[]
for i,y in enumerate(y_test):
a=y_train[result[i]]
b=y
if a==b:
print a,b
x=1
else: 
x=0
list.append(x)
print 1-(float(list.count(1))/len(x_test))

from pyflann import *

from numpy import *

from numpy.random import *

from sklearn.datasets import fetch_mldata

mnist=fetch_mldata('MNIST original')

# 学習用データN個，検証用データを残りの個数に設定

N = 60000

x_train, x_test = np.split(mnist.data, [N])

y_train, y_test = np.split(mnist.target, [N])

N_test = y_test.size

dataset = x_train

testset = x_test

flann = FLANN()

params = flann.build_index(dataset, target_precision=0.9, log_level = "info"); print params

result, dists = flann.nn_index(testset,1, checks=params["checks"]);

list=[]

for i,y in enumerate(y_test):

a=y_train[result[i]]

b=y

if a==b:

print a,b

x=1

else:

x=0

list.append(x)

print 1-(float(list.count(1))/len(x_test))

pyflannというライブラリで近似K近傍法を実行できます。高速です。

EX1.3 k近傍法のクロスバリデーション

#!/usr/bin/env python

import os
import matplotlib.pyplot as pl
import numpy as np
from sklearn.neighbors import KNeighborsClassifier as KNN
from sklearn.cross_validation import cross_val_score
from utils.util import DATA_DIR


def load_data():
"""Since the knnClassify3c.mat is the matlab v7.3 or later file
we have to load data from txt"""
train_file = os.path.join(DATA_DIR, 'knnClassify3c', 'knnClassify3cTrain.txt')
test_file = os.path.join(DATA_DIR, 'knnClassify3c', 'knnClassify3cTest.txt')
train = np.loadtxt(train_file,
dtype=[('x_train', ('f8', 2)),
('y_train', ('f8', 1))])
test = np.loadtxt(test_file,
dtype=[('x_test', ('f8', 2)),
('y_test', ('f8', 1))])
return train['x_train'], train['y_train'], test['x_test'], test['y_test']

x_train, y_train, x_test, y_test = load_data()

#plot train fig
pl.figure()
y_unique = np.unique(y_train)
markers = '*x+'
colors = 'bgr'
for i in range(len(y_unique)):
pl.scatter(x_train[y_train == y_unique[i], 0],
x_train[y_train == y_unique[i], 1],
marker=markers[i],
c=colors[i])
pl.savefig('knnClassifyDemo_1.png')

#plot test fig
pl.figure()
for i in range(len(y_unique)):
pl.scatter(x_test[y_test == y_unique[i], 0],
x_test[y_test == y_unique[i], 1],
marker=markers[i],
c=colors[i])
pl.savefig('knnClassifyDemo_2.png')

x = np.linspace(np.min(x_test[:, 0]), np.max(x_test[:, 0]), 200)
y = np.linspace(np.min(x_test[:, 1]), np.max(x_test[:, 1]), 200)
xx, yy = np.meshgrid(x, y)
xy = np.c_[xx.ravel(), yy.ravel()]

# use the knn model to predict
for k in [1, 5, 10]:
knn = KNN(n_neighbors=k)
knn.fit(x_train, y_train)
pl.figure()
y_predicted = knn.predict(xy)
pl.pcolormesh(y_predicted.reshape(200, 200))
pl.title('k=%s' % (k))
pl.savefig('knnClassifyDemo_k%s.png' % (k))

#plot train err and test err with different k
ks = [1, 5, 10, 20, 50, 100, 120]
train_errs = []
test_errs = []
for k in ks:
knn = KNN(n_neighbors=k)
knn.fit(x_train, y_train)
train_errs.append(1 - knn.score(x_train, y_train))
test_errs.append(1 - knn.score(x_test, y_test))
pl.figure()
pl.plot(ks, train_errs, 'bs:', label='train')
pl.plot(ks, test_errs, 'rx-', label='test')
pl.legend()
pl.xlabel('k')
pl.ylabel('misclassification rate')
pl.savefig('knnClassifyDemo_4.png')

#cross_validate
scores = []
for k in ks:
knn = KNN(n_neighbors=k)
score = cross_val_score(knn, x_train, y_train, cv=5)
scores.append(1 - score.mean())
pl.figure()
pl.plot(ks, scores, 'ko-')
min_k = ks[np.argmin(scores)]
pl.plot([min_k, min_k], [0, 1.0], 'b-')
pl.xlabel('k')
pl.ylabel('misclassification rate')
pl.title('5-fold cross validation, n-train = 200')

#draw hot-map to show the probability of different class
knn = KNN(n_neighbors=10)
knn.fit(x_train, y_train)
xy_predic = knn.predict_proba(xy)
levels = np.arange(0, 1.01, 0.1)
for i in range(3):
pl.figure()
pl.contourf(xy_predic[:, i].ravel().reshape(200, 200), levels)
pl.colorbar()
pl.title('p(y=%s | data, k=10)' % (i))
pl.savefig('knnClassifyDemo_hotmap_%s.png' % (i))
pl.show()

100

101

102

103

104

#!/usr/bin/env python

import os

import matplotlib.pyplot as pl

import numpy as np

from sklearn.neighbors import KNeighborsClassifier as KNN

from sklearn.cross_validation import cross_val_score

from utils.util import DATA_DIR

def load_data():

"""Since the knnClassify3c.mat is the matlab v7.3 or later file

we have to load data from txt"""

train_file = os.path.join(DATA_DIR, 'knnClassify3c', 'knnClassify3cTrain.txt')

test_file = os.path.join(DATA_DIR, 'knnClassify3c', 'knnClassify3cTest.txt')

train = np.loadtxt(train_file,

dtype=[('x_train', ('f8', 2)),

('y_train', ('f8', 1))])

test = np.loadtxt(test_file,

dtype=[('x_test', ('f8', 2)),

('y_test', ('f8', 1))])

return train['x_train'], train['y_train'], test['x_test'], test['y_test']

x_train, y_train, x_test, y_test = load_data()

#plot train fig

pl.figure()

y_unique = np.unique(y_train)

markers = '*x+'

colors = 'bgr'

for i in range(len(y_unique)):

pl.scatter(x_train[y_train == y_unique[i], 0],

x_train[y_train == y_unique[i], 1],

marker=markers[i],

c=colors[i])

pl.savefig('knnClassifyDemo_1.png')

#plot test fig

pl.figure()

for i in range(len(y_unique)):

pl.scatter(x_test[y_test == y_unique[i], 0],

x_test[y_test == y_unique[i], 1],

marker=markers[i],

c=colors[i])

pl.savefig('knnClassifyDemo_2.png')

x = np.linspace(np.min(x_test[:, 0]), np.max(x_test[:, 0]), 200)

y = np.linspace(np.min(x_test[:, 1]), np.max(x_test[:, 1]), 200)

xx, yy = np.meshgrid(x, y)

xy = np.c_[xx.ravel(), yy.ravel()]

# use the knn model to predict

for k in [1, 5, 10]:

knn = KNN(n_neighbors=k)

knn.fit(x_train, y_train)

pl.figure()

y_predicted = knn.predict(xy)

pl.pcolormesh(y_predicted.reshape(200, 200))

pl.title('k=%s' % (k))

pl.savefig('knnClassifyDemo_k%s.png' % (k))

#plot train err and test err with different k

ks = [1, 5, 10, 20, 50, 100, 120]

train_errs = []

test_errs = []

for k in ks:

knn = KNN(n_neighbors=k)

knn.fit(x_train, y_train)

train_errs.append(1 - knn.score(x_train, y_train))

test_errs.append(1 - knn.score(x_test, y_test))

pl.figure()

pl.plot(ks, train_errs, 'bs:', label='train')

pl.plot(ks, test_errs, 'rx-', label='test')

pl.legend()

pl.xlabel('k')

pl.ylabel('misclassification rate')

pl.savefig('knnClassifyDemo_4.png')

#cross_validate

scores = []

for k in ks:

knn = KNN(n_neighbors=k)

score = cross_val_score(knn, x_train, y_train, cv=5)

scores.append(1 - score.mean())

pl.figure()

pl.plot(ks, scores, 'ko-')

min_k = ks[np.argmin(scores)]

pl.plot([min_k, min_k], [0, 1.0], 'b-')

pl.xlabel('k')

pl.ylabel('misclassification rate')

pl.title('5-fold cross validation, n-train = 200')

#draw hot-map to show the probability of different class

knn = KNN(n_neighbors=10)

knn.fit(x_train, y_train)

xy_predic = knn.predict_proba(xy)

levels = np.arange(0, 1.01, 0.1)

for i in range(3):

pl.figure()

pl.contourf(xy_predic[:, i].ravel().reshape(200, 200), levels)

pl.colorbar()

pl.title('p(y=%s | data, k=10)' % (i))

pl.savefig('knnClassifyDemo_hotmap_%s.png' % (i))

pl.show()

補足

ディレクトリの構成から上手くいかない場合は適宜変更してください。
必要なライブラリはインストールされているものとします。

参考文献

Kevin P. Murphy (2012) 『Machine Learning: A Probabilistic Perspective (Adaptive Computation and Machine Learning series) (English Edition) 1st Edition』, The MIT Press
https://github.com/probml/pmtk3

投稿者プロフィール

株式会社Crosstab 代表漆畑充

当ブログは【アナリティクスやデータの総合的な活用を通じて「未知の知」の獲得に貢献する】株式会社Crosstabが作成しています。

筆者について
2007年より金融機関向けデータ分析業務に従事。与信及びカードローンのマーケテイングに関する数理モデルを作成。その後大手ネット広告会社にてアドテクノロジーに関するデータ解析を行う。またクライアントに対してデータ分析支援及び提言/コンサルティング業務を行う。
統計モデルの作成及び特にビジネスアウトプットを重視した分析が得意領域である。
お問い合わせは株式会社Crosstabまでお願いいたします。