ビールと統計学

仕事でベイズ最適化を試す機会があったが、中身をちゃんと理解できていなかったので、有名な『ガウス過程と機械学習』を買ってみた。あくまで個人的な勉強ノート。

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

作者:持橋大地,大羽成征
講談社

Amazon

単回帰
重回帰
線形回帰

単回帰

$N$ 個の観測地のペア（ $D={(x_1, y_1), (x_2, y_2), \cdots, (x_N, y_N) }$ ）があったとき、 $x$ から $y$ を以下のような一次式で予測することを単回帰という

$\displaystyle \hat{y} = a + bx \qquad (a, b \in \mathbb{R})$

切片 $a$ と係数 $b$ の決定(最小二乗法)

単回帰の予測値 $\hat{y}_n$ と実測値 $y_n$ の誤差の二乗を最小化するように係数 $a, b$ を決める方法。誤差の二乗を $E$ とすると、

$\displaystyle \begin{align} E = \sum_{n=1}^{N} (y_n - \hat{y}_n)^2 &= \sum_{n=1}^{N} \bigl(y_n - (a + bx_n)\bigr)^2 \\ &= \sum_{n=1}^{N} (y_n^2+a^2+b^2x_n^2-2ay_n + 2abx_n - 2bx_ny_n) \end{align}$

と書ける。

$E$ が最小となるとき、 $a, b$ についての微分が0となるので、

$\begin{cases} \displaystyle \frac{\partial E}{\partial a} = 2 \sum_{n=1}^{N} (a-y_n+bx_n)=0\\ \displaystyle \frac{\partial E}{\partial b} = 2 \sum_{n=1}^{N} (bx_n^2+ax_n-x_ny_n) = 0 \end{cases}$

である。

これを整理して、 $\sum_{n} =\displaystyle \sum_{n}^{N}$ と略記すると以下が得られる。

$\begin{align} \begin{cases} aN &+ b\sum_{n} x_n &= \sum_{n} y_n\\ a\sum_{n} x_n &+ b\sum_{n}x_n^2 &= \sum_{n} x_ny_n \end{cases} \end{align}$

この式について、行列形式で書いた式を単回帰の正規方程式と呼ぶ。

$\begin{pmatrix} N & \sum_{n} x_n \\ \sum_{n} x_n & \sum_{n} x_n^2 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum_{n} y_n \\ \sum_{n} x_n y_n \end{pmatrix}$

Pythonでやってみる

非常に単純な単回帰用のデータセット $D$ として以下の場合を考える。

x	y
-2	1
1	5
4	6

グラフに出力すると下のようになる。（凡例の答え、 $y=3.17+0.83x$ はscikit-learnのLinearRegression()で計算した値）

単回帰の正規方程式にデータを代入すると、 $N=3$ であるから以下のようになる。

$\begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} 3 & 3 \\ 3 & 21 \end{pmatrix} ^{-1} \begin{pmatrix} 12 \\ 27 \end{pmatrix}$

逆行列をnumpyのlinalg.inv()で計算してやると以下の通り。上記の答えと同じ切片と係数が計算できることがわかる。

import numpy as np

A = np.array([[3, 3], [3, 21]])
A_inv = np.linalg.inv(A) #逆行列を計算
b = np.array([[12], [27]])
w = A_inv @ b

print(w)
# [[3.16666667]
#  [0.83333333]]

重回帰

入力が $D$ 次元のベクトル $(x_1, x_2, \cdots, x_D)^{T}$ のとき、計算の都合上、最初の次元が必ず1となるように新しく定義したベクトルを $\mathbf{x}^{T} = (1, x_1, x_2, \cdots, x_D)$ とする。ここで、切片を $w_0$ 、重み係数を $w_1, \cdots, w_D$ とおき、以下を求めることを重回帰という。

$\displaystyle \hat{y} = w_0 + w_1x_1 + \cdots +w_Dx_D = \mathbf{w}^T \mathbf{x} \qquad \text{where} \quad \mathbf{w}^T = (w_0, \cdots, w_D)$

ここで $N$ 個のデータについて並べると以下となる。ただし、 $\mathbf{w}^T\mathbf{x} = \mathbf{x}^T \mathbf{w}$ の関係を用いている。

