ふんわり R-tips

ぜんぜんわからない、俺たちは雰囲気でRをやっている

「最小二乗法による回帰分析」で理解する機械学習の初歩

「最小二乗法による回帰分析」で理解する機械学習の初歩

機械学習を使って予測

機械学習が何かと問われれば、様々な答えがあるかと思いますが、私は「未知のデータに対して予測を立てる」ための技術だと理解しています。過去のデータをもとにして、感や経験ではなく科学的なプロセスで予測・判断を行うことが機械学習のゴールです。

機械学習のアルゴリズムには様々なものがあり、これを使えば終わりという単純なものではありません。未知のデータにどれほど当てはまるのかを見つけ出すために、この分野では新たなアルゴリズムやアルゴリズムのパラメータを調整する方法がどんどん改善されています。

機械学習ののなかで最もシンプルで、基礎となるのが「最小二乗法による回帰分析」です。

最小二乗法による回帰分析

私はビールが大好きで、いつでもどこでもずっと飲んでいます。仕事終わりの冷たいビールはおいしいですよね。生産者からすると、作りすぎると余って困るし、足りないと私が激怒します。いつでも新鮮なビールを飲むためには、数か月先の需要を読みながら生産を続ける必要があります。

月別のビールの需要を予測したい

ビールの需要は毎月変化しますが、ある月に売れるビールの量は比較的安定しています。

http://www.brewers.or.jp/data/doko-list.html

ビール酒造組合が公開しているデータから2013年の月別出荷量は、

数量(KL)
1 129.152
2 167.694
3 202.965
4 228.356
5 224.145
6 263.399
7 302.885
8 262.091
9 198.010
10 225.693
11 222.533
12 316.087
x <- 1:12
y <- c(129.152, 167.694, 202.965, 228.356, 224.145, 263.399, 302.885, 262.091, 198.010, 225.693, 222.533, 316.087)
plot(x, y, xlim = c(0, 12), ylim = c(100, 350))

f:id:phmpk:20170412201817p:plain

夏場は多く、2月は日数が少ないので少なく、1月と12月は集計タイミングで凸凹になっているんでしょうか。それとも忘年会で大量消費されるんでしょうか。

図に示したデータはガタガタですが、その月々に固有のノイズが加わって、実際の出荷量が得られていると仮定します。この仮定が正しければ、ノイズを取り払ったなめらかな曲線を見出すことで、来年度の出荷量を予測することが可能になります。ノイズが加わるので、正確な予測にはなりませんが、確率的には下の図のような曲線状の値になる可能性が最も高くなるでしょう。この予測を使うことで、1年前の出荷量をそのまま使うよりも正確に需要を読んで生産することができます。

f:id:phmpk:20170412201854p:plain

fit <- lm(y ~ poly(x, 4))
xx <- seq(0,12,0.1)
lines(xx, predict(fit, data.frame(x=xx)), col = 2)

二乗誤差で曲線のあてはまらなさを評価

ここでさらに、上の図で示したような曲線が4次関数で表されるものと仮定します。


f(x) = w_0 + w_1 x + w_2 x^2 + w_3 x^3 + w_4 x^4 \tag{1} \label{1}


\eqref{1}x に1~12の値を代入すると、その月の需要予測が計算されます。w_0, \cdots, w_4 は未知のパラメータで、この値を変化させて、図のデータに最もフィットするものを導き出します。

最もフィットするものを探すための手がかりが、どれだけフィットしているのかを表す指標の一つである「二乗誤差」です。


E = \frac{1}{2} \sum_{n=1}^{12}{(\sum_{m=0}^{4} w_m n^m - d_n)^2} \tag{2} \label{2}


\sum_{n=1}^{12} は、1月~12月までの合計を表します。\sum_{m=0}^{4} w_{m} n^{m} は、4次関数 f(x)x=n を代入したものです。d_n は、n 月の出荷量を表します。

\eqref{2} で計算される値は、予測される値と実際に観測された値の差を二乗した値の合計です。二乗する理由は、正と負の値が混じっていても、あてはまらなさを数字の大きさで測れるようにするためです。\frac{1}{2} が出てきているのは、次の操作を簡単にするためのテクニック (微分すると係数が消えます) で、特別な意味があるわけではありません。

最小二乗法の計算過程

E を最も小さくする w_0, \cdots, w_4 を計算するために、次の連立方程式を解きます。


\frac{\partial E}{\partial w_m} = 0 \quad (m = 0, \cdots , 4) \tag{3} \label{3}


