5. 線形回帰 （Linear Regression）#

Hide code cell source 
import numpy as np 
from numpy.random import default_rng
import plotly.express as px
import pandas as pd 
import plotly.graph_objs as go

SEED = 2022_11_10

5.1. データの準備#

まずはデモデータを作成しましょう。デモデータの作成は以下の記事を参考にしました。

@sakabe(株式会社ヘッドウォータース). “【機械学習】練習用データセットの作成【身長・体重】”. Qiita. 2019-05-23. https://qiita.com/sakabe/items/df7fac64e6ffbb50848c, (参照 2022-11-10).

# デモデータの生成
rng = default_rng(SEED)

m_height = rng.normal(171.7, 6.6, 1000)
f_height = rng.normal(158.3, 5.7, 1000)
m_bmi = rng.normal(23.12, 2, 1000)
f_bmi = rng.normal(20.82, 2, 1000)

m_weight = (m_height * 0.01) ** 2 * m_bmi
f_weight = (f_height * 0.01) ** 2 * f_bmi

male = pd.DataFrame()
female = pd.DataFrame()
male["height"] = m_height
male["weight"] = m_weight
male["sex"] = "male"

female["height"] = f_height
female["weight"] = f_weight
female["sex"] = "female"

df = pd.concat([male,female])
df
height weight sex
0 180.097661 79.823196 male
1 168.651353 62.357759 male
2 172.127079 68.235781 male
3 184.574798 82.567780 male
4 162.787830 55.248077 male
... ... ... ...
995 145.860850 49.764990 female
996 159.824910 61.096695 female
997 161.421906 44.226520 female
998 158.992977 48.269525 female
999 152.012887 55.610335 female

2000 rows × 3 columns

fig = px.scatter(
    df, 
    x="weight", y="height", 
    height=800, width=800,
    title="身長体重のdemo data",
    color="sex",
    #trendline="ols"
    )

fig

身長体重データを散布図にしました。この散布図に対してよくフィットするような直線を引くことができれば、データの大まかな傾向を見ることができます。

5.2. アルゴリズム#

直線は一次関数です。\(y=ax+b\)で表される直線を描きたいので、\(a\)と\(b\)の二つの定数が分かれば良いことになります。 ここで\(a\)は傾き、\(b\)はy切片です。

一般的に用いられる最小二乗法を用いて、これらの値を求めましょう。

\[a = \frac{\text{XとYの共分散}}{\text{Xの分散}} = \frac{S_{x,y}}{S_{x,x}} = {\frac{\displaystyle \frac{1}{n}\sum_{i=1}^{n}(X_i-\hat{X})(Y_i-\hat{Y})}{\displaystyle \frac{1}{n}\sum_{i=1}^{n}(X_i-\hat{X})^2}}\]

この式をPythonを使って算出します。NumPyのnp.cov関数を用いることで、分散共分散行列が求まります。\(X\)と\(Y\)について求める場合、この関数は次のような行列を返します。

\[\begin{split} \sum = \begin{pmatrix} s_{xx}&s_{xy}\\ s_{yx}&s_{yy}\\ \end{pmatrix} \end{split}\]

ここで\(\sum\)は分散共分散行列を表します。\(S_{xy}\)は\(x\)と\(y\)の共分散で、\(S_{xx}\)は\(x\)の分散です。また、\(S_{xy}\)と\(S_{yx}\)は同じ値になります。分散共分散行列は以下のように定義されます（求められます）。

\[\begin{split} \sum = E\left [\left\{\begin{pmatrix}x\\y\end{pmatrix} - E \begin{pmatrix}x\\y\end{pmatrix} \right\} \left\{\begin{pmatrix}x\\y\end{pmatrix} - E \begin{pmatrix}x\\y\end{pmatrix} \right\} ^\top\right ] \end{split}\]

後は\(a\)を求める式に当てはめるだけです。

5.3. 実験#

5.3.1. 実験#

cov_mat = np.cov(df["weight"],df["height"], ddof=0)
print("分散共分散行列 Σ:\n", cov_mat)

s_xx = cov_mat[0,0]
s_xy = cov_mat[0,1]
s_yy = cov_mat[1,1]
a = s_xy/s_xx
print(f"{a=}")

分散共分散行列 Σ:
 [[115.54049919  81.77930932]
 [ 81.77930932  82.26103428]]
a=0.7077977842593172

次に、y切片を求めます。

\(a\)が分かったので、\(x\)と\(y\)が分かれば\(y=ax+b\)を式変形することで\(b\)が求められます。
\(x\)と\(y\)にはそれぞれの平均値を利用しましょう。

# df["height"].mean() = a * df["weight"].mean() + b 
b = df["height"].mean() -( a * df["weight"].mean() )
b

122.09958131811504

よって一次関数は以下の式です。

\[y=0.7077977842593172 \times x + 121.62779834449698\]

これをグラフにしてみましょう。

