機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

線形回帰を最小二乗法で解く

これも機械学習の1つなんですね。統計学のものだと思っていました。

線形回帰とは何か

まず回帰分析というのは、目的変数yと説明変数xの間の関係を求めること。yとxの間にどのような関係があるか、つまり y=f(x)で表されるとき、fはどう書けるか?を求めるもの。具体的な例で言えば、ある人の体重(x)が60kgだった場合、身長(y)は何cmか?みたいな感じですかね。もちろん何のヒントもなしじゃ予測しようがありませんから、複数のサンプルが与えられているわけですね。100人分の身長と体重のデータから関係式を求め、体重しかわからない人の身長を予測する、と。さて、回帰分析の中でも線形回帰とは、yとxの関係が例えば


y=w_{0}+w_{1}x+w_{2}x^{2}\tag{1}

のように表されるものを指すようです。で、係数wをどう設定すればyとxの関係が最も適切になるか、を求めることが線形回帰の目標です。これ、いきなり x^{2}とかになってて、全然線形じゃなくね、って思ったんですけど、線形っていうのは係数wに対してなんだそうです。*1 ややこしい! TensorFlow入門 線形回帰と非線形回帰の問題を解いてみたとか、TensorFlow入門~最小二乗法による非線形回帰は多分この勘違いをしていそうですよね。非線形回帰は、例えばy=w^{x}みたいなのだそうです。(参考

最小二乗法で解く

ある入力xに対する観測値tの組み合わせが複数与えられており、これを元に新しい入力xに対する予測値yの関係式を求めていきます。この値の組み合わせは訓練データとも呼びます。ここでyとxの関係は、

\begin{align}
y=w_{0}+w_{1}x+w_{2}x^{2}+\cdots+w_{M}x^{M}\tag{2}
\end{align}

でモデル化します。Mは解析者が任意に決定する値。求めるものは係数wです。この係数wは、最小二乗法で解析的に求めることが可能です。導出は、最小二乗法の解の導出に書きました。今回は結論だけ。係数wは以下の式で求められます。


\mathbf{w}=(\mathbf{\Phi}^{T}\mathbf{\Phi})^{-1}\mathbf{\Phi}^{T}\mathbf{t}\tag{3}

ここで\mathbf{w}は係数wの列ベクトル、\mathbf{t}は訓練データtの列ベクトルです。\mathbf{x}はM次多項式で近似する場合、0乗項も含めて要素数M+1の列ベクトルです。\mathbf{x}を訓練データ数分、行で並べたものが\mathbf{\Phi}です。

入力xと観測値tの組み合わせは以下のグラフのようなものを用意しました。sin関数に適当なノイズを加えたものです。 f:id:opabinia2:20180124214331p:plain

これをM=3、9でそれぞれ近似してみます。

M=3の結果(係数は、w0:1.8106e-01 w1 :1.1390e+01 w2:-3.4349e+01  w3:2.2887e+01でした。) f:id:opabinia2:20180124214328p:plain

M=9の結果 f:id:opabinia2:20180124214335p:plain

ちゃんと近似できていますね!M=9だと、訓練データ数に対して多項式の表現能力が高すぎて過学習を起こしていることも確認できました。この過学習という現象は、訓練データを丸暗記してしまったような状態で、学習によって賢くなったとは言えないんですね。もっと複雑なモデルで機械学習を行うと厄介な問題になるとのこと。なお過学習は正則化最小二乗法によって防ぐことができます。(参考:正則化最小二乗法

Excelの近似式と比較

ちなみに、最小二乗法による近似式の計算はExcelの機能にもあります。グラフのプロットを右クリックして「近似曲線の追加」を押せば、多項式近似のグラフと関係式が表示できます。結果はこうなりました。

f:id:opabinia2:20180124214338p:plain

Pythonで求めたM=3のときの係数と全く同じ値が表示されました。この近似式の表示は仕事でもけっこう使うことが多いんですけど、計算の中身を知ることができると面白いですね! 重回帰分析も最小二乗法で解いてみました→重回帰分析を最小二乗法で解く

今回使用したソースコード。

import numpy as np
import matplotlib.pyplot as plt


# y = w0*x^0+....wM*x^M を、引数xの配列数分求める
def y(w, x, M):
    X = np.empty((M + 1, x.size))
    for i in range(M + 1):
        X[i, :] = x ** i
    return np.dot(w.T, X)


# ランダムシードを固定
np.random.seed(0)
# 多項式の最大べき乗数(x^0+...+x^M)
M = 3
# 訓練データ数
N = 10

# 訓練データの列ベクトル
x = np.linspace(0, 1, N).reshape(N, 1)
# 訓練データtの列ベクトル
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)
t = t.reshape(N, 1)

# 行列Phiを作成
Phi = np.empty((N, M+1))
for i in range(M+1):
    Phi[:, i] = x.reshape(1, N) ** i

# 係数wの列ベクトルを解析的に求める
w = np.linalg.solve(np.dot(Phi.T, Phi), np.dot(Phi.T, t))

# 求めた係数wを元に、新たな入力x2に対する予測値y2を求める
x2 = np.linspace(0, 1, 100)
y2 = y(w, x2, M)

# 結果の表示
plt.xlim(0.0, 1.0)
plt.ylim(-1.5, 1.5)
plt.scatter(x, t)
plt.plot(x2, y2.T)
plt.show()

*1:参考:線形性  式(1)の多項式をwの関数と捉えれば線形性が成り立っている