メインコンテンツへスキップ
見出し画像

＜学習シリーズ＞線形モデル編：単回帰分析

    KIYO

    １．概要

    １－１．緒言

    　本記事は”学習シリーズ”として自分の勉強備忘録用になります。

    　過去の記事で機械学習・AIの記事を多数作成しましたが、シンプルな線形モデルは外挿が比較的得意のためいろんな分野で使用されます。
    　本記事では一から線形モデルを学びなおしたいと思います。

    　今後記載していくモデルは下記の通りです。

    • 重回帰

    • ロジスティック回帰

    • Ridge回帰

    • Rasso回帰

    • ガウス過程回帰

    １－２．用語・記号の説明（全般）

    　本記事で使用する用語および記号は下記の通りです。

    【用語一覧】

    ●回帰 (regression)：実数値を予測する問題
    ●分類 (classification)：    カテゴリ（離散値）を予測する問題
    ●教師あり学習 (supervised learning)：機械学習でデータセットに対して正解値（ラベル）がある学習方法
    ●教師なし学習 (unsupervised learning)：    機械学習でデータセットに対して正解値（ラベル）が無い学習方法
    ●予測値 (predicted value)：    関数で計算された出力変数y（用語は下記参照）
    ●目標値 (target value)：    予測値がとるべき値（教師あり学習の正解値）
    ●    目的関数 (objective function)：機械学習において性能の良さの指標を表す関数です。一般的には”予想値と目標値の差異”から作成される関数であり、この関数を最小化することでよいモデルと判断します。
    ●多重共線性（multicollinearity）：    多変量解析（重回帰分析など）において、いくつかの説明変数間で線形関係（強い相関関係）があると共線性といい、共線性が複数認められる場合は多重共線性があると言います。
    ●二乗和誤差 (sum-of-squares error)：正解値yyと予測値y^\hat{y}の差の２乗(y−y^)2(y-\hat{y})^2です。これをモデルと正解値の誤差とも呼びます。
    ●    最小二乗法：二乗和誤差を最小化することでモデルの当てはまりを最適化する手法です。

    【記号一覧】

    ●f()f()：y=f(x)y=f(x)におけるf()であり関数と呼びます。中身は入力変数xを処理して新しい出力変数yを生成する計算式です。
    ●xx：y=f(x)y=f(x)におけるxです。複数は下記の通り複数あります。
    　名称１：説明変数（Explanatory variable）
    　名称２：独立変数（Independent variable）
    　名称３：外生変数（Exogenous variable）
    　名称４：入力変数（Input variable）
    　名称５：入力値（Input value）
    ●yy：y=f(x)y=f(x)におけるyです。複数は下記の通り複数あります。
    　名称１：目的変数（Response variable）
    　名称２：従属変数（Dependent variable）
    　名称３：被説明変数（Explained variable）
    　名称４：内生変数（endogenous variable）
    　名称５：出力変数（Output variable）
    　名称６：予測値/推論値（prediction value）
    ●wiw_i(weight)：重み（単回帰の場合は傾きとも言う）
    ●bb(bias)：バイアス（単回帰の場合は切片とも言う）

    １－３．サンプル用データ

    　サンプルデータはScikit-learnのToy datasetsを使用します。今回は回帰(Regression)のため"diabetes"を使用しました。単回帰は変数が1つだけの1次元モデルのため、特徴量として重要な"s5"だけ抽出しました。

    [IN]
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

#データセットの読み込み
diabetes = datasets.load_diabetes()
data, target = diabetes.data, diabetes.target
df_data = pd.DataFrame(data, columns=diabetes.feature_names)
df_data = df_data[['s5']] #単回帰として使用するため1次元のみ抽出
print(type(df_data), df_data.shape, type(target), target.shape)

#データの可視化
fig, ax = plt.subplots(1, 1, figsize=(10, 5), facecolor='w')
ax.scatter(df_data, target)
ax.set(xlabel='s5', ylabel='target', title='Correlation of Diabetes')
plt.grid()
plt.show()