$\begin{pmatrix} \hat{y}_1 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} \mathbf{w}^T\mathbf{x}_1 \\ \vdots \\ \mathbf{w}^T\mathbf{x}_N \end{pmatrix} = \begin{pmatrix} \mathbf{x}^T_1 \\ \vdots \\ \mathbf{x}^T_N \end{pmatrix}\mathbf{w} = \begin{pmatrix} 1 & x_{11} & \cdots & x_{1D} \\ \vdots & \vdots & & \vdots \\ 1 & x_{N1} & \cdots & x_{ND} \end{pmatrix}\mathbf{w}$

ここで、上式の一番右辺の左側の行列 $\mathbf{X}$ は計画行列と呼ばれ、あるn番目の水準における入力変数のベクトルを $\mathbf{x}_n=(x_{n1}, x_{n2}, \cdots, x_{nD})^{T}$ としたとき、 $\mathbf{x}_n^{T}$ を縦に並べたものである。

※要するに以下のようなデータセットが存在した場合に、入力の最初の次元(変数)として1を加えた入力変数のデータセットそのものである。

データセット $D$ ↓

$x_1$	$x_2$	$y$
1	2	4
-1	1	2
3	0	1

係数行列 $\mathbf{X}$ ↓

$\mathbf{X}= \begin{pmatrix} 1 & 1 & 2 \\ 1 & -1 & 1 \\ 1 & 3 & 0 \end{pmatrix}$

これを用いると重回帰は以下の式となる。

$\hat{\mathbf{y}} = \mathbf{X}\mathbf{w}$

重み係数ベクトル $\mathbf{w}$ の決定(最小二乗法)

ベクトル $\mathbf{x}$ の要素の二乗和は $\mathbf{x}^{T} \mathbf{x}$ と書けるので、誤差の二乗和は以下のようになる。 ※以下では、 $(\mathbf{AB})^{T} = \mathbf{B}^{T} \mathbf{A}^{T}$ の関係を用いて $\mathbf{y}$ について括りだしてまとめる変形を行っている

$\begin{align} \displaystyle E &= \sum_{n=1}^{N} (y_n - \hat{\mathbf{y}})^2 \\ & = (\mathbf{y} - \mathbf{Xw})^T (\mathbf{y} -\mathbf{Xw}) \\ & = \mathbf{y}^T\mathbf{y} - \mathbf{y}^T\mathbf{Xw} - (\mathbf{Xw})^T \mathbf{y} + (\mathbf{Xw})^T \mathbf{Xw} \\ & = \mathbf{y}^T\mathbf{y} - 2\mathbf{w}^T \mathbf{X}^T \mathbf{y} + \mathbf{w}^T \mathbf{X}^T \mathbf{Xw} \end{align}$

この $E$ が最小になるように、重みベクトル $\mathbf{w}$ で微分すると、以下のようになる(※)。

$\begin{align} \displaystyle \frac{\partial E}{\partial \mathbf{w}} &= -2\mathbf{X}^T\mathbf{y} + (\mathbf{X}^T\mathbf{X} + (\mathbf{X}^T\mathbf{X})^T)\mathbf{w} \\ &= -2\mathbf{X}^T\mathbf{y} + 2 \mathbf{X}^T \mathbf{Xw} \\ &= 0 \end{align}$

よって、重回帰モデルの正規方程式は

$\mathbf{X}^T\mathbf{Xw} = \mathbf{X}^T \mathbf{y}$

であるから、 $\mathbf{X}^T\mathbf{X}$ が逆行列を持つならば、重み係数 $\mathbf{w}$ は以下の式で求まる。

$\mathbf{w} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$

※ここではベクトルの一次式、二次形式の微分公式を使用している。

$\displaystyle \frac{\partial}{\partial \mathbf{w}} \mathbf{w}^T\mathbf{x} = (x_1, \cdots, x_D)^T = \mathbf{x} \\ \displaystyle \frac{\partial}{\partial \mathbf{w}}\mathbf{w}^T\mathbf{Aw} = (\mathbf{A} + \mathbf{A}^T)\mathbf{w}$

Pythonでやってみる

以下のようなデータセット $D$ を用いる

$x_1$	$x_2$	$y$
1	2	4
-1	1	2
3	0	1
-2	-2	-1

この重回帰をscikit-learnのLinearRegression()で計算すると、重み係数は以下のようになる。※ $w_0$ は切片として別に計算される。

import numpy as np
from sklearn.linear_model import LinearRegression

#データD
X = np.array([[1, 2],[-1, 1],[3,0],[-2,-2]])
y = np.array([4, 2, 1, -1])

