線形回帰をPythonで数式から逃げずに実装してみた

こんにちは、エンジニアのBBです。

皆さんは線形回帰と聞いてどのようなものか明確に説明できるでしょうか。
なんとなくのイメージは出来るけど……と、明後日のほうを見ながら言う方を対象にちょっと前までそんな感じだった僕が説明していきます。

線形回帰とはなにか

まずはこの図を見てください。

あなたは左下から右上にかけてまっすぐに線を引きたくなったのではないでしょうか。こんな具合に。

このいい感じに線を引く操作を一定のルールに従って行うことを線形回帰といいます。いい感じに引かれた線は、未知の値（横軸）に対しておおよその答え（縦軸）がどのあたりに位置するかを予想してくれたり、変な値を見つけることに役立ちます。今回は線形回帰について数式を解きつつ、pythonのオープンソース機械学習ライブラリscikit-learn（通称：sklearn）に入っている線形回帰用のツール、LinearRegressionの偉大さを感じていきたいと思います。

機械学習はこの手順が基本！

さて、いきなり話がそれますが機械学習における基本的な手順を説明させてください。手順を把握しておくことで、テクニカルワードや数式の位置づけを体系的に理解する道しるべになると思います。機械学習の手順は大雑把に分けて３つです。

データを用意する
モデルを訓練する
モデルを評価する

この３つを繰り返してモデルの精度を改善していきます。上記手順のそれぞれにそれっぽい名前の手法や言い回しがあり、慣れないうちはすぐ迷子になるので意識的に自分がどの位置にいるのかを確認しておくといいと思います。

１．データを用意する

機械学習にはもちろん解析対象になるデータが必要になります。冒頭の図に使ったx, yのデータもそうですし、画像や文章も解析対象データとなります。解析対象のデータは必ず数値や数値の並び（行列やベクトルなど）に落とし込む必要があります。現実世界の現象をいかに数値に落とし込むかが機械学習の成功のカギとなり、データ活用の全行程の約8割をつぎ込む重要手順となっています。

２．モデルを訓練する

機械学習の名前の通り「機械」に手順１で用意したデータを「学習させる（訓練する）」工程です。モデルと一口に言っても様々あり、解析対象によっては自作する必要があるかもしれません。今回のテーマである線形回帰もこのモデルの一つとなります。

３．モデルを評価する

モデルで計算した結果（これもモデルによって形式はイロイロ）を評価する工程です。評価が悪かった場合は、手順１に戻ってデータを加工したり、手順２でモデルのチューニングやモデル自体を変更したりと頑張ります。

線形回帰でボストンの住宅価格データをみてみよう

ここからは数式の説明を交えつつ、線形回帰をpythonを使って実装していきます。また、線形回帰用ツール、LinearRegressionとの計算結果と比較して「やっぱ出来合いのツールスゲーわ」という感じで締めくくります。先述に合わせ、手順３つでこれからすることをざっくり説明します。

１．データを用意する
sklearnが用意してくれているサンプルデータ、ボストンの住宅価格データを使います。
データは、住宅価格とその他その住宅に関するステータス１３項目がひとまとめになっており、それが約５００軒分あります。イメージとしては横軸にステータス１３項目（一気に図示できませんが。。）、縦軸に住宅価格をとって線形回帰をする感じです。住宅のステータス１３項目を使って住宅価格を求めるとすると、ステータス１３項目を説明変数、住宅価格を目的変数といいます。

２．モデルを訓練する
モデルは宣言した通り、線形回帰を実装します。比較対象のツールはsklearnのLinearRegressionを使います。

３．モデルを評価する
結果の評価についてはRMSEをいう指標を使いますが、細かいことは後で記載します。

１．データを用意する

今回は出来合いのデータを用いるのでこのステップは飛ばしてステップ２から始めましょう。

２．モデルを訓練する

ここでは主に式を交えつつの説明をしていきますが、ここでの最終目標をまず見ておきましょう。

(1) $\begin{eqnarray*} \vec{\hat{y}}= X\vec{w} \end{eqnarray*}$