\eqref{3} は、w_0, \cdots, w_4 のそれぞれの変数で偏微分した E が0になることを表しています。多数の変数を持つ関数に対して、特定の変数で微分するのが偏微分です。ある関数が最大または最小値をとるときには、微分係数が0になります。

\eqref{2}E\eqref{3} に代入して、


\frac{\partial}{\partial w_m} \frac{1}{2} \sum_{n=1}^{12}{(\sum_{m=0}^{4} w_{m'} n^{m'} - d_n)^2} = 0 \tag{4} \label{4}


m' が出てきているのは、\eqref{4}m に0~4の値を代入して5回実行するときに、異なる意味の m が混ざるのを避けるためです。最初に定義した \eqref{2}mm' に置き換えました。

次に使うテクニックが「合成関数の微分」です。\eqref{4} で二乗しているカッコ内の式を以下の関数に置き換えます。


f_n(w_m) = \sum_{m'=0}^{4}{w_{m'} n^{m'} - d_n} \tag{5} \label{5}


m' の値を変えて足し合わせていくときに、w_m が出てくるので、式 \eqref{5}w_m の関数とみなしています。この関数を代入して、次式を得ます。


\frac{\partial}{\partial w_m} \frac{1}{2} \sum_{n=1}^{12}{\{ f_n(w_m) \}^2} = 0 \tag{6} \label{6}


微分と \sum の順番を入れ替えて、


\sum_{n=1}^{12}\frac{\partial}{\partial w_m} \frac{1}{2} {\{ f_n(w_m) \}^2} = 0 \tag{7} \label{7}


「合成関数の微分」の公式を適用します。「 w_m で微分していた式」を、「 f_n で微分して、 f_nw_m で微分したものを掛けた式」にします。


\frac{\partial}{\partial w_m} \frac{1}{2} {\{ f_n(w_m) \}^2} = \frac{d}{d f_n} (\frac{1}{2} {f_n}^2) \times \frac{d f(w_m)}{d w_m} \tag{8} \label{8}


\frac{d}{d f_n} (\frac{1}{2} {f_n}^2) の値は f_n となります。一方、 \frac{d f(w_m)}{d w_m} の値は、 m' が0~1の値をとるので、 w_{m'} のどれか1つが w_m に一致し、そのときの係数 n^{m} が1次関数の微分となります。


\sum_{n=1}^{12}{f_n \times n^m} = 0 \tag{9} \label{9}


f_n を代入して、次式を得ます。