[OUT]
<class 'pandas.core.frame.DataFrame'> (442, 1) <class 'numpy.ndarray'> (442,)
    画像

    ２．単回帰（線形モデル）とは：y=ax+b

    　単回帰分析とは 1 つの入力変数xxから 1 つの出力変数yyを推測する処理です。入力値に重みをかけバイアスを足したものを回帰方程式、グラフを回帰直線と呼びます。モデルの形が直線になることから線形モデルとも呼ばれます。

    単回帰：y=ax+b単回帰：y = ax+b

    ●xx：入力変数
    ●yy(またはy^\hat y)：出力変数
    ●aa：傾き(機械学習では重みと呼びwwを使う）
    ●bb：切片（またはバイアス）

    　推定値y^\hat yとデータ値yiy_iの誤差を最小にすることで最適なパラメータaa, bbを決定できます。

    ３．直接解法の計算

    　深層学習（ディープラーニング）などではモデルの学習時に誤差逆伝搬を使用して学習させますが、単回帰の重みとバイアスは直接計算可能です。
    　最小二乗法を用いて誤差が最小になるパラメータを直接計算してみます。

    ３－１．結論

    　求めたいパラメータ(傾きaa, 切片bb)は下記で計算できます。

    a=1n∑i=1nxiyi−yˉxˉ1n∑i=1nxi2−xˉ2=Cov(X,Y)Var(X)a =\frac{\frac{1}{n}\sum_{i=1}^{n} x_iy_i - \bar{y}\bar{x}}{\frac{1}{n}\sum_{i=1}^{n} x_i^2 - \bar{x}^2}=\frac{Cov(X,Y)}{Var(X)}

    b=yˉ−axˉb=\bar{y}-a\bar{x}

    　参考までに展開した状態だと下記の通りです。

    a=n∑i=1nxiyi−∑i=1nxi∑i=1nyin∑i=1nxi2−(∑i=1nxi)2b=(∑i=1nxi2)∑i=1nyi−(∑i=1nxi)∑i=1nxiyin∑i=1nxi2−(∑i=1nxi)2\begin{aligned} a &= \frac{n \sum_{i=1}^{n} x_iy_i - \sum_{i=1}^{n} x_i \sum_{i=1}^{n} y_i}{n \sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2} \\ b &= \frac{ (\sum_{i=1}^{n} x_i^2)\sum_{i=1}^{n} y_i - (\sum_{i=1}^{n} x_i)\sum_{i=1}^{n} x_iy_i }{n \sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2} \end{aligned}

    ３－２．用語・数式の一覧

    　まずは各種変数を用語・数式を示します。統計の用語は下記に詳細を記述しております。

    入力変数x= (x1 x2 ⋯ xn) 入力変数x = \begin{pmatrix}x_1 \ x_2 \ \cdots \ x_n\end{pmatrix} 

    出力変数y= (y1 y2 ⋯ yn)出力変数y = \begin{pmatrix}y_1 \ y_2 \ \cdots \ y_n\end{pmatrix}

    　３－２－１．平均値

    　平均値とは「全データを合計値をデータ数Nで割ったもの」でありデータの重心を意味します。下記より、nxˉ=∑i=1nxi n\bar{x}=\sum_{i=1}^{n} x_i (要素数×平均値=各要素の合計値)への変換が可能です。

    xの平均値：xˉ=1n∑i=1nxi=1n(x1+x2+⋯+xn)xの平均値：\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i = \frac{1}{n} \begin{pmatrix} x_1 + x_2 + \cdots + x_n \end{pmatrix}

    yの平均値：yˉ=1n∑i=1nyi=1n(y1+y2+⋯+yn)yの平均値：\bar{y} = \frac{1}{n} \sum_{i=1}^{n} y_i = \frac{1}{n} \begin{pmatrix} y_1 + y_2 + \cdots + y_n \end{pmatrix}

    　３－２－２．分散

    　分散とは「平均値から各データの差の2乗を全て足した値」でありデータのばらつきを表します（下記式は標本分散となります）。

    xの分散σx2=1n∑i=1n(xi−xˉ)2=1n∑i=1nxi2−xˉ2xの分散\sigma_x^2 = \frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})^2=\frac{1}{n}\sum_{i=1}^{n} x_i^2 - \bar{x}^2

    　右辺の式導出は下記の通りです。

    σx2の式展開=1n∑i=1n(xi−xˉ)2=1n∑i=1n(xi2−2xixˉ+xˉ2)\sigma_x^2 の式展開= \frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})^2= \frac{1}{n}\sum_{i=1}^{n} (x_i^2 - 2x_i \bar{x} + \bar{x}^2)

    =1n(∑i=1nxi2−2xˉ∑i=1nxi+ ∑i=1nxˉ2)= \frac{1}{n}(\sum_{i=1}^{n} x_i^2 - 2\bar{x}\sum_{i=1}^{n} x_i + \sum_{i=1}^{n} \bar{x}^2)

    =1n(∑i=1nxi2−2nxˉ2+nxˉ2)=1n∑i=1nxi2−xˉ2= \frac{1}{n}(\sum_{i=1}^{n} x_i^2 - 2n\bar{x}^2 + n \bar{x}^2)=\frac{1}{n}\sum_{i=1}^{n} x_i^2 - \bar{x}^2

    　３－２－３．共分散

    　共分散とは「二組の対応するデータの間の関係を表す数値」であり、X の偏差 × Y の偏差の平均=(x-平均のx)(y-平均のy)の平均で計算できます。

    画像

    共分散行列Cov(x,y)=1n∑i=1n(xi−xˉ)(yi−yˉ)=1n∑i=1nxiyi−xˉyˉ共分散行列Cov(x,y) = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})= \frac{1}{n}\sum_{i=1}^{n} x_iy_i - \bar{x}\bar{y}

    Covの展開=1n((x1−xˉ)(y1−yˉ)+(x2−xˉ)(y2−yˉ)+⋯+(xn−xˉ)(yn−yˉ))Covの展開=\frac{1}{n} \begin{pmatrix} (x_1-\bar{x})(y_1-\bar{y}) + (x_2-\bar{x})(y_2-\bar{y}) + \cdots + (x_n-\bar{x})(y_n-\bar{y}) \end{pmatrix}

    右辺式の導出：1n∑i=1n(xi−xˉ)(yi−yˉ)=1n∑i=1n(xiyi−xiyˉ−xˉyi+xˉyˉ)右辺式の導出：\frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})= \frac{1}{n} \sum_{i=1}^{n} (x_iy_i - x_i \bar{y} - \bar{x}y_i + \bar{x}\bar{y})

    =1n(∑i=1nxiyi−yˉ∑i=1nxi−xˉ∑i=1nyi+nxˉyˉ)=1n(∑i=1nxiyi−yˉ(nxˉ)−xˉ(nyˉ)+nxˉyˉ)= \frac{1}{n}(\sum_{i=1}^{n} x_iy_i - \bar{y}\sum_{i=1}^{n} x_i - \bar{x}\sum_{i=1}^{n} y_i + n\bar{x}\bar{y})=\frac{1}{n}(\sum_{i=1}^{n} x_iy_i - \bar{y}(n\bar{x})- \bar{x}(n\bar{y}) + n\bar{x}\bar{y})

    =1n(∑i=1nxiyi−nxˉyˉ)=1n∑i=1nxiyi−xˉyˉ= \frac{1}{n}(\sum_{i=1}^{n} x_iy_i - n\bar{x}\bar{y})=\frac{1}{n}\sum_{i=1}^{n} x_iy_i - \bar{x}\bar{y}

    ３－３．重みとバイアスの導出

    　３－３－１：最小二乗和Eの計算

    　まずは最小化したい二乗和誤差を計算します。

    二乗和誤差E=∑i=1n(yi−y^)=∑i=1n(yi−(axi+b))2二乗和誤差E=\sum_{i=1}^{n}(y_i-\hat{y})=\sum_{i=1}^{n} (y_i - (ax_i + b))^2

    　求めたい傾きaと切片bを求めるには、aとbそれぞれで二乗和誤差Eが最小になる値、つまりEに対する偏微分が0になる値∂E∂a=0\frac{\partial E}{\partial a}=0, ∂E∂b=0\frac{\partial E}{\partial b}=0を求めます。

    (yi−(axi+b))2(y_i - (ax_i + b))^2 を分解してaとbそれぞれでまとめると下記の通りです。

    E=(yi−(axi+b))2=yi2−2yi(axi+b)+(axi+b)2E=(y_i - (ax_i + b))^2= y_i^2 - 2y_i(ax_i + b) + (ax_i + b)^2

    =yi2−2yi(axi+b)+a2xi2+2abxi+b2=y_i^2 - 2y_i(ax_i + b) + a^2 x_i^2 + 2abx_i + b^2

    　３－３－２．aでの偏微分を計算

    　二乗和誤差Eをaでまとめると下記の通りです。

    E=∑i=1n(xi2a2+(2bxi−2xiyi)a+(b2−2byi+yi2))E= \sum_{i=1}^{n} (x_i^2a^2+(2bx_i-2x_iy_i)a+ (b^2 - 2by_i + y_i^2))

    　上記をaで偏微分すると下記の通りです。

    ∂E∂a=∑i=1n(2axi2+2(bxi−xiyi))=2a∑i=1nxi2+2∑i=1n(bxi−xiyi)=2a∑i=1nxi2+2b∑i=1nxi−2∑i=1nxiyi\begin{aligned} \frac{\partial E}{\partial a} = \sum_{i=1}^{n} (2a x_i^2 + 2(bx_i - x_iy_i)) \\ =2a \sum_{i=1}^{n} x_i^2 + 2 \sum_{i=1}^{n} (bx_i - x_iy_i) \\ =2a \sum_{i=1}^{n} x_i^2 +2b\sum_{i=1}^{n}x_i-2\sum_{i=1}^{n} x_iy_i \end{aligned}

    　∂E∂a=0\frac{\partial E}{\partial a}=0より、aで整理すると下記の通りです。

    a=∑i=1nxiyi−b∑i=1nxi∑i=1nxi2a = \frac{\sum_{i=1}^{n} x_iy_i - b\sum_{i=1}^{n} x_i}{\sum_{i=1}^{n} x_i^2}

    　y=ax+by=ax+bより上式にb=y−axb=y-axを代入するとbが消せます。

    a=∑i=1nxiyi−(yˉ−axˉ)∑i=1nxi∑i=1nxi2a = \frac{\sum_{i=1}^{n} x_iy_i - (\bar{y} - a\bar{x})\sum_{i=1}^{n} x_i}{\sum_{i=1}^{n} x_i^2}

    a∑i=1nxi2=∑i=1nxiyi−yˉ∑i=1nxi+axˉ∑i=1nxia\sum_{i=1}^{n} x_i^2 = \sum_{i=1}^{n} x_iy_i - \bar{y}\sum_{i=1}^{n} x_i + a\bar{x}\sum_{i=1}^{n} x_i

    a(∑i=1nxi2)=∑i=1nxiyi−yˉ(∑i=1nxi)+axˉ(∑i=1nxi)a\left(\sum_{i=1}^{n} x_i^2\right) = \sum_{i=1}^{n} x_iy_i - \bar{y}\left(\sum_{i=1}^{n} x_i\right) + a\bar{x}\left(\sum_{i=1}^{n} x_i\right)

    　aで整理すると

    a=∑i=1nxiyi−yˉ∑i=1nxi∑i=1nxi2−xˉ∑i=1nxia = \frac{\sum_{i=1}^{n} x_iy_i - \bar{y}\sum_{i=1}^{n} x_i}{\sum_{i=1}^{n} x_i^2 - \bar{x}\sum_{i=1}^{n} x_i}

    a=∑i=1nxiyi−yˉnxˉ∑i=1nxi2−nxˉ2a = \frac{\sum_{i=1}^{n} x_iy_i - \bar{y}n\bar{x}}{\sum_{i=1}^{n} x_i^2 - n\bar{x}^2}

    a=1n∑i=1nxiyi−yˉxˉ1n∑i=1nxi2−xˉ2=Cov(X,Y)Var(X)a = \frac{\frac{1}{n}\sum_{i=1}^{n} x_iy_i - \bar{y}\bar{x}}{\frac{1}{n}\sum_{i=1}^{n} x_i^2 - \bar{x}^2}=\frac{Cov(X,Y)}{Var(X)}

    　３－３－３．bでの偏微分を計算

    　二乗和誤差Eをbでまとめると下記の通りです。

    E=∑i=1n(b2+(−2yi2axi)b+(yi2−2axiyi+a2xi2))E= \sum_{i=1}^{n}(b^2+(-2y_i2ax_i)b + (y_i^2 - 2ax_i y_i+ a^2 x_i^2 ))

    　上記をbで偏微分すると下記の通りです。

    ∂E∂b=∑i=1n(2b−2yi+2axi)=2∑i=1nb−2∑i=1nyi+2a∑i=1nxi\begin{aligned} \frac{\partial E}{\partial b} = \sum_{i=1}^{n} (2b - 2y_i + 2ax_i) \\ = 2 \sum_{i=1}^{n} b - 2 \sum_{i=1}^{n} y_i + 2a \sum_{i=1}^{n} x_i \end{aligned}

    =2nb−2∑i=1nyi+2a∑i=1nxi=2nb- 2 \sum_{i=1}^{n} y_i + 2a \sum_{i=1}^{n} x_i

    　∂E∂b=0\frac{\partial E}{\partial b}=0のためbで整理すると、

    b=1n∑i=1nyi−1na∑i=1nxi=yˉ−axˉb =\frac{1}{n} \sum_{i=1}^{n} y_i - \frac{1}{n}a \sum_{i=1}^{n} x_i=\bar{y}-a\bar{x}

    ３－４．コラム：重みとバイアスの導出（行列編）

    　参考として上記の正規方程式を行列でも解いてみます。

    ∂E∂a=∑i=1n(2axi2+2(bxi−xiyi))=2a∑i=1nxi2+2b∑i=1nxi−2∑i=1nxiyi=2(a∑i=1nxi2+b∑i=1nxi−∑i=1nxiyi)\begin{aligned}\frac{\partial E}{\partial a} = \sum_{i=1}^{n} (2a x_i^2 + 2(bx_i - x_iy_i)) \\=2a \sum_{i=1}^{n} x_i^2 +2b\sum_{i=1}^{n}x_i-2\sum_{i=1}^{n} x_iy_i \\=2(a \sum_{i=1}^{n} x_i^2 +b\sum_{i=1}^{n}x_i-\sum_{i=1}^{n} x_iy_i) \end{aligned}

    ∂E∂b=∑i=1n(2b−2yi+2axi)=2∑i=1nb−2∑i=1nyi+2a∑i=1nxi=2(∑i=1nb−∑i=1nyi+a∑i=1nxi)\begin{aligned}\frac{\partial E}{\partial b} = \sum_{i=1}^{n} (2b - 2y_i + 2ax_i) \\= 2 \sum_{i=1}^{n} b - 2 \sum_{i=1}^{n} y_i + 2a \sum_{i=1}^{n} x_i \\= 2(\sum_{i=1}^{n} b - \sum_{i=1}^{n} y_i + a \sum_{i=1}^{n} x_i) \end{aligned}

    ∂E∂a=0\frac{\partial E}{\partial a}=0, ∂E∂b=0\frac{\partial E}{\partial b}=0より

    a∑i=1nxi2+b∑i=1nxi=∑i=1nxiyia \sum_{i=1}^{n} x_i^2 +b\sum_{i=1}^{n}x_i=\sum_{i=1}^{n} x_iy_i

    a∑i=1nxi+∑i=1nb=∑i=1nyia \sum_{i=1}^{n} x_i + \sum_{i=1}^{n} b =\sum_{i=1}^{n} y_i

    　これを行列で記載すると単回帰の正規方程式の行列が得られます。

    [∑i=1nxi2∑i=1nxi∑i=1nxin][ab]=[∑i=1nxiyi∑i=1nyi]\begin{bmatrix} \sum_{i=1}^{n} x_i^2 & \sum_{i=1}^{n} x_i \\ \sum_{i=1}^{n} x_i & n \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^{n} x_iy_i \\ \sum_{i=1}^{n} y_i \end{bmatrix}

    　両辺に左辺にある逆行列をかけるとパラメータが得られます。

    [ab]=1(n∑i=1nxi2−(∑i=1nxi)2)[n−∑i=1nxi−∑i=1nxi∑i=1nxi2][∑i=1nxiyi∑i=1nyi]\begin{bmatrix} a \\ b \end{bmatrix} = \frac{1}{\left(n \sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2\right)} \begin{bmatrix} n & -\sum_{i=1}^{n} x_i \\ -\sum_{i=1}^{n} x_i & \sum_{i=1}^{n} x_i^2 \end{bmatrix} \begin{bmatrix} \sum_{i=1}^{n} x_iy_i \\ \sum_{i=1}^{n} y_i \end{bmatrix}

    　上記計算より、連立方程式で求めた解と同じになることが確認できます。

    [ab]=1(n∑i=1nxi2−(∑i=1nxi)2)[n∑i=1nxiyi−∑i=1nxi∑i=1nyi∑i=1nxi2∑i=1nyi−∑i=1nxi∑i=1nxiyi]\begin{bmatrix} a \\b \end{bmatrix}= \frac{1}{\left(n \sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2\right)} \begin{bmatrix} n\sum_{i=1}^{n} x_iy_i - \sum_{i=1}^{n} x_i \sum_{i=1}^{n} y_i \\ \sum_{i=1}^{n} x_i^2 \sum_{i=1}^{n} y_i-\sum_{i=1}^{n} x_i \sum_{i=1}^{n} x_iy_i \end{bmatrix}

    【参考：逆行列の計算】

    [ABCD]−1=1AD−BC[D−B−CA]\begin{bmatrix} A & B \\ C & D \end{bmatrix}^{-1} = \frac{1}{AD - BC} \begin{bmatrix} D & -B \\ -C & A \end{bmatrix}

    ４．単回帰モデルの実装：スクラッチ編

    　動作を理解するためにPythonのライブラリを使用せずスクラッチで計算します。目視で確認できるようにデータ数を1桁に調整しました。

    [IN]
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sns
import os

