線形回帰を最小二乗法で解く - 機械学習に詳しくなりたいブログ

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

線形回帰とは何か

まず回帰分析というのは、目的変数 $y$ と説明変数 $x$ の間の関係を求めること。 $x$ と $y$ の間にどのような関係があるか、つまり $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{\Phi}$ は $\mathbf{x}$ を訓練データ数分、行で並べたものです。ここで $\mathbf{x}$ は、例えばM次多項式で近似する場合、0乗項も含めて要素数M+1の列ベクトルです。つまり $\mathbf{x} = (x^{0},x^{1},\ldots,x^{M-1})^{T}$ で、例えば $x=1$ なら $\mathbf{x} = (1,1,2,4,\ldots,1^{M-1})^{T}$ です。

入力 $x$ と観測値 $t$ の組み合わせは以下のグラフのようなものを用意しました。sin関数に適当なノイズを加えたものです。横軸が $x$ で縦軸が $t$ です。sin関数に適当なノイズを加えたものです。

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

M=3の結果（係数は、w0：1.8106e-01　w1 ：1.1390e+01　w2：-3.4349e+01　 w3：2.2887e+01でした。）

M=9の結果

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

Excelの近似式と比較

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

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

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

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の関数と捉えれば線形性が成り立っている