fig = px.scatter(
    df, 
    x="weight", y="height", 
    height=800, width=800,
    title="身長体重のdemo data, データ全体を元にした回帰直線",
    color="sex",
    #trendline="ols"
    )
linear_x = np.linspace(df.weight.min(), df.weight.max())
predicted_y = a*linear_x +b
fig.add_scatter(x=linear_x,y=predicted_y, mode="lines", name="Trace Line")

fig.show()

今回はmaleとfemaleの二つのクラスがあることを考慮せずに回帰直線を引きました。もしも男性と女性それぞれに対して回帰直線を引きたい場合は、それぞれのクラスに属するデータを抽出してから上述の流れで\(a\)と\(b\)を求めましょう。

fig2 = px.scatter(
    df, 
    x="weight", y="height", 
    height=800, width=800,
    title="身長体重のdemo data, クラス毎の回帰直線",
    color="sex",
    trendline="ols",
    trendline_color_override= "green"
    )

fig2.show()

5.3.2. 相関係数#

\[ r =\frac{[x \text { と } y \text { の共分散 }]}{[x \text { の標準偏差 }][y \text { の標準偏差 }]} =\frac{S_{xy}}{S_{xx}S_{yy}} =\frac{\displaystyle \frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\sqrt{\displaystyle \frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)^2} \sqrt{\displaystyle \frac{1}{n} \sum_{i=1}^n\left(y_i-\bar{y}\right)^2}} =\frac{\displaystyle \sum_{i=1}^n\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\sqrt{\displaystyle \sum_{i=1}^n\left(x_i-\bar{x}\right)^2} \sqrt{\displaystyle \sum_{i=1}^n\left(y_i-\bar{y}\right)^2}} \]

\(x\)と\(y\)の相関係数である\(r\)は、この二つの変数の直線的な関係の強さを示す指標です。\(r\)は標本の相関係数を指す記号であり、母集団（全体）の相関係数は\(\rho\)で表します。
相関係数は1~-1の値をとります。おおよその目安は以下の通りです。

\[\begin{split} \begin{array}{|c|c|} \hline \text { 相関係数 } r & \begin{array}{c} \text { 相関の強さ } \\ (※ p \text { 值 }<0.05) \end{array} \\ \hline 0.7 \leq r \leq 1.0 & \text { 強い正の相関 } \\ \hline 0.4 \leq r \leq 0.7 & \text { 正の相関 } \\ \hline 0.2 \leq r \leq 0.4 & \text { 弱い正の相関 } \\ \hline-0.2 \leq r \leq 0.2 & \text { (ほとんど相関がない } \\ \hline-0.4 \leq r \leq-0.2 & \text { 弱い負の相関 } \\ \hline-0.7 \leq r \leq-0.4 & \text { 負の相関 } \\ \hline-1.0 \leq r \leq-0.7 & \text { 強い負の相関 } \\ \hline \end{array} \end{split}\]

この例の相関係数\(r\)を求めてみましょう。

r = s_xy / (s_xx**0.5 * s_yy**0.5)
print(f"相関係数{r=}")

相関係数r=0.8388402011676876

強い正の相関があることが分かります。

https://atarimae.biz/archives/7966

https://agirobots.com/var-cov-matrix/

https://www.albert2005.co.jp/knowledge/statistics_analysis/multivariate_analysis/single_regression

https://qiita.com/sakabe/items/df7fac64e6ffbb50848c

https://qiita.com/innovation1005/items/b712ce54a7a697a9bf03

最後に、相関係数Rとデータの関係を確認しておきましょう。

def build_r_animation_frame():
    sample_df = []
    for i in range(-10,11):
        mean = np.array([0, 0])
        cov = np.array([
            [1, 0.1*i],
            [0.1*i, 1]]) 
        r = cov[0,1] / (cov[0,0]**0.5 * cov[1,1]**0.5)
        x, y = rng.multivariate_normal(mean, cov, 2000).T
        sample_df.append(pd.DataFrame(dict(x=x,y=y,r=r,)))
    return pd.concat(sample_df)

r_animation_df = build_r_animation_frame()
r_animation = px.scatter(
    r_animation_df, 
    x="x",y="y",
    trendline="ols", trendline_color_override="red", 
    width=800,height=800, 
    animation_frame="r",
    title="相関係数rとデータの関係（r=-1:負の相関, r=0:無相関, r=1:正の相関）",
    range_x=[-3,3], range_y=[-3,3],
    )
#r_animation.layout.updatemenus[0]["buttons"][0]["args"][1]["transition"]["duration"] = 0
r_animation.show()

5.4. Scikit-learnによる実行#

linear regressionを実装したパッケージはたくさんありますが、今後のノートブックと統一するためにscikit-learn (以降sklearn)を利用します。

sklearnは以下のステップで実行可能です。

  1. model=sklearn.Model(hyperparams)によるインスタンス生成

  2. model.fit(X_train,y_train)による訓練

  3. model.predict(X_test) による予測 or model.transform(X_test)による入力行列の変形

  4. Optional: model.score(X_test, y_test)による評価値の算出

教師データの場合、これ以前にデータをX_train,X_test, y_train,y_testに分ける必要があります。通常のsklearnではValidation setを用意することなく、Training setのみを使って訓練を行います。必要であれば、Training set (X_train, y_train)を更にX_train, X_val, y_train, y_valに分けることでValidation setを作成します。これにはtrain_test_split関数を使うと簡単です。

例題として、iris datasetのsepal length (cm)とsepal width (cm)を利用し、回帰直線と回帰係数を求めます。

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

load_irisはas_frame=TrueでDataFrameを返します。target(label)の有無などの複数のdfを持ったdictを返すので、”frame”(※全て入ったdf)を指定します。

dataset_df = load_iris(as_frame=True)["frame"] 
print(dataset_df.shape)
display(dataset_df.head())

(150, 5)
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm) target
0 5.1 3.5 1.4 0.2 0
1 4.9 3.0 1.4 0.2 0
2 4.7 3.2 1.3 0.2 0
3 4.6 3.1 1.5 0.2 0
4 5.0 3.6 1.4 0.2 0 
corr = dataset_df.corr()
display(corr)
px.imshow(corr, text_auto=True, title="相関係数")
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm) target
sepal length (cm) 1.000000 -0.117570 0.871754 0.817941 0.782561
sepal width (cm) -0.117570 1.000000 -0.428440 -0.366126 -0.426658
petal length (cm) 0.871754 -0.428440 1.000000 0.962865 0.949035
petal width (cm) 0.817941 -0.366126 0.962865 1.000000 0.956547
target 0.782561 -0.426658 0.949035 0.956547 1.000000 
dataset_df.sort_values('sepal length (cm)', inplace=True)

