Rで機械学習する

2021-05-15 (土) 20:12:25 (134d) | Topic path: Top / 機械学習 / Rで機械学習する

はじめに

『Rによるバイオインフォマティクスデータ解析』を参考にして、機械学習を行います。

準備

Rのインストールについては,次のページを見てください.

データ

ここでは、標準で使用できるirisデータセットを使います。

data(iris)

このデータセットは、アヤメの種類(Species)をがく片の長さ(Sepal.Length)、幅(Lepal.Width)、花びらの長さ(Petal.Length)、幅(Petal.Width)によって分類する問題です。 長さと幅は連続値、種類はsetosa, versicolor, virginicaのいずれかをとる離散値です。

このデータセットには、setosa, versicolor, virginicaという3種類のアヤメについて、それぞれ50個ずつ、合計150個のデータが含まれています。 ランダムに10個のデータを選択して、見てみましょう。

iris[sort(sample(1:150,10)),]
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
4            4.6         3.1          1.5         0.2     setosa
22           5.1         3.7          1.5         0.4     setosa
65           5.6         2.9          3.6         1.3 versicolor
97           5.7         2.9          4.2         1.3 versicolor
100          5.7         2.8          4.1         1.3 versicolor
108          7.3         2.9          6.3         1.8  virginica
116          6.4         3.2          5.3         2.3  virginica
122          5.6         2.8          4.9         2.0  virginica
136          7.7         3.0          6.1         2.3  virginica
146          6.7         3.0          5.2         2.3  virginica

機械学習

機械学習 (machine learning) は、コンピューター(機械)に何かを学習させる技術の総称です。

ここでは、機械学習を次のように定義します。

入力ベクトル [math]X[/math] と出力値 [math]y[/math] の組 [math](X, y)[/math]の集合 [math]D[/math] が与えられたとき、すべての [math](X, y) \in D[/math] について次の式を満たす関数 [math]f[/math] を [math]D[/math] から求める。

\[f(X) = y\]

ここで、 [math]D[/math] の要素 [math](X, y)[/math] を事例 (example, instance) といいます。 つまり、[math]D[/math] は事例の集合です。 また、 [math]y[/math] は教師信号 (supervisory signal) と呼ばれます。

教師なし学習と教師あり学習

機械学習には教師なし学習 (unsupervised learning)、教師あり学習 (supervised learning)、強化学習 (reinforcement learning) があります。

教師信号が与えられない(わからない)機械学習を教師なし学習といい、教師信号が与えられる機械学習を教師あり学習といいます。 強化学習についてはここでは説明しません。

教師なし学習

クラスタリング

クラスタリングは、分類対象のデータ集合をいくつかのグループに分割するものです。 分割された部分データ集合をクラスターといいます。

クラスタリングの手法には、主に階層的アプローチ分割最適化アプローチがあります ここでは、後者の分割最適化アプローチの[math]k[/math]平均法だけを説明します。 前者の階層的アプローチについては、以下のページを参照してください。

[math]k[/math]平均法

[math]k[/math]平均法 ([math]k[/math]-means)は、データの集合を [math]k[/math] 個のクラスターに分割します。 クラスター数 [math]k[/math] は最初に決めておかなければなりません。 アルゴリズムの詳細は、次のページを参照してください。

まず、irisからSpeciesを除いた1列目から4列目だけを取り出します。

X <- iris[1:4]

4次元データですが、このうちのPetal.LengthとPetal.Widthだけを取り出して表示します。

x3 <- X$Petal.Length
x4 <- X$Petal.Width
plot(x3, x4)
iris.png

Rで[math]k[/math]平均法を用いるには、kmeans関数を用い、分割するデータとクラスター数 [math]k[/math] を指定します。 今回はカテゴリーが3つだとわかっているので、[math]k=3[/math] とします。

km <- kmeans(X, 3)

作成されたオブジェクトのcluster変数には割り当てられたクラスター番号が,centers変数には各クラスターの中心座標が格納されます.