from sklearn import datasets
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

#データセットの読み込み
diabetes = datasets.load_diabetes()
data, target = diabetes.data, diabetes.target
df_data = pd.DataFrame(data, columns=diabetes.feature_names)
df_data = df_data[['s5']] #単回帰として使用するため1次元のみ抽出

#データの分割：参照用のデータ数を少なくする
x_train, x_test, y_train, y_test = train_test_split(df_data, target, test_size=0.98, random_state=0)
print(f'x_train: {x_train.shape}, x_test: {x_test.shape}, y_train: {y_train.shape}, y_test: {y_test.shape}')

#Excel処理用としてデータを結合／出力
array_dataset = np.concatenate([x_train, y_train.reshape(-1, 1)], axis=1)
df_dataset = pd.DataFrame(array_dataset, columns=['x', 'y'])
display(df_dataset)

#Excelファイルの出力
if not os.path.exists('data'):
    os.mkdir('data') #ディレクトリが存在しない場合のみ作成
    print(f'ディレクトリを作成しました: {os.path.abspath("data")}')

df_dataset.to_excel('data/dataset_diabetes.xlsx', index=False)

#データの可視化
fig, ax = plt.subplots(1, 1, figsize=(10, 5), facecolor='w')
ax.scatter(x_train, y_train)
ax.set(xlabel='s5', ylabel='target', title='Correlation of Diabetes')
plt.grid()
plt.show()