これです。このパラメータ $\vec{w}$ を求めることがゴールです。 $X$ は物件のステータス、 $\vec{\hat{y}}$ は予想された物件価格を表し、「ステータスから物件価格を予想する」式となります。このような関係を $\vec{\hat{y}}$ はパラメータ $\vec{w}$ に対して線形であるといいます。式(1)は行列、ベクトルが使われているので正直、慣れていないと何を言っているのかわからないと思いますので次の３ステップで考えていきましょう。

ステータス１項目、１軒分のデータを使う
ステータス１項目、全軒分のデータを使う
ステータス１３項目、全軒分のデータを使う

徐々に、データ構造を複雑にしていっているのがわかると思います。

１．ステータス１項目、１軒分のデータを使う

ステップ１では式(1)はこんな感じになります。

(2) $\begin{eqnarray*} \hat{y}= w_{0} + w_{1}x \end{eqnarray*}$

式(2)が意味するところは $\hat{y}$ （予想の住宅価格）は $x$ （ステータス１）に $w_{1}$ を掛けたものと定数 $w_{0}$ を足したものということです。これを図に起こすとこうなります。

そもそも点が１つしかなくこの点さえ通っていればなんでもよいので、いい感じの線は決められそうにありません。なので式(2)の出番はなさそうです。

２．ステータス１項目、全軒分のデータを使う

ステップ２では式(1)はこんな感じになります。

(3) $\begin{eqnarray*} \vec{\hat{y}}=w_{0} + w_{1}\vec{x} \end{eqnarray*}$

式に使われている文字の頭に矢印がっくついてがなにやら不穏な感じになっていますがこれはベクトルといいます。身構える必要はなく、単純に全軒分のデータに対してステップ１の式を書くのが面倒なのでこのように表記した単なる省略記法とかそんな感じにとらえてもらえるとわかりやすいと思います。実際はこんな感じの式が隠れているのです。

$\begin{eqnarray*} \begin{bmatrix} \hat{y}_{1} \\ \vdots \\ \hat{y}_{n}\end{bmatrix} = \begin{bmatrix}w_{0} + w_{1}x_{1} \\ \vdots \\w_{0} + w_{1}x_{n} \end{bmatrix}\nonumber \end{eqnarray}$

文字をベクトルとして表記するとこんな感じになります。

$\begin{eqnarray*} \vec{x} = \begin{bmatrix} x_{1} \\ \vdots \\ x_{506} \end{bmatrix}\nonumber , \vec{y} = \begin{bmatrix} y_{1} \\ \vdots \\ y_{506} \end{bmatrix}\nonumber \end{eqnarray}$

今回使うデータは506軒分ありますので $n = 506$ となります。506個の $x$ に対して、予想住宅価格 $\hat{y}$ も506個ありますのでベクトルとなります。ここまで省略すると欲が出てきます。式(3)には $w_{0}$ と $w_{1}$ という似たような表記が出てきています。まあ、こちらのさじ加減なのですが、この2つも下のようなベクトルを使って省略することができます。

$\begin{eqnarray*} \vec{w}=\binom{w_{0}}{w_{1}} \nonumber\end{eqnarray}$

省略するためにここで $\vec{x}$ にちょっと下に示すような細工をしておきます。

$\begin{eqnarray*} X = \begin{bmatrix} 1 & x_{1}\\ \vdots & \vdots \\ 1 & x_{506} \end{bmatrix} \nonumber\end{eqnarray}$

このようなベクトルのレベルアップ版を行列といいます。読んで字の通り行と列を持っています。この $X$ と $\vec{w}$ を使って式を変形するのですが、行列とベクトルの掛け算はちょっと複雑です。こちらの図がわかりやすかったので拝借します。

引用：Math03 行列と座標変換

この図にあるような順で掛け算が行われるため、追加した１を並べた列が $w_{0}$ を作ってくれます。
式(3)はこんな感じに変形されます。

