NumPyを使って線形モデルのパラメータを最小二乗法で推定する

線形モデルと最小二乗法について調べることがありましたので、まとめておきます。

今回は青本の2章を参考に、繰り返し測定、1元配置、多項式の3つのモデルで例示します。また、NumPyを使った行列計算でナイーブに最小二乗法を実装します。

自然科学の統計学 (基礎統計学)

作者: 東京大学教養学部統計学教室
出版社/メーカー: 東京大学出版会
発売日: 1992/08/01
メディア: 単行本
購入: 26人クリック: 308回
この商品を含むブログ (22件) を見る

線形モデル

ある観測値について未知母数の線形式＋誤差を仮定するモデルを、線形モデルと言います。誤差は正規分布 $N(0,\sigma^2)$ を仮定します。

繰り返し測定

飲料中のある成分の含有量を繰り返し測定した場合に以下の観測値が得られたとします。

No.	1	2	3	4	5	6	7	8	9
含有量(%)	37.55	37.74	37.69	37.71	37.73	37.81	37.53	37.58	37.74

こうした繰り返し実験では、以下のようにモデル化できます。

$$ y_i=\mu+\epsilon_i (i=1, \cdots , n) $$

1元配置

例えば、4種類の脂肪で揚げたドーナツの脂肪吸収量(g)を、各6回計測し、以下の値が得られたとします。

脂肪1	脂肪2	脂肪3	脂肪4
64	78	75	55
72	91	93	66
68	97	78	49
77	82	71	64
56	85	63	70
95	77	76	68

このようにa通りの条件で各r回の独立した実験を行った場合などは、1元配置のモデルとして以下のように表現できます。なお $1, \cdots , i , \cdots , a$ 、 $1, \cdots , j , \cdots , r$ です。

$$ y_{ij}=\mu+\alpha_i+\epsilon_{ij} $$

多項式

液体中のある成分の含有量によって曇点の変化を観測した実験を想定します。散布図にプロットしてみると、曲線を描いていることが確認できます。

含有量(%)	0	1	2	3	4	5	6	7	8	0	2	4	6	8	10	0	3	6	9
曇点(℃)	22.1	24.5	26.0	26.8	28.2	28.9	30.0	30.4	31.4	21.9	26.1	28.5	30.3	31.5	33.1	22.8	27.3	29.8	31.8

f:id:ohke:20180616095332p:plain

2次多項式を仮定すると、以下の式でモデル化できます。

$$ y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\epsilon_i $$

行列による一般化

観測値 $\vec{y}=(y_1, \cdots , y_n)^T$ 、未知母数 $\vec{\theta}=(\theta_1, \cdots , \theta_p)^T$ 、誤差 $\vec{\epsilon}=(\epsilon_1, \cdots , \epsilon_n)^T$ 、係数行列 ${\bf X}$ (n行p列のデザイン行列)を用いると、線形モデルは以下で一般化して表現できます。

$$ \vec{y} = {\bf X}\vec{\theta}+\vec{\epsilon} \Leftrightarrow \begin{pmatrix} y_{1} \\ \vdots \\ y_{n} \end{pmatrix}= \begin{pmatrix} x_{1,1} & \ldots & x_{1,p} \\ \vdots & \cdots & \vdots \\ x_{n,1} & \ldots & x_{n,p} \end{pmatrix} \begin{pmatrix} \theta_{1} \\ \vdots \\ \theta_{p} \end{pmatrix} + \begin{pmatrix} \epsilon_{1} \\ \vdots \\ \epsilon_{n} \end{pmatrix} $$

最小二乗法

観測値 $\vec{y}$ とその期待値の偏差二乗和(以下式)を最小にする解 $\hat{\vec{\theta}}$ を求めます。

$$ S(\vec{\theta})=||\vec{y} - {\bf X} \vec{\theta}||^2=(\vec{y} - {\bf X} \vec{\theta})^{T}(\vec{y} - {\bf X} \vec{\theta}) $$

繰り返し測定のモデルでは、 $S(\vec{\theta})=\sum(y_i-\mu)^2$
1元配置のモデルでは、 $S(\vec{\theta})=\sum\sum(y_{i,j}-\mu-\alpha_i)^2$
多項式のモデルでは、 $S(\vec{\theta})=\sum(y_i-\beta_0-\beta_1x_i-\beta_2x_i^2)^2$