X_train,X_test,y_train,y_test = train_test_split(dataset_df[['sepal length (cm)']],dataset_df['petal length (cm)'], shuffle=False)
for arr in [X_train,X_test,y_train,y_test]:
    print(f"{arr.shape=}, {type(arr)}")

arr.shape=(112, 1), <class 'pandas.core.frame.DataFrame'>
arr.shape=(38, 1), <class 'pandas.core.frame.DataFrame'>
arr.shape=(112,), <class 'pandas.core.series.Series'>
arr.shape=(38,), <class 'pandas.core.series.Series'>

model = LinearRegression()
model.fit(X_train, y_train)

LinearRegression()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
X_train_range = np.linspace(X_train.min(), X_train.max())
predicted_y_train_range = model.predict(X_train_range)

fig= px.line(
    x=X_train_range.squeeze(),
    y=predicted_y_train_range, 
    title=f"regression line (回帰係数={model.coef_[0]}, 切片={model.intercept_})",
    )
fig.add_trace(
    go.Scatter(
        x=X_train.to_numpy().squeeze(), 
        y=y_train, 
        mode="markers", 
        name="original data point",
        ),
    )
fig.show()

/Users/mriki/.pyenv/versions/miniforge3-4.10.3-10/envs/datasci/lib/python3.10/site-packages/sklearn/base.py:450: UserWarning:

X does not have valid feature names, but LinearRegression was fitted with feature names

example_df = X_train.copy()
example_df["petal length (cm)"] = y_train.copy()
example_df["predicted petal length (cm)"] = model.predict(X_train)
example_df["kind"] = "training"
_example_df = X_test.copy()
_example_df["petal length (cm)"] = y_test.copy()
_example_df["predicted petal length (cm)"] = model.predict(X_test)
_example_df["kind"] = "test"
example_df = pd.concat([example_df, _example_df])
example_df
sepal length (cm) petal length (cm) predicted petal length (cm) kind
13 4.3 1.1 0.383410 training
42 4.4 1.3 0.619521 training
38 4.4 1.3 0.619521 training
8 4.4 1.4 0.619521 training
41 4.5 1.3 0.855632 training
... ... ... ... ...
122 7.7 6.7 8.411181 test
118 7.7 6.9 8.411181 test
117 7.7 6.7 8.411181 test
135 7.7 6.1 8.411181 test
131 7.9 6.4 8.883403 test

150 rows × 4 columns

fig= px.line(
    example_df,
    x="sepal length (cm)",
    y= "predicted petal length (cm)",
    color="kind" 
    #title=f"regression line (回帰係数={model.coef_[0]}, 切片={model.intercept_})",
    )
fig.add_trace(
    go.Scatter(
        x=example_df["sepal length (cm)"], 
        y=example_df["petal length (cm)"],
        mode="markers", 
        name="original data point",
        ),
    )
fig.show()

5.4.1. Palmer Penguin Datasetを使った演習#

  1. bill_length_mmを説明変数としてbill_depth_mmを予測する回帰直線を書け。

  2. \(y=ax+b\)の形で式を表示せよ。

  3. 訓練データ、テストデータに対する決定係数を表示せよ。

rianrajagede/penguin-python