(4) $\begin{eqnarray*} \begin{split} \vec{\hat{y}}&=\begin{bmatrix} w_{0}+w_{1}x_{1} \\ \vdots \\ w_{0}+w_{1}x_{506} \end{bmatrix}\\ &=\begin{bmatrix} 1 & x_{1}\\ \vdots & \vdots\\ 1 & x_{506} \end{bmatrix} \begin{bmatrix} w_{0} \\ w_{1} \end{bmatrix}\\ &= X\vec{w} \end{split} \end{eqnarray*}$

式の変形が完了しましたので今度は $\vec{w}$ を求める作業に入っていきましょう。 $\vec{w}$ を決定するには式(4)で引かれる線が最も「いい感じの線」になることが条件ですが、「いい感じの線」とはそもそも何でしょうか。「いい感じ」を考えるにあたり最小二乗法という手法を引っ張ってきます。最小二乗法を今回に当てはめると、 $\vec{w}$ をいろいろ変えて全ての「計測された住宅価格 $y$ 」と「予測された住宅価格 $\hat{y}$ 」の差を２乗した値の和が最も小さくなった時に「最もいい感じの線」が決まります。上記を式にするとこうなります。

(5) $\begin{eqnarray*} \begin{split} E&=\frac{1}{2n}\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}\\ &=\frac{1}{2n}\sum_{i=1}^{n}(y_{i}-\vec{w}^{T}\vec{x}_{i})^{2} \end{split} \end{eqnarray*}$

式(5)に出てきた $n$ は先述と同じくデータの数を示しています。506軒分のデータがありますので今回は $n=506$ となります。 $\sum_{i=1}^{n}$ は右にくっついてる式を $i=$ １～ $n$ まで足し続けなさいという意味です。シグマと読みます。式の右側はある数値の２乗なので常に正の値になります。ある数値は実数値の前提です。ちなみに機械学習の教科書等ではある数値 $a$ が実数であることを示すとき $a\in \mathbb{R}$ のように表記されることがあります。不意をついて現れるため慣れない人にとっては恐怖の対象だったかもしれませんが、こういった意味で使われているだけなのでもう怖がる必要はありません。

ここで興味があるのは式(5)ではなく、式(5)の左側が最も小さくなる $\vec{w}$ のみです。そんなもんわかるかいと思われるかもしれませんが、右側は $\vec{w}$ の２次関数になっているので求めることができます。２次関数はある $\vec{w}$ で最小値となるのですが、２次関数の形がわからない方はこちらのサイトが役に立ちます。このサイトの式の入力欄に「x^2」と入力するとよいでしょう。もし式(5)を分解して $\vec{w}$ の２乗項と $\vec{w}$ の項が混在しているのが不安で本当に最小値をとる $\vec{w}$ は一つに求まるのか僕のように確信が持てない方は上記サイトに「x^2+x」と入力し、xの係数を気が済むまで変更してみるとよいと思います。

$\vec{w}$ に関しての関数である式(5)が最小値をとることはわかりましたが、ではどのように求めるのかというと皆さん大好きな微分を敢行します。微分はある $\vec{w}$ でのグラフの傾きに相当しますので、最小値をとる $\vec{w}$ での微分は０になるわけです。式(5)の微分は以下の通りです。

(6) $\begin{eqnarray*} \begin{split} \frac{dE}{d\vec{w}}&=\frac{1}{d\vec{w}}\left[\frac{1}{2n}\left \{ \left ( y_{1}-\vec{w}^{T}\vec{x}_{1} \right )^{2}+\dots + \left ( y_{n}-\vec{w}^{T}\vec{x}_{n} \right )^{2}\right \}\right]\\ &=\frac{1}{2n}\left \{ 2\left ( y_{1}-\vec{w}^{T}\vec{x}_{1} \right )\left(- \vec{x}_{1}\right )+\dots + 2\left ( y_{n}-\vec{w}^{T}\vec{x}_{n} \right )\left(- \vec{x}_{i}\right )\right \}\\ &=-\frac{1}{n}\sum_{i=1}^{n}\left \{ \left( y_{i}-\vec{w}^{T}\vec{x}_{i}\right ) \vec{x}_{i} \right \} \\ &=-\frac{1}{n}X^{T}\left(\vec{y}-\vec{\hat{y}} \right ) \end{split} \end{eqnarray*}$