$\hat{\vec{\theta}}$ は最良線形不偏推定量(BLUE)となるのが、最小二乗法の原理です。

正規方程式

最小二乗法の解 $\hat{\vec{\theta}}$ を求める方程式は以下となります。これを正規方程式といいます。

$$ {\bf X}^T{\bf X}\hat{\vec{\theta}}={\bf X}^T\vec{y} \Leftrightarrow \hat{\vec{\theta}}=({\bf X}^T{\bf X})^{-1}{\bf X}^T\vec{y} $$

ただし、上の通り ${\bf X }$ のランクが未知母数の次元数よりも小さくなる場合は、不定解となります。後ほど見ていきますが、1元配置の例がこれに当たります。

Pythonによる実装

NumPyによる実装を示します。

import numpy as np

繰り返し測定

デザイン行列 ${\bf X}$ は1行6列で、要素は全て1となります
$\hat{\mu}=37.68$ となりますが、これは $\bar{y}$ と一致します

# 9個のデータ
y = np.array([37.55, 37.74, 37.69, 37.71, 37.73, 37.81, 37.53, 37.58, 37.74])

# 1行9列のデザイン行列
X = np.array([np.ones(y.shape)]).T

theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
# array([37.67555556])

1元配置

1元配置の場合は上の式をそのまま使えません。というのも、全ての方程式において $\mu+\alpha_i$ の形で現れており、どう操作しても $\mu$ と $\alpha_i$ を分離できず、逆行列 ${\bf X}^T$ を計算できないのです(ランク落ち)。

そこで純粋な逆行列の代わりに、一般逆行列 ${\bf X}^{-}$ を採用した $\hat{\vec{\theta}}=({\bf X}^T{\bf X})^{-}{\bf X}^T\vec{y}$ を解きます。
一般逆行列については、こちら↓の投稿が参考になります。

https://www.slideshare.net/wosugi/ss-79624897

zellij.hatenablog.com

一般逆行列はnumpy.linalg.pinvで計算できます
デザイン行列は24行(r*a)5列(a+1)
- $\mu$ に対応する1列目は全て1
- $\alpha_1$ に対応する2列目は1〜6行目が1、 $\alpha_2$ に対応する3列目は7〜12行目が1, .... として、それ以外は0
結果は $\hat{\vec{\theta}}=(\mu, \alpha_1, \alpha_2, \alpha_3, \alpha_4)=(59, 13, 26, 17, 3)$

# 1元配置
y = np.array([
    64, 72, 68, 77, 56, 95,
    78, 91, 97, 82, 85, 77,
    75, 93, 78, 71, 63, 76,
    55, 66, 49, 64, 70, 68
]).T

# 24行5列
X = np.array([
    np.ones(24),
    np.concatenate([np.ones(6), np.zeros(18)]),
    np.concatenate([np.zeros(6), np.ones(6), np.zeros(12)]),
    np.concatenate([np.zeros(12), np.ones(6), np.zeros(6)]),
    np.concatenate([np.zeros(18), np.ones(6)])
]).T

# 行列式 np.linalg.det(X.T @ X) の値は0なので逆行列が存在しない
#theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y

# 逆行列の代わりに一般逆行列を使って推定する
theta_hat = np.linalg.pinv(X.T @ X) @ X.T @ y
# array([59., 13., 26., 17.,  3.])

多項式

デザイン行列は19行3列
- $\beta_0$ に対応する1列目は全て1
- $\beta_1$ に対応する2列目は $x$ , $\beta_2$ に対応する3列目は $x^2$ とする
結果は $\hat{\vec{\theta}}=(\beta_0, \beta_1, \beta_2)=(22.561, 1.668, -0.068)$

# 19個のデータ
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 2, 4, 6, 8, 10, 0, 3, 6, 9])
y = np.array([22.1, 24.5, 26.0, 26.8, 28.2, 28.9, 30.0, 30.4, 31.4, 21.9, 26.1, 28.5, 30.3, 31.5, 33.1, 22.8, 27.3, 29.8, 31.8])

# 19行3列のデザイン行列
X = np.array([np.ones(y.shape), x, x**2]).T

theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
# array([22.56123063,  1.66802044, -0.06795836])

まとめ

今回は線形モデルとそのパラメータを求める最小二乗法について、3つの例で見ていきました。またNumPyで正規方程式を解く実装を行いました。