[OUT]
x_train: (8, 1), x_test: (434, 1), y_train: (8,), y_test: (434,)

　　x　　　　	y
0	0.014823	311.0
1	0.007837	122.0
2	0.084495	243.0
3	0.129019	248.0
4	-0.029528	91.0
5	0.079121	281.0
6	-0.018118	142.0
7	0.073410	295.0
    画像

    ４－１．Excelで計算

    　単回帰はExcelでも手動で簡単に計算できます。計算式は前述の通りです。傾き(重み)は2パターンで計算してみました。

    a=1n∑i=1nxiyi−yˉxˉ1n∑i=1nxi2−xˉ2=Cov(X,Y)Var(X)a =\frac{\frac{1}{n}\sum_{i=1}^{n} x_iy_i - \bar{y}\bar{x}}{\frac{1}{n}\sum_{i=1}^{n} x_i^2 - \bar{x}^2}=\frac{Cov(X,Y)}{Var(X)}

    b=yˉ−axˉb=\bar{y}-a\bar{x}

    画像
    画像
    画像

    ４－２．Pythonで計算

    　Pythonで行列・ベクトルなどの計算をする場合はNumpyやPandasが便利であるためこれらを使用して計算します。
    　NumpyとPandasで計算する場合の注意点は下記の通りです。

    1. 分散はデフォルトだとNumpyは標本分散、Pandasは不偏分散で出力

    2. np.cov()はddof=0で標本分散、ddof=1で不偏分散で出力(default:ddof=1)

    [IN]
xarray_train = x_train.values #DataFrameをNumpy配列に変換
print(f'x_train:{xarray_train.shape}, y_train:{y_train.shape}, xarray_train:{type(xarray_train)}, y_train:{type(y_train)}')

#統計量を計算
x_mean = xarray_train.mean() #xの平均値
y_mean = y_train.mean() #yの平均値
x_var = xarray_train.var() #xの標本分散※Numpyのvar()は不偏分散

#共分散を計算
datasets = np.concatenate([xarray_train, y_train.reshape(-1, 1)], axis=1) #列方向に結合してデータセット作成
Cov = np.cov(xarray_train.T, y_train, ddof=0) #ddof=0で標本分散, ddof=1で不偏分散
Cov_xy = Cov[0, 1] #xとyの共分散

slope = Cov_xy / x_var #傾き
intercept = y_mean - slope * x_mean #切片
print(f'xの平均値: {x_mean:.5f}, yの平均値: {y_mean}, xの標本分散: {x_var:.5f}, 共分散: {Cov_xy:.5f}')
print(f'傾き: {slope:.2f}, 切片: {intercept:.3f}')

[OUT]
x_train:(8, 1), y_train:(8,), xarray_train:<class 'numpy.ndarray'>, y_train:<class 'numpy.ndarray'>

xの平均値: 0.04263, yの平均値: 216.625, xの標本分散: 0.00280, 共分散: 2.85528
傾き: 1020.85, 切片: 173.104

    【補足：np.cov()における共分散行列の出力確認】
    　np.cov()で共分散行列を出力できます。出力は下記の通りであり,
    XYの共分散はcovxycov_{xy}やcovyxcov_{yx}から抽出できます。

    Σ=[covxxcovxycovyxcovyy]=[xの分散XYの共分散YXの共分散yの分散]=[sxx2covxycovyxsyy2] \Sigma =\begin{bmatrix} cov_{xx} & cov_{xy} \\ cov_{yx} & cov_{yy} \end{bmatrix} =\begin{bmatrix} xの分散 & XYの共分散 \\ YXの共分散 & yの分散 \end{bmatrix} =\begin{bmatrix} s^{2}_{xx} & cov_{xy} \\ cov_{yx} & s^{2}_{yy} \end{bmatrix} 

    　注意点としてnp.cov()には引数"ddof"があり、この値で標本分散と不偏分散を選択できますがデフォルト値:ddof=1(不偏分散)のため、今回は"ddof=0"を指定しました。
    　出力の形式及び結果は下記の通りです。

    [IN]
xarray_train = x_train.values #DataFrameをNumpy配列に変換
print(f'x_train:{xarray_train.shape}, y_train:{y_train.shape}, xarray_train:{type(xarray_train)}, y_train:{type(y_train)}')

#統計量を計算
x_mean = xarray_train.mean() #xの平均値
y_mean = y_train.mean() #yの平均値
x_var = xarray_train.var() #xの標本分散※Numpyのvar()は不偏分散