式(6)は $\vec{w}$ をちょっと動かしたときの変化量でもあり、グラフの傾きです。最小値をとる $\vec{w}$ の左では微分はマイナスになり、右側では微分はプラスになります。

なので適当に決めた $\vec{w}$ に対して式(6)を何度も何度も引き算していくと最小値をとる $\vec{w}$ に収束していく訳ですね。ちょうどボールの中にパチンコ玉を入れたかのように自分で移動していくイメージです。このとき使う式が式(7)です。

(7) $\begin{eqnarray*} \begin{split} \vec{w}&\leftarrow \vec{w}-\alpha \frac{dE}{d\vec{w}}\\ &=\vec{w}+\frac{\alpha}{n}X^{T}\left(\vec{y}-\vec{\hat{y}} \right ) \end{split} \end{eqnarray*}$

ここで突然現れた $\alpha$ ですが、 $\vec{w}$ をいかにちょっとずつ進めるかの度合いになります。 $\alpha$ を大きくすれば少ない計算量で $\vec{w}$ を収束できますが、答えの精度が落ちることが感覚的にわかりますね。長くなりましたがようやく $\vec{w}$ を求めるための式がわかりました。
あとはプログラムに落として計算させてみましょう。


#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression
from matplotlib import pyplot as plt

boston = load_boston()
np.random.seed(1234)

if __name__ == '__main__':
    # ###################### 下準備 #########################
    # 説明変数：13項目のうち1つだけ使う
    x = boston.data[:, 5]
    # 目的変数：住宅価格のデータを使う
    y = boston.target

    # ############## 計算式から回帰を実行 ####################
    # w0の項を追加するためにxに小細工をする
    large_x = np.array([[1, v] for v in x])
    # ランダムな値をw0, w1の初期値に設定する
    w = np.random.random(len(large_x[0]))
    # ちょっと変えるための度合いを決める
    alpha = 0.01
    # ------wをちょっとずつ変えていくループを実行------
    # ループ回数100回、1000回、10000回、100000回でwを保存する
    w_100 = []
    w_1000 = []
    w_10000 = []
    for i in range(100000):
        # 予測値y_hatを計算する。
        y_hat = np.dot(large_x, w)  # 式２－１
        # 損失関数の微分を計算する
        dw = np.dot((y - y_hat), large_x) / len(large_x)  # 式２－５
        # wをちょっとだけ動かす
        # print(i, 'w =', w)
        w += alpha * dw  # 式２－６
        if i == 100-1:
            # listはnp.copyしないとコピー元と同期してしまうので注意！
            w_100 = np.copy(w)
        elif i == 1000-1:
            w_1000 = np.copy(w)
        elif i == 10000-1:
            w_10000 = np.copy(w)

    # ######### skleanのツールLinearRegressionを使って回帰 ###########
    # ものすごく簡単...
    lr = LinearRegression(fit_intercept=True)
    lr.fit(large_x, y)
    pred = lr.predict(large_x)

    # ########################### 作図 #############################
    # ----------------散布図を作成-------------------
    plt.scatter(x, y, s=1)

    # ----------------計算式の結果-------------------
    plt.plot(x, np.dot(large_x, w_100), label='org_100')
    plt.plot(x, np.dot(large_x, w_1000), label='org_1000')
    plt.plot(x, np.dot(large_x, w_10000), label='org_10000')
    plt.plot(x, np.dot(large_x, w), c='black', linewidth=3.0, label='org_100000')

    # ---------------sklearnの結果------------------
    plt.plot(x, pred, c='y', linestyle='dashdot', label='sklearn')

    plt.legend()
    plt.show()

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

import numpy as np

from sklearn.datasets import load_boston

from sklearn.linear_model import LinearRegression

from matplotlib import pyplot as plt

boston = load_boston()

np.random.seed(1234)

if __name__ == '__main__':

# ###################### 下準備 #########################