km$cluster
  [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [39] 3 3 3 3 3 3 3 3 3 3 3 3 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [77] 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 2 2 1 2 2 2 2 2 2 1
[115] 1 2 2 2 2 1 2 1 2 1 2 2 1 1 2 2 2 2 2 1 2 2 2 2 1 2 2 2 1 2 2 2 1 2 2 1
km$centers
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1     5.901613    2.748387     4.393548    1.433871
2     6.850000    3.073684     5.742105    2.071053
3     5.006000    3.428000     1.462000    0.246000

クラスタリングの結果を表示します。 ここでは、色と形をクラスタリングによって得られたクラスター番号によって指定しています。

c <- km$cluster
plot(x3, x4, col=c, pch=c)
iris_kmeans.png

実際のSpeciesを色と形で表すと次のようになります。

y <- as.numeric(iris$Species)
plot(x3, x4, col=y, pch=y)
iris_species.png

[math]k[/math]平均法は、最初のクラスター番号がランダムに割り当てられるため、クラスター番号の順序には意味がありません。

教師あり学習

教師あり学習には、予測したいもの(教師信号 [math]y[/math])が数値である回帰 (regression) と、予測したいものがカテゴリーである分類 (classification) があります。

予測したいものを目的変数(統計の分野では従属変数)、予測に用いるものを説明変数(統計の分野では独立変数)といいます。

回帰分析

線形回帰

回帰分析の中で最もシンプルなものが線形回帰 (linear regression) です。 線形回帰では、次の式で予測モデルを表します。 \[ \hat{y} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p \] ここで、[math]\hat{y}[/math] は目的変数の予測値、[math]x_1, \dots, x_p[/math] は説明変数を表します。

説明変数が1つだけの線形回帰を単回帰、説明変数が2つ以上のものを重回帰といいます。

ここでは、花びらの長さ(Petal.Length)を説明変数、幅(Petal.Width)を目的変数にします。

x <- iris$Petal.Length
y <- iris$Petal.Width

これを散布図で見てみると、次のように分布しています。

plot(x, y)
iris_petal.png

回帰分析を行うには、lm関数モデル式を入力します。 モデル式は、記号 ~ を使って「目的変数 ~ 説明変数」というように表します。 説明変数が複数ある場合は、説明変数を + で繋げてしています。

r <- lm(y~x)

学習された結果を見てみると、次のようになります。

r
Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
    -0.3631       0.4158  

Coefficientsは係数を表し、Interceptは定数項の値(切片)を表します。 つまり、花びらの長さを [math]x[/math]、がくの長さを [math]y[/math] とすると、学習された式は次のように書けます。 \[ \hat{y} = -0.3631 + 0.4158 x \]

学習したモデル式を図示すると次のようになります。

plot(x, y)
abline(r)
iris_linear.png

実際には、複数の説明変数を用いた重回帰分析で、係数が0から離れた値にならないように正則化という方法を用います。 詳細は次のページを参照してください。

回帰分析には、分類分析用のSVMランダム・フォレストを回帰分析に用いる方法などもあります。

分類分析

ロジスティック回帰

ロジスティック回帰 (logistic regression) は、カテゴリーが2つだけのときに使える分類分析の方法です。

一方のカテゴリーに所属する確率 [math]p[/math] のロジット (logit) [math]\ln \left( \frac{p}{1-p} \right)[/math] を線形回帰によって予測します。 ロジットの予測式が \[\ln \left( \frac{\hat{p}}{1-\hat{p}} \right) = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \] のとき、このカテゴリーに所属する確率の予測式は次のように表されます。 \[ \hat{p} = \frac{1}{1+\mathrm{e}^{-(\beta_0 + \beta_1 x_1 + \dots + \beta_p x_p)}} \]

ここでは、irisの種類を目的変数、花びらの長さを説明変数とします。 irisには3種類のアヤメがありますが、ここではversicolorとvirginicaを使います。

x <- iris[51:150,]$Petal.Length
y <- as.numeric(iris[51:150,]$Species) - 2

as.numeric関数で数値に変換するとversicolorが2、virginicaが3になるので、それぞれ2を引いて0と1にしています。

これを図示すると次のようになります。

plot(x, y)
iris_vers_virg.png

ロジスティック回帰を使うには、glm関数を用い、family = binomialを指定します。

lr <- glm(y~x, family='binomial')

学習された結果を見てみると、次のようになります。

lr
Call:  glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
    -43.781        9.002  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    138.6 
Residual Deviance: 33.43 	AIC: 37.43

このモデルの [math]y = 1[/math] となる確率は次の式で表される。 \[ p = \frac{1}{1 + \mathrm{e}^{-(-43.781 + 9.002 x)}} \]

学習したモデルを図示するために、3から7まで0.1刻みのデータ列を用意し、このデータに対して学習したモデルに確率を予測させてグラフにします。

x_ <- seq(3, 7, 0.1)
y_ <- predict(lr, data.frame(x=x_), type='response')
plot(x, y)
lines(x_, y_)
iris_lr.png

このように、ロジスティック回帰のモデルはロジスティック関数(シグモイド関数)になります。 [math]x = 5.0[/math] のとき、この花がvirginicaである確率は [math]p = 0.774[/math] と予測されます。

決定木

決定木は、説明変数に対する条件を用いて事例集合を分割することによってラベル(カテゴリー)を予測する木を学習します。 実物を見た方が理解しやすいのでさっそく決定木を学習します。

決定木を学習するには、treeパッケージtree関数を用います。 まず、パッケージをインストールして、使えるようにします。 (インストールは1度だけで構いません。)

install.packages("tree")
library(tree)

使えるようになったら、決定木を学習します。

dt <- tree(Species~., iris)

ここでは、モデル式の指定方法がこれまでと違っています。 第2引数では、データを指定します。 第1引数では、目的変数と説明変数にそのデータの列を用いたモデル式を指定します。 説明変数が「.」(ドット)のときは、「目的変数以外の全ての列」という意味になります。

次に、学習された決定木を表示します。

dt
node), split, n, deviance, yval, (yprob)
      * denotes terminal node

 1) root 150 329.600 setosa ( 0.33333 0.33333 0.33333 )  
   2) Petal.Length < 2.45 50   0.000 setosa ( 1.00000 0.00000 0.00000 ) *
   3) Petal.Length > 2.45 100 138.600 versicolor ( 0.00000 0.50000 0.50000 )  
     6) Petal.Width < 1.75 54  33.320 versicolor ( 0.00000 0.90741 0.09259 )  
      12) Petal.Length < 4.95 48   9.721 versicolor ( 0.00000 0.97917 0.02083 )  
        24) Sepal.Length < 5.15 5   5.004 versicolor ( 0.00000 0.80000 0.20000 ) *
        25) Sepal.Length > 5.15 43   0.000 versicolor ( 0.00000 1.00000 0.00000 ) *
      13) Petal.Length > 4.95 6   7.638 virginica ( 0.00000 0.33333 0.66667 ) *
     7) Petal.Width > 1.75 46   9.635 virginica ( 0.00000 0.02174 0.97826 )  
      14) Petal.Length < 4.95 6   5.407 virginica ( 0.00000 0.16667 0.83333 ) *
      15) Petal.Length > 4.95 40   0.000 virginica ( 0.00000 0.00000 1.00000 ) *