#共分散を計算
Cov = np.cov(xarray_train.T, y_train, ddof=0) #ddof=0で標本分散, ddof=1で不偏分散
Cov_xy = Cov[0, 1] #xとyの共分散

slope = Cov_xy / x_var #傾き
intercept = y_mean - slope * x_mean #切片
print(f'xの平均値: {x_mean:.5f}, yの平均値: {y_mean}, xの標本分散: {x_var:.5f}, 共分散: {Cov_xy:.5f}')
print(f'傾き: {slope:.2f}, 切片: {intercept:.3f}')


#数値の確認
print(f'共分散行列(標本分散)\n {Cov}') #ddof=0
print(f'共分散行列(不偏分散)\n {np.cov(xarray_train.T, y_train)}') #ddof=1

#Seabornでヒートマップで共分散行列を可視化
plt.figure(figsize=(6, 6), facecolor='w')
plt.rcParams['font.size'] = 16
#表示：値、x軸ラベル、y軸ラベル、非表示：カラーバー
sns.heatmap(Cov, annot=True, cmap='autumn_r',
            xticklabels=['x', 'y'], 
            yticklabels=['x', 'y'], 
            cbar=False,
            fmt='.4f')

#数値の上に文字を追加
plt.text(0.5, 0.4, '$Var(x):Cov_{xx}$', ha='center', va='center', color='black', fontsize=20)
plt.text(0.5, 1.4, '$Cov_{xy}$', ha='center', va='center', color='black', fontsize=20)
plt.text(1.5, 0.4, '$Cov_{yx}$', ha='center', va='center', color='black', fontsize=20)
plt.text(1.5, 1.4, '$Var(y):Cov_{yy}$', ha='center', va='center', color='black', fontsize=20)

[OUT]
共分散行列(標本分散)
 [[   0.003    2.855]
 [   2.855 6397.234]]

共分散行列(不偏分散)
 [[   0.003    3.263]
 [   3.263 7311.125]]
    画像

    ５．単回帰モデルの実装：ライブラリ編

    　より簡単にモデルを作成するためにライブラリを使用します。

    ５－１．Scikit-learn：LinearRegression

    　Scikit-learnを使用して単回帰を実装します。Scikit-learnのAPIはシンプルであり、ほとんどのモデルは下記で実装できます。

    1. 使用する機械学習モデルをインポート

    2. モデルのインスタンス化

    3. 学習：model.fit()

    4. 学習後のパイパーパラメータ確認

    5. 推論：model.predict()

    [IN]
from sklearn.linear_model import LinearRegression

#単回帰モデルの実装
linear = LinearRegression() #インスタンス化
linear.fit(x_train, y_train) #学習
print(f'傾き: {linear.coef_[0]:.2f}, 切片: {linear.intercept_:.3f}') #傾き、切片の表示
y_pred = linear.predict(x_test) #予測

[OUT]
傾き: 1020.85, 切片: 173.104

    　出力も確認しました。

    [IN]
#データの可視化
fig, ax = plt.subplots(1, 1, figsize=(10, 5), facecolor='w')
ax.scatter(x_train, y_train, label='データ')
ax.plot(x_test, y_pred, color='red', lw=1.0, label='回帰直線')
ax.set(xlabel='s5', ylabel='target', title='Correlation of Diabetes')
ax.text(-0.05, 200, f'y = {linear.coef_[0]:.2f}x + {linear.intercept_:.3f}', ha='center', va='center', color='red', fontsize=14)
plt.grid(), plt.legend()
plt.show()