# 説明変数：13項目のうち1つだけ使う

x = boston.data[:, 5]

# 目的変数：住宅価格のデータを使う

y = boston.target

# ############## 計算式から回帰を実行 ####################

# w0の項を追加するためにxに小細工をする

large_x = np.array([[1, v] for v in x])

# ランダムな値をw0, w1の初期値に設定する

w = np.random.random(len(large_x[0]))

# ちょっと変えるための度合いを決める

alpha = 0.01

# ------wをちょっとずつ変えていくループを実行------

# ループ回数100回、1000回、10000回、100000回でwを保存する

w_100 = []

w_1000 = []

w_10000 = []

for i in range(100000):

# 予測値y_hatを計算する。

y_hat = np.dot(large_x, w) # 式２－１

# 損失関数の微分を計算する

dw = np.dot((y - y_hat), large_x) / len(large_x) # 式２－５

# wをちょっとだけ動かす

# print(i, 'w =', w)

w += alpha * dw # 式２－６

if i == 100-1:

# listはnp.copyしないとコピー元と同期してしまうので注意！

w_100 = np.copy(w)

elif i == 1000-1:

w_1000 = np.copy(w)

elif i == 10000-1:

w_10000 = np.copy(w)

# ######### skleanのツールLinearRegressionを使って回帰 ###########

# ものすごく簡単...

lr = LinearRegression(fit_intercept=True)

lr.fit(large_x, y)

pred = lr.predict(large_x)

# ########################### 作図 #############################

# ----------------散布図を作成-------------------

plt.scatter(x, y, s=1)

# ----------------計算式の結果-------------------

plt.plot(x, np.dot(large_x, w_100), label='org_100')

plt.plot(x, np.dot(large_x, w_1000), label='org_1000')

plt.plot(x, np.dot(large_x, w_10000), label='org_10000')

plt.plot(x, np.dot(large_x, w), c='black', linewidth=3.0, label='org_100000')

# ---------------sklearnの結果------------------

plt.plot(x, pred, c='y', linestyle='dashdot', label='sklearn')

plt.legend()

plt.show()

↓実行結果

ちょっと動かす回数を変化させて図にしてみましたが、いかがでしょうか。
回数が多くなるにつれ少しずつ「いい感じの線」になってきていることが感覚的にもわかると思います。
最終的に１０００００回実行した結果はsklearnのツールを使った線と重なっていますので、よい結果が出ているのだと思います。

３．ステータス１３項目、全軒分のデータを使う

いよいよステップ３です。
ここまでお付き合いいただいた方はうんざりしているころだと思いますが、安心してください！
ここで使う式(1)は式(4)として既にステップ２で出てきてしまっています！もう式を追わなくてもよいのですよかった！
さて、使う式が同じだというのならステップ２とステップ３はどう違うのでしょうか。
答えは $\vec{w}$ と $X$ の中身になります。といっても $x$ の量に応じて $\vec{w}$ も増えるだけですが。

$\begin{eqnarray*} X=\begin{bmatrix} 1 &x_{11}&\dots &x_{1p} \\ \vdots &\vdots &\ddots &\vdots \\ 1 &x_{n1}&\dots &x_{np} \end{bmatrix}, \vec{w}=\begin{bmatrix} w_{0} \\ \vdots \\ w_{p} \end{bmatrix}\nonumber \nonumber\end{eqnarray}$

ここで新たな添え字 $p$ が登場しますがこれは $x$ に含まれるステータスの数になります。要するにこの場合は１３です。
ステップ３をプログラムに起こすとほぼステップ２と同じですが $X$ に値を格納する箇所が異なっていますのでコードを載せておきます。また、使用するステータスが増えましたので「いい感じの線」を作図によって目視することができませんので作図の箇所も不要です。


    # ###################### 下準備 #########################
    # 説明変数：13項目の全項目を使う
    # 　　　　　この時点で小細工をしてしまう
    large_x = np.array([np.concatenate(([1], v)) for v in boston.data])

# ###################### 下準備 #########################

# 説明変数：13項目の全項目を使う