\sum^{12}_{n=1}{\{ (\sum^{4}_{m'=0}{w_{m'} n^{m'} - t_n) n^m} \}} = 0 \tag{10} \label{10}


n^{m} でカッコを展開して、


\sum^{12}_{n=1}{\{ \sum^{4}_{m'=0}{w_{m'} n^{m'} n^m - d_n n^m} \}} = 0 \tag{11} \label{11}


\sum^{12}_{n=1} でカッコを展開して、


\sum^{4}_{m'=0} w_{m'} \sum^{12}_{n=1} n^{m'} n^m - \sum^{12}_{n=1} d_n n^m = 0 \tag{12} \label{12}


次の縦ベクトルと行列 \mathbf{w}, \mathbf{d}, \mathbf{\Theta} を使って、式を書き換えます。


\mathbf{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_4 \end{pmatrix}, \mathbf{d} = \begin{pmatrix} d_1 \\ d_2 \\ \vdots \\ d_{12} \end{pmatrix}



\mathbf{\Theta} = \begin{pmatrix} 
  1^0 & 1^1 & \cdots & 1^4 \\ 
  2^0 & 2^1 & \cdots & 2^4 \\
  \vdots & \vdots & \ddots & \vdots \\
  12^0 & 12^1 & \cdots & 12^4
\end{pmatrix}


\sum を使った和の計算は、次の行列で表された式となります。


\mathbf{w}^{\mathrm{T}} \mathbf{\Theta}^{\mathrm{T}} \mathbf{\Theta} - \mathbf{d}^{\mathrm{T}} \mathbf{\Theta} = 0 \tag{13} \label{13}


\mathrm{T} は行列の転置を表します。縦ベクトルが横ベクトルに、行列は左上から右下にかけて対角線で折り返したものになります。式を変形して、


\mathbf{w}^{\mathrm{T}} \mathbf{\Theta}^{\mathrm{T}} \mathbf{\Theta} = \mathbf{d}^{\mathrm{T}} \mathbf{\Theta} \tag{14} \label{14}


転置行列の公式 (\mathrm{AB})^{\mathrm{T}} = \mathrm{A}^{\mathrm{T}} \mathrm{B}^{\mathrm{T}}, (\mathrm{ABC})^{\mathrm{T}} = \mathrm{C}^{\mathrm{T}} \mathrm{B}^{\mathrm{T}} \mathrm{A}^{\mathrm{T}} を使い、両辺の転置行列を取って、


(\mathbf{\Theta}^{\mathrm{T}} \mathbf{\Theta}) \mathbf{w} = \mathbf{\Theta}^{\mathrm{T}} \mathbf{d} \tag{15} \label{15}


両辺に左から (\mathbf{\Theta}^{\mathrm{T}} \mathbf{\Theta})^{-1} をかけて、


\mathbf{w} = (\mathbf{\Theta}^{\mathrm{T}} \mathbf{\Theta})^{-1} \mathbf{\Theta}^{\mathrm{T}} \mathbf{d} \tag{16} \label{16}


これで最小二乗法による回帰分析の導出は完了です。\eqref{16} で近似ではない解析的な値が導出できたので、Rで値を計算してみます。solveは逆行列を計算するための関数で、行列の積は%*%演算子 (*は要素同士の積を計算するときに使います) です。

t <- matrix(c(1^0, 2^0, 3^0, 4^0, 5^0, 6^0, 7^0, 8^0, 9^0, 10^0, 11^0, 12^0,
              1^1, 2^1, 3^1, 4^1, 5^1, 6^1, 7^1, 8^1, 9^1, 10^1, 11^1, 12^1,
              1^2, 2^2, 3^2, 4^2, 5^2, 6^2, 7^2, 8^2, 9^2, 10^2, 11^2, 12^2,
              1^3, 2^3, 3^3, 4^3, 5^3, 6^3, 7^3, 8^3, 9^3, 10^3, 11^3, 12^3,
              1^4, 2^4, 3^4, 4^4, 5^4, 6^4, 7^4, 8^4, 9^4, 10^4, 11^4, 12^4), 
            12, 5)
w <- solve(t(t) %*% t) %*% t(t) %*% y
w
#             [,1]
# [1,] 164.2635707
# [2,] -56.3221728
# [3,]  35.1758205
# [4,]  -5.1476140
# [5,]   0.2243251

Rで回帰分析を行うlm関数の結果と比較してみましょう。lmの引数y ~ poly(x,4)で、目的変数 y を説明変数 x の4次多項式で回帰することを指定しています。計算結果はオブジェクトfitに格納されます。係数は、fit$coefficientsで参照できます。

fit <- lm(y ~ poly(x, 4))
fit$coefficients
# (Intercept) poly(x, 4)1 poly(x, 4)2 poly(x, 4)3 poly(x, 4)4 
#   228.58417   113.33279   -59.27529    73.70531    68.82609 

値が異なります。以下に示すように、関数の形はよく似ていて、ほぼ重なっています。係数のわずかな違いは、lmが実際には最小二乗法ではなく近似で解を計算しているために生じます。

fit2 <- function(x){return(w[1]+w[2]*x+w[3]*x^2+w[4]*x^3+w[5]*x^4)}
par(new=T)
plot(fit2, 0, 12, xlim = c(0, 12), ylim = c(100, 350), col = 3, ann = F)

f:id:phmpk:20170412201921p:plain

最小二乗法による回帰の次のステップ

こうして求まった曲線がうまくいかないと思ったときに、

  • 4次関数という仮定が誤っているのではないか? 5次関数のほうがいいかも?

  • 1年間のデータだけでなく、複数年のデータを使って回帰したほうがいいかも?

  • 最小二乗法以外の回帰アルゴリズムを使ったほうがいいかも?

といったことを考えていくことが、機械学習を使ったデータサイエンスになります。

まとめ

機械学習の初歩である最小二乗法を使った回帰の例を示しました。実際にここで書いたような数式の変換を手で動かすことはありません。Rでは、回帰を行うための関数が用意されています。最小二乗法のポイントは、重みを定めたときの「あてはまらなさ」を示す関数を用意して、その関数の値が最も小さくなるような解を求めていることです。計算のポイントは微分したときに値がゼロになる重みを、偏微分という方法を使うことです。こうすることで、他の変数の値を、とりあえず無視することができます。