[OUT]
    画像

    ５－２．Pytorch：nn.Linear()

    　Pytorchで単回帰を実装します。実装フローは下記の通りです。

    1. （線形）モデルをインスタンス化

    2. 各値を標準化(平均0, 分散1)※これが無いと収束しない可能性がある

    3. 各種パラメータ作成(評価関数、最適化関数、学習率など)

    4. 学習モードに切り替える

    5. 推論値を計算(今回はバッチ学習)

    6. 損失関数を使用して"正解値"と"推論値"から誤差を求める

    7. loss.backward()で各パラメータの勾配を求める

    8. optimizer.step()(または手動）によりパラメータを更新する

    9. zero_grad()によりパラメータの勾配を初期化する

    10. 5~9をループ計算する

    [IN]
import torch
import torch.nn as nn
import torch.optim as optim

#データセットの読み込み
diabetes = datasets.load_diabetes()
data, target = diabetes.data, diabetes.target
df_data = pd.DataFrame(data, columns=diabetes.feature_names)
df_data = df_data[['s5']] #単回帰として使用するため1次元のみ抽出

#データの分割：参照用のデータ数を少なくする
x_train, x_test, y_train, y_test = train_test_split(df_data, target, test_size=0.98, random_state=0)

#単回帰モデルの実装
linear = nn.Linear(1, 1) #インスタンス化

#パラメータの設定
lr = 0.05 #学習率
iteration = 1000 #学習回数

criterion = nn.MSELoss() #損失関数の定義
optimizer = optim.SGD(linear.parameters(), lr=lr) #最適化手法の定義

#x_train:dataFrame, y_train:Numpyのため形が異なる
x_train_tensor, y_train_tensor = torch.tensor(x_train.values, dtype=torch.float32), torch.tensor(y_train, dtype=torch.float32) #Numpy配列をTensorに変換
print(f'x_train_tensor:{x_train_tensor.shape}, y_train_tensor:{y_train_tensor.shape}, x_train_tensor:{type(x_train_tensor)}, y_train_tensor:{type(y_train_tensor)}')


#標準化
x_mean = x_train_tensor.mean()
y_mean = y_train_tensor.mean()
x_std = x_train_tensor.std()
y_std = y_train_tensor.std()
print(f'x_mean:{x_mean:.3f}, y_mean:{y_mean:.3f}, x_std:{x_std:.3f}, y_std:{y_std:.3f}')
x_scale = (x_train_tensor- x_mean) / x_std
y_scale = (y_train_tensor - y_mean) / y_std
print(f'x_scale:{x_scale.shape}, y_scale:{y_scale.shape}, xの平均:{x_scale.mean():.1f}, yの平均:{y_scale.mean():.1f}, xの標準偏差:{x_scale.std()}, yの標準偏差:{y_scale.std()}')

#パラメータの表示
for name, param in linear.named_parameters():
    print(f'学習前|パラメータ名: {name}, 値: {param.detach().numpy()}')

#学習
linear.train() #学習モードに設定

print('=' * 50)
for iter in range(iteration):
    optimizer.zero_grad() #勾配の初期化
    output = linear(x_scale) #予測)
    loss = criterion(output, y_scale.unsqueeze(1)) #損失の計算
    loss.backward() #勾配の計算
    optimizer.step() #パラメータの更新

    if (iter + 1) % 100 == 0:
        print(f'iter: {iter + 1}, loss: {loss.item():.4f}')

print('=' * 50)


#パラメータの表示
for name, param in linear.named_parameters():
    print(f'学習後|パラメータ名: {name}, 値: {param.detach().numpy()}')
    

#標準化したモデルのパラメータを元に戻す
params_linaer = {'weight': linear.weight, 'bias': linear.bias}
weight_scale, bias_scale = params_linaer['weight'], params_linaer['bias']
weight_pytorch = weight_scale * y_std / x_std #標準化を戻す
bias_pytorch = bias_scale * y_std + y_mean - weight_pytorch * x_mean #標準化を戻す
weight_pytorch, bias_pytorch = weight_pytorch.detach().numpy(), bias_pytorch.detach().numpy()
print(f'weight_pytorch:{weight_pytorch}, bias_pytorch:{bias_pytorch}')

[OUT]
x_train_tensor:torch.Size([8, 1]), y_train_tensor:torch.Size([8]), x_train_tensor:<class 'torch.Tensor'>, y_train_tensor:<class 'torch.Tensor'>
x_mean:0.043, y_mean:216.625, x_std:0.057, y_std:85.505
x_scale:torch.Size([8, 1]), y_scale:torch.Size([8]), xの平均:-0.0, yの平均:0.0, xの標準偏差:1.0, yの標準偏差:1.0

学習前|パラメータ名: weight, 値: [[0.91]]
学習前|パラメータ名: bias, 値: [-0.466]

==================================================
iter: 100, loss: 0.4763
iter: 200, loss: 0.4763
iter: 300, loss: 0.4763
iter: 400, loss: 0.4763
iter: 500, loss: 0.4763
iter: 600, loss: 0.4763
iter: 700, loss: 0.4763
iter: 800, loss: 0.4763
iter: 900, loss: 0.4763
iter: 1000, loss: 0.4763
==================================================

学習後|パラメータ名: weight, 値: [[0.675]]
学習後|パラメータ名: bias, 値: [0.]

weight_pytorch:[[1020.852]], bias_pytorch:[[173.104]]

    　最後に可視化しました。可視化時は標準化の値を元の値に変換しました。

    [IN]
#推論
linear.eval() #推論モードに設定
_y_pred = linear(x_scale) #予測
y_pred = _y_pred * y_std + y_mean #標準化を戻す

#データの可視化
fig, ax = plt.subplots(1, 1, figsize=(10, 5), facecolor='w')
ax.scatter(x_train, y_train, label='データ')
ax.plot(x_train, y_pred.detach().numpy(), color='red', lw=1.0, label='回帰直線')
ax.set(xlabel='s5', ylabel='target', title='Correlation of Diabetes')
ax.text(0.06, 150, f'y = {float(weight_pytorch):.2f}x + {float(bias_pytorch):.3f}', ha='center', va='center', color='red', fontsize=14)
plt.grid(), plt.legend()
plt.show()

[OUT]
    画像

    ６．補足編：最尤推定による単回帰計算

    　単回帰の最適パラメータ推定、つまり誤差最小化ために最小二乗法を使用しました。最小二乗法とは別に確率論を用いた最尤推定でも計算可能です。
    　一般的には最小二乗法の理解だけで充分のため補足の章として説明します。なお、最尤推定を理解したい方は別記事に記載しております。

    ６－１．サンプルデータ準備

    　データは４章で使用したDiabetesデータセットを使用し、目視で確認できるようにデータ数を1桁に調整しました。

    [IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib
import os

from sklearn import datasets
from sklearn.model_selection import train_test_split

#データセットの読み込み
diabetes = datasets.load_diabetes()
data, target = diabetes.data, diabetes.target
df_data = pd.DataFrame(data, columns=diabetes.feature_names)
df_data = df_data[['s5']] #単回帰として使用するため1次元のみ抽出

#データの分割：参照用のデータ数を少なくする
x_train, x_test, y_train, y_test = train_test_split(df_data, target, test_size=0.98, random_state=0)
print(f'x_train: {x_train.shape}, x_test: {x_test.shape}, y_train: {y_train.shape}, y_test: {y_test.shape}')
print(f'x_train: {type(x_train)}, x_test: {type(x_test)}, y_train: {type(y_train)}, y_test: {type(y_test)}')

#出力用としてデータを結合／出力
array_dataset = np.concatenate([x_train, y_train.reshape(-1, 1)], axis=1)
df_dataset = pd.DataFrame(array_dataset, columns=['x', 'y'])
display(df_dataset)

#データの可視化
fig, ax = plt.subplots(1, 1, figsize=(8, 5), facecolor='w')
ax.scatter(x_train, y_train)
ax.set(xlabel='s5', ylabel='target', title='Correlation of Diabetes')
plt.grid()
plt.show()

[OUT]
x_train: (8, 1), x_test: (434, 1), y_train: (8,), y_test: (434,)
x_train: <class 'pandas.core.frame.DataFrame'>, x_test: <class 'pandas.core.frame.DataFrame'>, y_train: <class 'numpy.ndarray'>, y_test: <class 'numpy.ndarray'>


	x	　　　　　y
0	0.014823	311.0
1	0.007837	122.0
2	0.084495	243.0
3	0.129019	248.0
4	-0.029528	91.0
5	0.079121	281.0
6	-0.018118	142.0
7	0.073410	295.0
    画像

    ６－２．正規分布と尤度

    　確率モデルを扱う上で理解が必須である正規分布と尤度について簡単に説明します。詳細は別記事参照のこと。

    　６－２－１．正規分布

    　正規分布とは連続値を扱う確率密度分布の一つであり自然現象などにおけるばらつきを表現できます。平均値μ\mu, 分散σ2\sigma^2を用いて下記数式で表現でき、μ=0\mu=0, 分散σ2=1\sigma^2=1の場合を標準正規分布と呼びます。

    N(μ,σ2)=f(x)=12πσ2exp⁡(−(x−μ)22σ2)\mathcal{N}(\mu, \sigma^2)=f(x)=\dfrac{1}{\sqrt{2\pi\sigma^2}}\exp(-\dfrac{(x-\mu)^ 2}{2\sigma^ 2})

    [IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib
import os
from scipy import stats

def plot_line(x, y, ax, label, x_label, y_label, title):
    ax.plot(x, y, label=label)
    ax.set(xlabel=x_label, ylabel=y_label, title=title)
    ax.grid(), ax.legend()
    plt.show()

def normal_distribution(x, mu, sigma):
    return 1 / (np.sqrt(2 * np.pi) * sigma) * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) # 正規分布

x = np.linspace(-5, 5, 1000)
y = normal_distribution(x, 0, 1) # 平均0, 標準偏差1の正規分布(標準正規分布)

fig, ax = plt.subplots(figsize=(10, 6), facecolor="w")
plot_line(x, y, ax, "標準正規分布", "x", "y", "標準正規分布")

[OUT]
    画像

    　６－２－２．尤度

    　尤度とはデータに対するパラメータの尤もらしさを表します。正規分布はパラメータとして平均値μ\mu, 分散σ2\sigma^2を持っていますがデータxxに対して尤もらしいパラメータを表現します。

    N(μ,σ2)=f(x)=12πσ2exp⁡(−(x−μ)22σ2)(−∞<x<∞)\mathcal{N}(\mu, \sigma^2)=f(x)=\dfrac{1}{\sqrt{2\pi \sigma^2}} \exp \left(-\frac{(x - \mu)^2} {2\sigma^2} \right) \hspace{20px} (-\infty < x < \infty)

    　例としてデータ[10, 14, 20, 22]があり正規分布のパラメータμ\mu, σ2\sigma^2を変化させたとき尤もらしいパラメータが存在します。
    　下図の例では右下の図、つまり正規分布の中心にデータが集まっており確率的にデータの出現数が多い場所が尤もらしいと判断できます。

    画像
    [IN]
def plot_line(x, y, ax, label, x_label, y_label, title, varbose=True):
    if varbose:
        ax.plot(x, y, label=label)
    else:
        ax.scatter(x, y, label=label, color="r", edgecolors="k")
        ax.set(xlabel=x_label, ylabel=y_label, title=title)
        ax.grid(), ax.legend()

def calc_likelihood(data, mu, sigma):
    likelihoods = stats.norm.pdf(data, mu, sigma)
    likelihood = np.prod(likelihoods)
    return likelihood
    
    
data = np.array([10, 14, 20, 22])
mu, sigma = data.mean(), data.std()
print(f'平均: {mu}, 標準偏差: {sigma}')

x = np.linspace(-30, 30 , 100)
y1 = normal_distribution(x, 0, 1) # 平均0, 標準偏差1の正規分布(標準正規分布)
y2 = normal_distribution(x, 5, 5) # 平均5, 標準偏差5の正規分布
y3 = normal_distribution(x, 16.5, 0.5) # 平均16.5, 標準偏差0.5の正規分布
y4 = normal_distribution(x, mu, sigma) # 平均mu, 標準偏差sigmaの正規分布

xs = [x, x, x, x]
ys = [y1, y2, y3, y4]
labels = ["標準正規分布", r'$\mu=5, \sigma=5$', r'$\mu=16.5, \sigma=0.5$', r'$\mu=calc, \sigma=calc$']

fig, axs = plt.subplots(2, 2, figsize=(12, 10), facecolor="w")
axs = axs.ravel() #データ1 
plot_line(xs[0], ys[0], axs[0], labels[0], "x", "y", labels[0], varbose=True)
plot_line(data, np.zeros_like(data), axs[0], 'data', "x", "y", labels[0], varbose=False) #データ2 
plot_line(xs[1], ys[1], axs[1], labels[1], "x", "y", labels[1], varbose=True)
plot_line(data, np.zeros_like(data), axs[1], 'data', "x", "y", labels[1], varbose=False) #データ3 
plot_line(xs[2], ys[2], axs[2], labels[2], "x", "y", labels[2], varbose=True)
plot_line(data, np.zeros_like(data), axs[2], 'data', "x", "y", labels[2], varbose=False) #データ4 
plot_line(xs[3], ys[3], axs[3], labels[3], "x", "y", labels[3], varbose=True)
plot_line(data, np.zeros_like(data), axs[3], 'data', "x", "y", labels[3], varbose=False)

print(rf'$\mu=0, \sigma=1$の尤度 : {calc_likelihood(data, 0, 1):.10f}')
print(rf'$\mu=5, \sigma=5$の尤度 : {calc_likelihood(data, 5, 5):.10f}')
print(rf'$\mu=16.5, \sigma=0.5$の尤度 : {calc_likelihood(data, 16.5, 0.5):.10f}')
print(rf'$\mu=calc, \sigma=calc$の尤度 : {calc_likelihood(data, mu, sigma):.10f}')

[OUT]
平均: 16.5, 標準偏差: 4.769696007084728
$\mu=0, \sigma=1$の尤度 : 0.0000000000
$\mu=5, \sigma=5$の尤度 : 0.0000000002
$\mu=16.5, \sigma=0.5$の尤度 : 0.0000000000
$\mu=calc, \sigma=calc$の尤度 : 0.0000066235

（図は上に記載）

    　この尤もらしさを関数で表現したのが尤度（尤度関数）LLです。複数データの場合は各尤度（正規分布）の総乗となります。なお正規分布の縦軸は確率密度であり正規分布の中心ほどデータの出やすさが大きくなるため、尤度が最大（全データが同時に出現確率が最大）のパラメータが最適なパラメータとなります。

    L(μ,σ)=∏i=1nf(x)=∏i=1n12πσ2exp⁡(−(x−μ)22σ2)L(\mu, \sigma) = \prod_{i=1}^n f(x) = \prod_{i=1}^n \dfrac{1}{\sqrt{2\pi\sigma^2}}\exp(-\dfrac{(x-\mu)^ 2}{2\sigma^ 2})

    ６－３．単回帰のモデル設計：誤差項と尤度関数

    　単回帰のモデルを傾きaa, 切片bbとして考えますが、これにϵ\epsilonを追加します。ϵ\epsilonは誤差項（モデルのばらつき）です。

    yi=axi+b+ϵiy_i = ax_i + b + \epsilon_i

    　ここで誤差項ϵ\epsilonは平均0, 分散σ2\sigma^ 2の正規分布に従うと仮定します。この誤差はyi=axi+b+ϵiy_i = ax_i + b + \epsilon_iでのばらつきであるため、推論値を中心とする正規分布となります（追って可視化します）。

    ϵi∼N(0,σ2)=f(ϵi)=12πσ2exp⁡(−ϵi22σ2)\epsilon_i \sim N(0, \sigma^2)=f(\epsilon_i)=\dfrac{1}{\sqrt{2\pi\sigma^2}}\exp(-\dfrac{\epsilon_{i}^ 2}{2\sigma^ 2})

    　単回帰モデルをϵi=yi−(axi+b)\epsilon_i=y_i -(ax_i + b)に変形して正規分布f(ϵi)f(\epsilon_i)に代入すると各データにおける確率密度関数（正規分布）f(yi)f(y_i)が得られます。

    f(yi)=12πσ2exp⁡(−(yi−(axi+b))22σ2)f(y_i) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(y_i - (ax_i + b))^2}{2\sigma^2} \right)

    　確率密度関数f(yi)f(y_i)を尤度とすると、パラメータaa, bbを用いて尤度関数L(a,b)L(a,b)は下記の通り表せます。

    L(a,b)=∏i=1nf(yi)L(a, b) = \prod_{i=1}^n f(y_i)

    　最後に計算しやすいように対数を取り対数尤度を取得します。

    log⁡L(a,b)=−n2log⁡(2πσ2)−12σ2∑i=1n(yi−(axi+b))2\log L(a, b) = -\frac{n}{2} \log(2 \pi \sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n} (y_i - (ax_i + b))^2

    ６－４．誤差項εの確率

    　最尤推定は「尤度関数を最大化するパラメータを推定」するため、今回の例では「誤差項εの確率が最大となるパラメータを推定」することになります。ここで理解しにくい”誤差項εの確率が最大”について説明します。
    　まずは結論となるイメージ図を紹介します。

    画像
    画像
    画像

    　単回帰モデル：yi=axi+b+ϵiy_i = ax_i + b + \epsilon_iの誤差項ϵi \epsilon_iはデータのばらつきやノイズを表現する役割を果たしており、そのばらつきの確率モデルを正規分布として表現しています。
    　推論値yiy_iのばらつき（分布）内にデータtit_iが存在すれば”分布内にデータが存在する確率（密度分布）が高い->尤度が大きい”ことになります。各データにおける確率の積を取り、それが最大になるということは「推論値yiy_iで想定したばらつき・ノイズϵi \epsilon_iにデータtit_iが（同時に）存在する確率が最も高い」ことになります。
    　つまりデータに対して尤もらしいことを表現しており、尤度関数が最大になるパラメータを見つける（最尤推定）ことで最適なパラメータを決定できます。

    　推論値のばらつきの分布f(ϵi)f(\epsilon_i)が単回帰モデルやデータに対してどのように存在しているのかを可視化するコードを記載しました。なお単回帰モデルのパラメータは最適値を使用しており、かつ上図の分布に塗りつぶしを加えて更に分布が分かりやすいように表現しました。

    [IN]
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from typing import Tuple
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