# 　　　　　この時点で小細工をしてしまう

large_x = np.array([np.concatenate(([1], v)) for v in boston.data])

さて、残念なこと全項目を読み込んでプログラムを実行するとエラーが続発してしまい、うまく計算できなかったのではないでしょうか。
どうやら使用しているステータスのスケールががまちまちだった為、 $\vec{w}$ がうまく収束してくれなかった項目があったことが原因のようです。ここでいうスケールというのはあるステータスAは０～１に分布しているのに、ステータスBは１００～２００に分布していたりすることを指します。個別に回帰するのであれば $\alpha$ やループ回数を調整してやればよさそうですが、全部の項目に対応できる値となると見当もつきません。。

３．モデルを評価する

前章の最後の最後で手詰まりが起きてしまい、解析がうまくいきませんでしたね！この章では計算結果を評価するはずでしたが……細かいことは気にせずどうにか計算できるようにしてから結果を解析していきましょう。
ステップ３では項目間のスケールの違いがエラーの原因となっていましたのでどうにかしたいものです。こんな時に役に立つのが正規化です。
ここでの正規化というのはデータの平均値からのばらつき（偏差）を、ばらつきの平均値（標準偏差）で割ってやることで、項目間のスケールが違っても平等に評価できるようにする処理を指します。標準化とも呼ばれ、式でいうと以下のようになります。

(8) $\begin{eqnarray*} {x}'=\frac{x-\bar{x}}{\sigma } \end{eqnarray*}$

ここで $\bar{x}$ を $x$ の平均値、 $\sigma$ を標準偏差といいます。この処理を全項目について適応するには以下のコードを追加します。


    # ###################### 下準備 #########################
    # 説明変数：13項目の全項目を使う
    # 　　　　　この時点で小細工をしてしまう
    large_x = np.array([np.concatenate(([1], v)) for v in boston.data])

    # 目的変数：住宅価格のデータを使う
    y = boston.target

    # ---------------- 説明変数を標準化する -----------------
    large_x_ = np.mean(large_x, 0)
    large_xs = np.std(large_x, 0)
    for i in range(1, len(large_x[0])):
        large_x[:, i] = (large_x[:, i] - large_x_[i]) / large_xs[i]

# ###################### 下準備 #########################

# 説明変数：13項目の全項目を使う

# 　　　　　この時点で小細工をしてしまう

large_x = np.array([np.concatenate(([1], v)) for v in boston.data])

# 目的変数：住宅価格のデータを使う

y = boston.target

# ---------------- 説明変数を標準化する -----------------

large_x_ = np.mean(large_x, 0)

large_xs = np.std(large_x, 0)

for i in range(1, len(large_x[0])):

large_x[:, i] = (large_x[:, i] - large_x_[i]) / large_xs[i]

これで何とか無事、 $\vec{w}$ を求めることができましたので、さっそく結果の評価をしていきましょう。今回はRMSE（Root Mean Square Error）：平均二乗誤差を使って結果を評価していきます。比較対象はライブラリsklearnのツールを使って計算した結果のRMSEです。RMSEは以下の式で計算できます。

(9) $\begin{eqnarray*} RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(\hat{y}_{i}-y_{i} \right )^{2}} \end{eqnarray*}$

RMSEは予想した住宅価格と本来の住宅価格の差に関連する値のようですね。予測と実際の値に差異が少なければ０に近づき、差異が大きいほど０から離れていきます。そのためRMSEが小さいほど「いい感じの線」であるといえます。では早速プログラムを使って計算してみましょう！


import numpy as np
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression

boston = load_boston()
np.random.seed(1234)