#重回帰
model_lr = LinearRegression()
model_lr.fit(X, y)

w0 = model_lr.intercept_
W = model_lr.coef_

print(w0)
# 1.2018779342723005
print(W)
# [-0.01643192  1.20892019]

一方で、計画行列を用いた行列計算においても同様の結果を得ることができる。

import numpy as np

X = np.array([[1,1,2],[1,-1,1],[1,3,0],[1,-2,-2]])
print(X)
# [[ 1  1  2]
#  [ 1 -1  1]
#  [ 1  3  0]
#  [ 1 -2 -2]]
y = np.array([4, 2, 1, -1])

X_T = X.T #転置

#計算
X_T_X_inv = np.linalg.inv(X_T @ X)
w = X_T_X_inv @ X_T @ y

print(w)
# [ 1.20187793 -0.01643192  1.20892019]

線形回帰

線形回帰(1変数)

単回帰のように1つの入力変数 $x$ と出力変数 $y$ のより複雑な関係性を表現したい場合には、 $x$ の関数 $\phi_1(x), \phi_2(x), \cdots, \phi_{H}(x)$ を $H$ 個用意して、表現力を向上させた線形回帰モデルを用いる。

$\hat{y} = w_0 + w_1\phi_1(x) + \cdots + w_H\phi_H(x)$

※ここで、例えば $\phi_(x)$ に $x^2$ や $\sin(x)$ を用いることで、複雑な関係性を表現できる。

※「線形」回帰は入力 $\mathbf{x}$ と出力 $y$ の関係性が線形なのではなくて、重み $w_0, \cdots, w_H$ について線形和となっている、という意味である。

ここで $N$ 個のデータ（ $D={(x_1, y_1), (x_2, y_2), \cdots, (x_N, y_N) }$ ）について並べると、行列形式で以下のように書ける。

$\begin{pmatrix} \hat{y}_1 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} 1 & \phi_1(x_1) & \cdots & \phi_H(x_1) \\ \vdots & \vdots & & \vdots \\ 1 & \phi_1(x_N) & \cdots & \phi_H(x_N) \end{pmatrix} \begin{pmatrix} w_0 \\ \vdots \\ \vdots \\ w_H \end{pmatrix}$

となり、すなわち $\hat{\mathbf{y}} = \mathbf{\Phi w}$ となる。このとき $\mathbf{\Phi}$ が計画行列になる。

※ここでは $x_1, x_2, \cdots, x_N$ は入力変数 $x$ の具体的なデータ値であり、あくまで1変数 $x$ と $y$ の関係を表現しようとしていることに注意する。

線形回帰(多変数)

上記の線形回帰について、入力変数がベクトル $\mathbf{x}$ となった場合にも容易に拡張することができる。

$\phi_0(x) \equiv 1$ と定義すると、以下のようになる。

$\begin{pmatrix} \hat{y}_1 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} \phi_0(\mathbf{x}) & \phi_1(\mathbf{x}_1) & \cdots & \phi_H(\mathbf{x}_1) \\ \vdots & \vdots & & \vdots \\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_H(\mathbf{x}_N) \end{pmatrix} \begin{pmatrix} w_0 \\ \vdots \\ \vdots \\ w_H \end{pmatrix}$

※このとき $\mathbf{x}$ は入力変数の集合であることに注意する。

特徴ベクトルと基底関数

ここで、 $\mathbf{\phi}(\mathbf{x}) = (\phi_0(\mathbf{x}), \phi_1(\mathbf{x}), \cdots,\phi_H(\mathbf{x}))^{T}$ を特徴ベクトルと呼ぶ。線形回帰では、入力変数のベクトル $\mathbf{x}$ を $\mathbf{\phi}(\mathbf{x})$ によって線形変換して、以下の線形モデルを適用していることに相当する。

$\hat{y} = \mathbf{w}^{T} \mathbf{\phi}(\mathbf{x})$

※線形変換の結果、 $\mathbf{\phi}(\mathbf{x})$ が $\mathbf{x}$ よりも高次元になる場合もある。例えば、 $\mathbf{x} = x$ であり、 $\mathbf{\phi}(\mathbf{x}) = (1, x, x^2, \sin{x})^{T}$ のような場合である。

重回帰モデルは線形回帰モデルにおける $\mathbf{\phi}(\mathbf{x}) = \mathbf{x}$ となるような特別な場合であると捉えることができる。