def plot_data_and_normal_distribution_3D(x_data: np.ndarray, y_data: np.ndarray, slope: float, intercept: float, sigma: float) -> None:
    # 3Dプロットの準備
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # 観測データを散布図でプロット
    for i, (x, y) in enumerate(zip(x_data, y_data)):
        ax.scatter(x, y, 0, color='red', edgecolor='black', linewidth=1, zorder=10)
        ax.text(x, y, 0, f'({x:.3f}, {y:.1f})', fontsize=5)

    # 線形回帰モデルのパラメータを用いて、回帰直線を表示
    y_preds = slope*x_data + intercept #y =ax+b
    ax.plot(x_data, y_preds, np.zeros_like(x_data), color='red', linewidth=0.5, ls = '--')

    #各データポイントにおける正規分布 (誤差項ε)を表示
    for i in range(len(x_data)):
        y_pred = slope * x_data[i] + intercept
        y_normal = np.linspace(y_pred-4*sigma, y_pred+4*sigma, 100) #正規分布の横軸 (y軸)
        z_normal = norm.pdf(y_normal, loc=y_pred, scale=sigma) #正規分布の確率密度関数 
        ax.plot(np.full_like(y_normal, x_data[i]), y_normal, z_normal, color='blue') #正規分布をyz平面で青でプロット 
        
        #正規分布の頂点表示 ：単回帰の線上に薄い青色の線を引く
        ax.plot([x_data[i], x_data[i]], [y_preds[i], slope * x_data[i] + intercept], [0, np.max(z_normal)], color='blue', alpha=0.5, lw=0.3, ls='--')

        # 正規分布（誤差項）をyz平面で薄い青で塗りつぶす
        verts = [list(zip(np.full_like(y_normal, x_data[i]), y_normal, z_normal))]
        ax.add_collection3d(Poly3DCollection(verts, facecolors='blue', edgecolors='blue', alpha=0.1))

    ax.set_xlabel('s5', fontsize=10)
    ax.set_ylabel('target(diabetes)', fontsize=10)
    ax.set_zlabel(r'誤差項\epsilonの確率密度関数', fontsize=10)
    ax.tick_params(labelsize=6)

    plt.show()