if __name__ == '__main__':
    # ###################### 下準備 #########################
    # 説明変数：13項目の全項目を使う
    # 　　　　　この時点で小細工をしてしまう
    large_x = np.array([np.concatenate(([1], v)) for v in boston.data])

    # 目的変数：住宅価格のデータを使う
    y = boston.target

    # ---------------- 説明変数を標準化する -----------------
    large_x_ = np.mean(large_x, 0)
    large_xs = np.std(large_x, 0)
    for i in range(1, len(large_x[0])):
        large_x[:, i] = (large_x[:, i] - large_x_[i]) / large_xs[i]

    # ############## 計算式から回帰を実行 ####################
    # ランダムな値をw0, w1の初期値に設定する
    w = np.random.random(len(large_x[0]))
    # ちょっと変えるための度合いを決める
    alpha = 0.01
    # ------wをちょっとずつ変えていくループを実行------
    # ループ回数100回、1000回、10000回、100000回でy_hatを保存する
    y_hat =[]
    y_hat_100 = []
    y_hat_1000 = []
    y_hat_10000 = []
    for i in range(100000):
        # 予測値y_hatを計算する。
        y_hat = np.dot(large_x, w)  # 式２－１
        # 損失関数の微分を計算する
        dw = np.dot((y - y_hat), large_x) / len(large_x)  # 式２－５
        # wをちょっとだけ動かす
        # print(i, 'w =', w)
        w += alpha * dw  # 式２－６
        if i == 100-1:
            # listはnp.copyしないとコピー元と同期してしまうので注意！
            y_hat_100 = np.copy(y_hat)
        elif i == 1000-1:
            y_hat_1000 = np.copy(y_hat)
        elif i == 10000-1:
            y_hat_10000 = np.copy(y_hat)

    # ######### skleanのツールLinearRegressionを使って回帰 ###########
    # ものすごく簡単...
    lr = LinearRegression(fit_intercept=True)
    lr.fit(large_x, y)
    pred = lr.predict(large_x)

    # ######################## RMSEを計算する ########################
    # -------------- 計算式の結果のループ回数ごとのRMSE --------------
    err_100 = y_hat_100 - y
    rmse_100 = np.dot(err_100, err_100)/len(y)
    err_1000 = y_hat_1000 - y
    rmse_1000 = np.dot(err_1000, err_1000) / len(y)
    err_10000 = y_hat_10000 - y
    rmse_10000 = np.dot(err_10000, err_10000) / len(y)
    err_100000 = y_hat - y
    rmse_100000 = np.dot(err_100000, err_100000) / len(y)
    print('rmse_100:', rmse_100)
    print('rmse_1000:', rmse_1000)
    print('rmse_10000:', rmse_10000)
    print('rmse_100000:', rmse_100000)

    # ----------------- sklearnの結果で計算したRMSE ------------------
    err_lr = pred - y
    rmse_lr = np.dot(err_lr, err_lr) / len(y)
    print('rmse_lr:', rmse_lr)

    # 出力結果はこちら
    # rmse_100: 96.77816336037456
    # rmse_1000: 22.148644640929934
    # rmse_10000: 21.897780844972825
    # rmse_100000: 21.897779217687496
    # rmse_lr: 21.897779217687496

import numpy as np

from sklearn.datasets import load_boston

from sklearn.linear_model import LinearRegression

boston = load_boston()

np.random.seed(1234)

if __name__ == '__main__':

# ###################### 下準備 #########################

# 説明変数：13項目の全項目を使う

# 　　　　　この時点で小細工をしてしまう

large_x = np.array([np.concatenate(([1], v)) for v in boston.data])

# 目的変数：住宅価格のデータを使う

y = boston.target

# ---------------- 説明変数を標準化する -----------------

large_x_ = np.mean(large_x, 0)

large_xs = np.std(large_x, 0)

for i in range(1, len(large_x[0])):

large_x[:, i] = (large_x[:, i] - large_x_[i]) / large_xs[i]

# ############## 計算式から回帰を実行 ####################

# ランダムな値をw0, w1の初期値に設定する

w = np.random.random(len(large_x[0]))

# ちょっと変えるための度合いを決める

alpha = 0.01

# ------wをちょっとずつ変えていくループを実行------

# ループ回数100回、1000回、10000回、100000回でy_hatを保存する

y_hat =[]

y_hat_100 = []

y_hat_1000 = []

y_hat_10000 = []

for i in range(100000):

# 予測値y_hatを計算する。

