『Python Machine Learning』Chapter.3をRでやってみた(前編)

@sfchaos氏がデータサイエンティスト養成読本機械学習入門編の振り返りと補足 - sfchaos blogで紹介している以下の本を写経している。

Python Machine Learning

単純パーセプトロンの解説・実装から始まり、scikit-learnによるクラス分類器の紹介、前処理、次元削減、モデル評価、アンサンブル学習と順を追って機械学習の手法が学べるので良書だと思う。Chapter.3のscikit-learnによるクラス分類器の紹介では以下のアルゴリズム毎にIrisのクラス分類の決定境界をmatplotlibで可視化しており直感的で分かりやすい。

• パーセプトロン
• ロジスティック回帰
• サポートベクターマシン
• 決定木
• ランダムフォレスト
• k近傍法

Rでもcaretとggplot2を使用して同じことができそうだったのでやってみた。

1. データの準備

Irisデータのロード

原著のPythonコードは以下の通り。Xはpetal lengthとpetal width、yは目的変数(species)

from sklearn import datasets
import numpy as np
iris = datasets.load_iris()
X = iris.data[:, [2, 3]]
y = iris.target

同じことをRだと以下のようになる

X <- iris[, c('Petal.Length','Petal.Width')]  # or iris[, 3:4](pythonの配列の添字は0始まりなので1ずつずれる)
y <- iris$Species  # or iris[, 5]

データの分割

データを学習用と検証用に分割してモデルの評価を行なう。原著のpythonコードは以下の通り。

from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)

Rの場合、caret::createDataPartitionで分割したレコードの行番号が取得できるので以下のように書くと同じ処理になる。（Pythonの方がシンプル…）

library(caret)
set.seed(71)
test <- createDataPartition(y, p=0.3, list=FALSE)

X_train <- X[-test, ]
y_train <- y[-test]
X_test <- X[test, ]
y_test <- y[test]

標準化

scikit-learn, caretともに学習用データセットの平均と標準偏差を用いて検証用データセットを標準化する関数が用意されている。 （手動でやると面倒なので地味に便利）

原著のコードは以下の通り。

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)

Rではcaret::preProcessで標準化できるが行名を削除しないといけないのが面倒。。。

preProc <- preProcess(X_train)
X_train_std <- predict(preProc, newdata=X_train)
X_test_std  <- predict(preProc, newdata=X_test)
rownames(X_train_std) <- NULL
rownames(X_test_std)  <- NULL

2. 決定境界プロット用の関数の作成

決定境界をプロットする関数を原著では以下のように定義している。

from matplotlib.colors import ListedColormap
import matplotlib.pyplot as plt

def plot_decision_regions(X, y, classifier,test_idx=None, resolution=0.02):
    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution), 
                           np.arange(x2_min, x2_max, resolution))

    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)

    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot all samples
    X_test, y_test = X[test_idx, :], y[test_idx]
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, 
                    c=cmap(idx), marker=markers[idx], label=cl)

    # highlight test samples
    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]
        plt.scatter(X_test[:, 0], X_test[:, 1], c='', alpha=1.0, 
                    linewidth=1, marker='o', s=55, label='test set')

meshgrid()やravel(), reshape(), contourf()については以下の書籍が参考になる。

データサイエンティスト養成読本 機械学習入門編 (Software Design plus)

若干見た目は変わるけどR版はggplot2で以下のようにした。テストデータセットは塗りつぶししないようにしている。

plot_decision_regions <- function(X, y, classifier,
                                  test_idx = NA, resolution = 0.02) {

  # setup marker generator and color map
  markers <- c(21, 22, 24, 23, 25) # ○, □, △, ◇, ▽
  ### 色はggplot2におまかせ
  
  # plot the decision surface
  x1_min <- min(X[, 1]) - 1
  x1_max <- max(X[, 1]) + 1
  x2_min <- min(X[, 2]) - 1
  x2_max <- max(X[, 2]) + 1
  
  xx <- expand.grid(seq(x1_min, x1_max, resolution),
                    seq(x2_min, x2_max, resolution))
  
  Z <- predict(classifier, newdata=xx)
  
  regions <- data.frame(x=xx[, 1], y=xx[, 2], z=Z)
  p <- ggplot(data=regions) + 
    geom_tile(aes(x, y, fill=z), alpha=0.4) +
    stat_contour(aes(x=x, y=y, z=as.integer(z)), size = 0.2, colour="grey50")
  
  # plot all samples
  points <- data.frame(x=X[, 1], y=X[, 2], z=y)
  if(is.na(test_idx)){
    p <- p +
      geom_point(data=points, aes(x=x, y=y, fill=z, shape=z), size=5) +
      scale_shape_manual(values = markers)
  } else {
    p <- p +
      geom_point(data=points[-test_idx, ], aes(x=x, y=y, fill=z, shape=z), size=5) +
      # highlight test samples
      geom_point(data=points[test_idx, ], aes(x=x, y=y, shape=z), size=5) +
      scale_shape_manual(values = markers)     
  }

  return(p)
}

長くなったので実際の決定境界のプロットは別記事にする。