# データ
x_data = x_train.values.ravel() #dataframe ->np.arrayから1次元配列に変換
y_data = y_train #コピー 

# 最小二乗法を用いて線形回帰モデルのパラメータを計算
slope, intercept = np.polyfit(x_data, y_data, 1)
print(f'単回帰の傾き: {slope:.3f}, 切片: {intercept:.3f}')

# 散布図と各データポイントにおける正規分布を表示
sigma = 20 #標準偏差を仮定 (適当に設定)
plot_data_and_normal_distribution_3D(x_data, y_data, slope, intercept, sigma)


[OUT]
単回帰の傾き: 1020.851, 切片: 173.104

    　結果を見ると左から４つ目のデータは推論値のばらつき（誤差項）に対して大きくずれておりますがそれ以外は軒並み分布の中央に近いことが確認できます。よって各データの確率密度（尤度）が高く、尤度関数（尤度の総乗）が最大化され、パラメータとして尤もらしいことが確認できます。

    画像
    画像
    画像

    ６－５．最尤推定による解析解の計算

    　求めた対数尤度関数を微分することで最適解を求めます。結論として最小二乗法と最尤推定で得られる最適パラメータは同じになります。

    log⁡L(a,b)=−n2log⁡(2πσ2)−12σ2∑i=1n(yi−(axi+b))2\log L(a, b) = -\frac{n}{2} \log(2 \pi \sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n} (y_i - (ax_i + b))^2

    　対数尤度関数を微分します。

     ∂log⁡L(a,b)∂a=1σ2∑i=1nxi(yi−(axi+b))=0∂log⁡L(a,b)∂b=1σ2∑i=1n(yi−(axi+b))=0\begin{aligned}  \frac{\partial \log L(a, b)}{\partial a} = \frac{1}{\sigma^2} \sum_{i=1}^n x_i (y_i - (ax_i + b))=0 \\ \frac{\partial \log L(a, b)}{\partial b} = \frac{1}{\sigma^2} \sum_{i=1}^n (y_i - (ax_i + b))=0 \end{aligned}

    　偏微分 ∂log⁡L(a,b)∂a \frac{\partial \log L(a, b)}{\partial a} , ∂log⁡L(a,b)∂b\frac{\partial \log L(a, b)}{\partial b}の連立方程式を解きます。

    {∑i=1nxi(yi−(axi+b))=∑i=1nxiyi−a∑i=1nxi2−b∑i=1nxi=0∑i=1n(yi−(axi+b))=∑i=1nyi−a∑i=1nxi−nb=0\begin{cases} \sum_{i=1}^n x_i (y_i - (ax_i + b)) =\sum_{i=1}^nx_iy_i -a\sum_{i=1}^nx_i^2 -b\sum_{i=1}^nx_i= 0 \\ \sum_{i=1}^n (y_i - (ax_i + b)) =\sum_{i=1}^ny_i-a\sum_{i=1}^nx_i-nb =0 \end{cases}

    よって最適解は下記の通りです。

    a=n∑i=1nxiyi−(∑i=1nxi)(∑i=1nyi)n∑i=1nxi2−(∑i=1nxi)2a = \frac{n\sum_{i=1}^n x_i y_i - (\sum_{i=1}^n x_i)(\sum_{i=1}^n y_i)}{n\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2}

    b=(∑i=1nyi)∑i=1nxi2−(∑i=1nxiyi)(∑i=1nxi)n∑i=1nxi2−(∑i=1nxi)2b = \frac{(\sum_{i=1}^n y_i)\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i y_i)(\sum_{i=1}^n x_i)}{n\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2}

    参考資料

    今後の教材用

    あとがき

    　目標は３月中に全部作成したい。

    　

     
     

    KIYO

     
     
    普段は製造業で企画/開発/設計しております。記事はプログラミング・機械学習、IoT関係の記事をメインで作成し、なるべく１つの記事で知りたいことを網羅していきます。内容は学術的より実装・アウトプット（ほしくなるもの）を重視して作成しています。 面白そうな仕事があればやりたいです！

    あなたへのおすすめ