一方で、見方を変えると、以下の線形モデル

$\hat{y} = \mathbf{w}^{T} \mathbf{\phi}(\mathbf{x}) = w_0\phi_0(\mathbf{x}) + w_1\phi_1(\mathbf{x}) + \cdots + w_H\phi_H(\mathbf{x})$

は、 $\mathbf{x}$ と $\hat{y}$ の関係を表す関数を、基底関数と呼ばれる、ベースとなる関数 $\hat{y} = \phi_0(\mathbf{x}), \quad \hat{y} = \phi_1(\mathbf{x}), \cdots$ の重み付き和として考えていることになる。

重み係数ベクトル $\mathbf{w}$ の決定(最小二乗法)

線形回帰モデルの解は、重回帰モデルと同じ形であらわされるため、以下となる。

$\mathbf{w} = (\mathbf{\Phi}^T\mathbf{\Phi})^{-1}\mathbf{\Phi}^T\mathbf{y}$

Pythonでやってみる

より実践的に、気象庁のHPで公開されている、世界の二酸化炭素(CO2)濃度の月別トレンドのグラフを線形回帰モデルで予測することを考える。

（出典：https://www.jma.go.jp/jma/kishou/books/hakusho/2020/csvindex.htm）

CO2濃度のトレンドには、周期的な（12か月単位の）季節変動がありつつも、全体としては上向きの曲線を示している。そこで、横軸である基準月からの月数を入力変数 $x$ とし、特徴ベクトルとして以下を考えた。全体のトレンドを2次関数で表現し、周期変動をcosで表現することを狙っている。

$\displaystyle \mathbf{\phi}(\mathbf{x}) = \Bigl(1, x, x^2, \cos \Bigl(\frac{\pi}{6} x + 2\Bigl) \Bigl)^T$

※12か月周期の変動を表すため、 $\pi x/6$ とし、周期変動のうちCO2濃度がピークとなる月とのずれを考慮して、 $\pi x/6+2$ と設定している。

import numpy as np
import pandas as pd

#計画行列を定義
def DesignMatrix(x):
    N = x.shape[0] 
    Phi_x = np.c_[np.ones(N).T, x, x**2, np.cos(np.pi/6*x+2)]
    return Phi_x

df = pd.read_csv('co2_trend.csv')
x = np.array(df['month']).reshape((-1,1))
y = np.array(df['co2']).reshape((-1, 1))

Px = DesignMatrix(x)
Px_T = Px.T #転置

Px_T_Px_inv = np.linalg.inv(Px_T @ Px)
w = Px_T_Px_inv @ Px_T @ y

print(w)
# [[ 3.44455477e+02]
#  [ 1.01239502e-01]
#  [ 1.22992739e-04]
#  [-2.26796433e+00]]

#予測値を計算
y_pred = Px @ w

結果として、線形回帰モデルは以下のようになった。

$\displaystyle \hat{y} = 344.5 + 0.101x + 0.000123 x^2 - 2.268 \cos \Bigl( \frac{\pi}{6} x +2 \Bigl)$

また、モデルと実際のデータを比較するとモデルが良好にトレンドを予測できていることがわかる。

ビールと統計学

勉強したことについてのらくがき帳

ガウス過程と機械学習ー1章、線形回帰モデル

単回帰

単回帰

切片 $a$ と係数 $b$ の決定(最小二乗法)

Pythonでやってみる

重回帰

重回帰

重み係数ベクトル $\mathbf{w}$ の決定(最小二乗法)

Pythonでやってみる

線形回帰

線形回帰(1変数)

線形回帰(多変数)

特徴ベクトルと基底関数

重み係数ベクトル $\mathbf{w}$ の決定(最小二乗法)

Pythonでやってみる

単回帰

単回帰

切片と係数の決定(最小二乗法)

Pythonでやってみる

重回帰

重回帰

重み係数ベクトルの決定(最小二乗法)

Pythonでやってみる

線形回帰

線形回帰(1変数)

線形回帰(多変数)

特徴ベクトルと基底関数

重み係数ベクトルの決定(最小二乗法)

Pythonでやってみる

切片 $a$ と係数 $b$ の決定(最小二乗法)

重み係数ベクトル $\mathbf{w}$ の決定(最小二乗法)

重み係数ベクトル $\mathbf{w}$ の決定(最小二乗法)