y_hat = np.dot(large_x, w) # 式２－１

# 損失関数の微分を計算する

dw = np.dot((y - y_hat), large_x) / len(large_x) # 式２－５

# wをちょっとだけ動かす

# print(i, 'w =', w)

w += alpha * dw # 式２－６

if i == 100-1:

# listはnp.copyしないとコピー元と同期してしまうので注意！

y_hat_100 = np.copy(y_hat)

elif i == 1000-1:

y_hat_1000 = np.copy(y_hat)

elif i == 10000-1:

y_hat_10000 = np.copy(y_hat)

# ######### skleanのツールLinearRegressionを使って回帰 ###########

# ものすごく簡単...

lr = LinearRegression(fit_intercept=True)

lr.fit(large_x, y)

pred = lr.predict(large_x)

# ######################## RMSEを計算する ########################

# -------------- 計算式の結果のループ回数ごとのRMSE --------------

err_100 = y_hat_100 - y

rmse_100 = np.dot(err_100, err_100)/len(y)

err_1000 = y_hat_1000 - y

rmse_1000 = np.dot(err_1000, err_1000) / len(y)

err_10000 = y_hat_10000 - y

rmse_10000 = np.dot(err_10000, err_10000) / len(y)

err_100000 = y_hat - y

rmse_100000 = np.dot(err_100000, err_100000) / len(y)

print('rmse_100:', rmse_100)

print('rmse_1000:', rmse_1000)

print('rmse_10000:', rmse_10000)

print('rmse_100000:', rmse_100000)

# ----------------- sklearnの結果で計算したRMSE ------------------

err_lr = pred - y

rmse_lr = np.dot(err_lr, err_lr) / len(y)

print('rmse_lr:', rmse_lr)

# 出力結果はこちら

# rmse_100: 96.77816336037456

# rmse_1000: 22.148644640929934

# rmse_10000: 21.897780844972825

# rmse_100000: 21.897779217687496

# rmse_lr: 21.897779217687496

ループ回数が多ければ多いほど徐々に値が小さくなっていき、１００００回を超えたあたりからほぼ値が変わらなくなっています。 $\vec{w}$ が収束していっていることが実感できますね。さらに１０００００回の結果はsklearnで計算した結果とぴたりと同じになっています！
本来はこの結果をもとに良い悪いの評価を行い、データの変更や、モデルの調整を行います。ちょうどこの章の冒頭で行った正規化がデータの変更に対応しています。

まとめ

今回は線形回帰についてまとめてきました。個人的には計算を繰り返すことで $\vec{w}$ が収束していく様子が実感できたのが非常に面白かったです。この $\vec{w}$ が収束する過程が機械学習の学習の部分です。sklearnに代表されるツールは便利ですが、実際に数式をコードに落とし込む工程は勉強になることが多く苦労する価値はあると感じています。
軽い気持ちで書いてきた今回の記事ですが、長文になってしまいました。ここまでお付き合いいただいた方はありがとうございます。線形回帰で躓いている方の助けになれば幸いです！

TOEICの呪縛

Variational Auto Encoder 〜外れ値検知への応用〜

線形回帰をPythonで数式から逃げずに実装してみた

線形回帰とはなにか

機械学習はこの手順が基本！

１．データを用意する

２．モデルを訓練する

３．モデルを評価する

線形回帰でボストンの住宅価格データをみてみよう

１．データを用意する

２．モデルを訓練する

１．ステータス１項目、１軒分のデータを使う

２．ステータス１項目、全軒分のデータを使う

３．ステータス１３項目、全軒分のデータを使う

３．モデルを評価する

まとめ

最近の記事

顔認証（DeepFace vs GPT4o）

時系列解析再訪

πの新しい公式

中心極限定理

DBSCANを用いたクラスタリング

レイアウト解析

Learning Chainer with Examples 〜深層学習への入門〜

行列の条件数

LSTM

Bayes推定

Auto-Kerasで深層学習のモデルを自動で生成する

時系列解析再訪

Google Vision APIでOCR

アーカイブ

カテゴリー