テキストではわかりにくいので、図示します。

plot(dt)
text(dt)
iris_dt.png

最初に、一番上の条件「Petal.Lengthの値が2.45より小さいか」によって事例を分けます。 Petal.Lengthの値が2.45よりも小さい事例は左に、そうでない事例は右に振り分けられます。 左に振り分けたれた事例のカテゴリーは全てsetosaと予測されます。 右に振り分けられた事例は、2番目の条件「Petal.Widthの値が1.75より小さいか」によってさらに振り分けられます。

このように、決定木は学習されるルールがとてもわかりやすいという特徴を持っています。

最後に、table関数を用いて、予測と正解を表にします。

p <- predict(dt, iris, type="class")
table(p, iris$Species)
p            setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         47         1
  virginica       0          3        49

setosaは50個全て正しく分類できていますが、versicolorは50個中3個をvirginicaと、virginicaは50個中1個をversicolorと予測してしまっています。

決定木についての詳細は、次のページを参照してください。

SVM

SVM (Support Vector Machine) は、高次元空間において事例を分離する超平面(線形判別関数)を解析的に求めます。

SVMを用いるには、e1071パッケージのsvm関数を用います。 そこで、パッケージをインストールして、使えるようにします。

install.packages("e1071")
library(e1071)

使えるようになったら、SVMで予測モデルを学習します。

svm <- svm(Species~., iris)

学習された結果を見てみると、次のようになります。

Call:
svm(formula = Species ~ ., data = iris)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 

Number of Support Vectors:  51

SVMにはいくつかのタイプがあり、デフォルトではC-SVMという分類用のSVMになります。 costは、C-SVMのパラメーターです。 SVMは、カーネル (kernel) を使ってデータを高次元に写像して学習することができ、デフォルトではRBF (radial basis function) が使われます。 また、SVMはロジスティック回帰と同じカテゴリーが2つだけのときに使える手法ですが、機械学習用のライブラリーの多くはカテゴリーが3以上のときにはカテゴリーを2つずつ取り出してモデルを作成し、その結果をまとめることができます。

SVMが学習したモデルを図示するのは難しいです。 そこで、予測結果を図示します。

p = as.numeric(predict(svm,iris))
plot(x3, x4, col=p, pch=p)
iris_svm.png

最後に、table関数を用いて、予測と正解を表にします。

p <- predict(svm, iris)
table(p, iris$Species)
p            setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         48         2
  virginica       0          2        48

setosaは50個全て正しく分類できていますが、versicolorは50個中2個をvirginicaと、virginicaも50個中2個をversicolorと予測してしまっています。

SVMについての詳細は、次のページを参照してください。

ランダム・フォレスト

ランダム・フォレスト (random forest) は、データをランダムにサンプリングしてそこから決定木をランダムに作ることを繰り返し行い、たくさんの決定木の多数決によって予測を行います。 このような形で多数決によって予測する手法をアンサンブル学習 (ensemble learning) といいます。

ランダム・フォレストを用いるには、randomForestパッケージのrandomForest関数を用います。 そこで、パッケージをインストールして、使えるようにします。

install.packages("randomForest")
library(randomForest)
rf <- randomForest(Species~., iris)

ランダム・フォレストでは、説明変数の重要度を求めることができます。

varImpPlot(rf)
iris_rf_imp.png

これをみると、Petal.WidthとPetal.LengthはSepal.LengthとSepal.Widthに比べて重要であることがわかります。

ランダム・フォレストが学習したモデルを図示するのも難しいです。 そこで、予測結果を図示します。

p <- as.numeric(predict(rf, iris))
plot(x3, x4, col=p, pch=p)
iris_rf.png

最後に、table関数を用いて、予測と正解を表にします。

p <- predict(rf, iris)
table(p, iris$Species)
p            setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         50         0
  virginica       0          0        50

ランダム・フォレストは全ての事例を正しく分類できました。

過学習

ランダム・フォレストは全ての事例を正しく分類できましたが、一般的には、分類精度(汎化能力)が高いことはいいことであるとは限りません。

例えば、irisデータの23行目のSpeciesが間違ってversicolorになっていたとします。

iris[23,5] <- iris[51,5]

このデータを図示すると次のようになります。

y <- as.numeric(iris$Species)
plot(x3, x4, col=y, pch=y)
iris_overfit.png

明らかに、左下のデータのカテゴリーが間違っていますが、汎化能力が高いランダム・フォレストで学習すると、このデータもカテゴリーを間違ったまま学習してしまいます。

rf <- randomForest(Species~., iris)
p <- predict(rf, iris)
table(p, iris$Species)
p            setosa versicolor virginica
  setosa         49          0         0
  versicolor      0         51         0
  virginica       0          0        50

このように、汎化能力が高すぎるために、データ中のノイズまでも学習してしまい、学習された予測モデルが真のモデルからズレてしまうことを過学習 (overfit) といいます。

過学習を防ぐには、データを訓練用と検証用に分け、訓練用データから学習した予測モデルの予測精度を検証用データで評価して過学習しないよう機械学習の手法やそのパラメーターを選択します。

まとめ

機械学習は、与えられたデータから予測モデルを学習します。

機械学習には、教師信号が与えれらない教師なし学習と教師信号が与えられる教師あり学習(とここでは説明していない強化学習)があります。

教師あり学習の教師信号が連続値の場合を回帰といい、離散値(カテゴリー)の場合を分類といいます。

機械学習には複数の手法があり、目的やデータに応じて使い分けます。

機械学習の汎化性能が高すぎると、データに含まれるノイズも学習してしまう過学習という問題が生じることに気をつける必要があります。

演習

次のようにすると、irisデータを80%の訓練データ iris_train と残り20%の検証データ iris_valid に分けることができる。

data(iris)
idx <- sample(nrow(iris), as.integer(nrow(iris)*0.8))
iris_train <- iris[idx,]
iris_test <- iris[-idx,]

iris_trainから学習したモデルを用いてiris_validに対して予測を行い、予測精度を確認せよ。

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS