線形回帰を理解してPythonで1からコードを書いてみよう

線形回帰について現時点で私が理解していることと、その過程で書いたPythonのコードをまとめます。数学はガバガバですがご了承ください。

## 回帰分析

回帰分析とは、原因となったデータと結果となったデータを分析し、その二つのデータの関係性を導き出すための手法です。線形と非線形に分けることができます。

### 線形回帰

身長と体重の二つの変数からなるデータを考えてみましょう。以下のデータがあったとします。

| 身長 | 体重 |
| --- | --- |
| 170 | 63 |
| 175 | 67 |
| 164 | 52 |
| 157 | 50 |
| 180 | 75 |

これを散布図で表してみます

![scatter.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/8b39ec2d-32f0-59cf-57b1-828a93a19c6f.png)

この図に対して、**いい感じな直線**を引いてみましょう。

![simple_linear.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/9ae8ebb6-949b-9dab-7ab4-cbd72a6f9942.png)

できました。
この**いい感じな直線**を計算で求めるのが線形回帰となります。

<br>

### 非線形回帰

本稿で扱うのは線形回帰ですが、非線形も軽く触れておきます。
先ほどはいい感じな直線が引けそうなデータでしたが、その様なデータばかりではありません。例えば下の図。

![simple_nonlinear.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/acaeb0e8-5cea-60da-12e4-747f7604d9a9.png)

この様なデータの場合、直線ではなく**曲線**を引くことを考えます。

![sin.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/330a4e07-6bd8-d10c-d7db-b41239b9f956.png)

こんな感じです。~~非線形回帰では、**いい感じの曲線**を計算で求めます。~~

（画像は多分$\sin x$なので非線形と言えますが、例えば適当な$m$次の項までの線型結合によって求められていた場合などはこれは線形回帰になります。コメントをいただいていたのに長らく放置していてすいません。）

<br>

### 単回帰と重回帰

上で挙げた二つの例はいずれも**単回帰**です。単回帰とは、原因となるデータが一次元の場合を指します。上の線形回帰の例では、「身長」という一種類のデータだけで体重を予測していたので、単回帰となります。
一方で、原因となるデータが多次元の場合は**重回帰**となります。例を挙げると、身長の他に「睡眠時間」「摂取カロリー」「運動量」などが加わった多次元のデータから体重を予測するなどです。

<br>

それでは実際に線形回帰をやっていきましょう

## 単回帰

まずは単回帰からです。最小二乗法を使っていい感じの直線を求めていきます。
最小二乗法では、線を引いた時、その線とデータの差の二乗和が最小になる様に線を引きます。図で見てみましょう。

![jijou.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/be21ea13-4bae-3c9a-41af-9ec5ee15fa2c.png)

この緑の線の長さの和が最小になる直線を求めようと言うことですね。

<br>

では計算していきます。

求めたい直線を

$f(k) = ak + b$

と置きます。 また、原因となったデータ$x$と、結果となったデータ$y$を

$x = (x_1 \ x_2\ \cdots \ x_n)$
$y = (y_1 \ y_2\ \cdots \ y_n)$

とします。因みに、原因となったデータを**説明変数**、結果となったデータを**目的変数**と呼びます。以下、本稿でもこの表現を用います。

次に、最小化したい関数$E$を以下の様に定義します。これが直線と点の距離の二乗和です。

```math
E = \sum_{i=1}^{n} \big(f(x_i) - y_i\big)^2 = \sum_{i=1}^{n} \big(ax_i + b - y_i\big)^2
```

そしてこの$E$を最小にする$a$と$b$を求めます。$E$は$a$と$b$の二次関数で、二つの係数はいずれも正なので、偏微分をして傾きが0になる値を探せばいいですね。

```math
\begin{align}
\frac{\partial E}{\partial a} &= 2\sum_{i=1}^{n} x_i \big(ax_i + b - y_i\big)=0 \, \cdots ① \, \\
\frac{\partial E}{\partial b} &= 2\sum_{i=1}^{n} \big(ax_i + b - y_i\big)=0 \, \cdots \, ②
\end{align}
```

この連立方程式を解いていきます。まず両辺を2で割り、展開してから$y_i$の項を移行します。

```math
\begin{align}
a\sum_{i=1}^{n}x_i^2 + \sum_{i=1}^{n}bx_i &= \sum_{i=1}^{n}x_iy_i \ \cdots ① \\
a\sum_{i=1}^{n}x_i + bn &= \sum_{i=1}^{n}y_i \ \cdots ②
\end{align}
```

②より、$b$が求まります。

```math
b = -a\frac{\sum_{i=1}^{n}x_i}{n} + \frac{\sum_{i=1}^{n}y_i}{n} = -a\bar{x} + \bar{y}
```

$b$を$x$と$y$の平均で表すことができました。これを①に代入しして色々いじります。

```math
\begin{align}
a\sum_{i=1}^{n}x_i^2 + \sum_{i=1}^{n}(-a\bar{x} + \bar{y})x_i &= \sum_{i=1}^{n}x_iy_i \\
a\sum_{i=1}^{n}x_i^2 - a\bar{x}\sum_{i=1}^{n}x_i &= \sum_{i=1}^{n}x_iy_i -  \bar{y}\sum_{i=1}^{n}x_i\\
\end{align}
```

両辺をnで割り、$a$を求めます。

```math
\begin{align}
a(\bar{x^2}-\bar{x}^2) &= \bar{xy} - \bar{x}\bar{y} \\
a &= \frac{\bar{xy} - \bar{x}\bar{y}}{\bar{x^2}-\bar{x}^2} = \frac{\sigma_{xy}}{\sigma_x^2}
\end{align}
```

$a$を$xy$の共分散と$x$の分散で表すことができました。まとめるとこうなります。

```math
\begin{align}
a &= \frac{\sigma_{xy}}{\sigma_x^2} \\
b &= -a\bar{x} + \bar{y}
\end{align}
```

これを $f(x) = ax + b$ に代入すれば完成ですね。

### Pythonで実装

以上の内容をPythonで実装していきましょう

```python
import numpy as np
import matplotlib.pyplot as plt

class simple_linear:
    # 回帰係数と切片の計算
    def fit(self, X, y):
        self.X = X  # 説明変数
        self.y = y  # 目的変数
        self.a = np.cov(X, y)[0, 1] / np.var(X)  # 回帰係数
        self.b = -self.a * X.mean() + y.mean()  # 切片
    
    # 予測
    def predict(self, X):
        return self.a*X + self.b
    
    # プロット
    def plot(self):
        plt.scatter(self.X, self.y)
        Xlim = np.array([self.X.min(), self.X.max()])
        plt.plot(Xlim, self.predict(Xlim), color='Red')
```

単回帰を行うクラスです。メソッドは*scikit-learn*に倣って書きました。単回帰ということで```plot()```というメソッドも書いてみました。クラスについてよく分からない方は中の関数だけ見て下さい。

```fit()```で先ほどの計算が行なわれます。```np.var()```で分散、```np.cov()```で共分散が得られます。また、```np.cov()```は共分散行列というものを返すメソッドで、(x, y)を引数に入れると以下の様な行列（ndarray）が返されます。

```math
\begin{pmatrix}
\sigma_{xx} & \sigma_{xy} \\
\sigma_{yx} & \sigma_{yy}
\end{pmatrix}
```

$\sigma_{xx}$、$\sigma_{yy}$はxとyの分散です。欲しいのは共分散$\sigma_{xy}$なので、インデックスとして[0,1]を指定しています（[1,0]でも可）。

実行してみるとこんな感じです。

```python
# 適当なデータ生成
X = np.arange(200)
y = X + np.random.randn(200)*30

reg_simple = simple_linear()
reg_simple.fit(X, y)
reg_simple.plot()
```

![run.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/9305796d-6f0b-8bda-c683-5135484c964c.png)

## 重回帰

今度は重回帰をしましょう。単回帰と同じ様に計算を進めてもできると思いますが、計算が煩雑になるため、ある前提を考えます。それは、説明変数と目的変数が全て標準化されているというものです。標準化とは、平均が0、標準偏差が1になる様にデータをスケーリングすることです。具体的には、全データから平均値を引き、標準偏差で割ることで得られます。こうすることで切片の値を0にできるため、計算がし易くなります。

それでは計算をしていきましょう。二つの変数が標準化されていると言う前提の元で進めていきます。

まず求めたい関数を、

$f(k_1, k_2 \cdots k_m) = a_1k_1 + a_2k_2 + \cdots + a_mk_m$

と置きます。$m$は説明変数の次元数です。重回帰なので、説明変数とそれぞれの係数は多次元になります。

また、説明変数$X$、目的変数$y$を

```math
X = \begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1m} \\
x_{21} & x_{22} & \cdots & x_{2m} \\
\vdots & \vdots & \ddots & \vdots \\
x_{n1} & x_{n2} & \cdots & x_{nm} \\
\end{pmatrix} \\
```

```math
y = \begin{pmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n
\end{pmatrix}
```

とし、$X_1 \ X_2 \ \cdots \ X_n$ を

```math
X_1 = \begin{pmatrix}
x_{11} \\
x_{12} \\
\vdots \\
x_{1m} \\
\end{pmatrix}\

X_2 = \begin{pmatrix}
x_{21} \\
x_{22} \\
\vdots \\
x_{2m} \\
\end{pmatrix}\

\cdots \

X_n = \begin{pmatrix}
x_{n1} \\
x_{n2} \\
\vdots \\
x_{nm} \\
\end{pmatrix}\
```

とします。$X$は、データ数$n$ $\times$ 次元数$m$ の説明変数全体の行列で、$X_i$ は、ある目的変数$y_i$に対する説明変数を表します。

次に、最小化したい関数$E$を以下の様に定義します。

```math
E = \sum_{i=1}^{n}(f(X_i) - y_i)^2 = \sum_{i=1}^{n}(a_1x_{i1} + a_2x_{i2} + \cdots + a_mx_{im}- y_i)^2
```

そしてこの$E$を最小にする $a_1 \ a_2 \cdots \ a_m$ を求めます。$E$は各回帰係数　$a_1 \ a_2 \cdots \ a_m$　の二次関数でそれらの係数はいずれも正なので、単回帰同様、偏微分をして傾きが0になる値を探します。

```math
\begin{align}
\frac{\partial E}{\partial a_1} &= 2\sum_{i=1}^{n}x_{i1}(a_1x_{i1} + a_2x_{i2} + \cdots + a_mx_{im}- y_i) = 0 \\
\frac{\partial E}{\partial a_2} &= 2\sum_{i=1}^{n}x_{i2}(a_1x_{i1} + a_2x_{i2} + \cdots + a_mx_{im}- y_i) = 0 \\
\vdots \\
\frac{\partial E}{\partial a_m} &= 2\sum_{i=1}^{n}x_{im}(a_1x_{i1} + a_2x_{i2} + \cdots + a_mx_{im}- y_i) = 0
\end{align}
```

この連立方程式を解いていきます。こちらも単回帰同様、まず両辺を2で割り、展開してから$y_i$の項を移行します。

```math
\begin{align}
a_1\sum_{i=1}^{n}x_{i1}x_{i1} + a_2\sum_{i=1}^{n}x_{i1}x_{i2} + \cdots + a_m\sum_{i=1}^{n}x_{i1}x_{im} &= \sum_{i=1}^{n}x_{i1}y_i \\
a_1\sum_{i=1}^{n}x_{i2}x_{i1} + a_2\sum_{i=1}^{n}x_{i2}x_{i2} + \cdots + a_m\sum_{i=1}^{n}x_{i2}x_{im} &= \sum_{i=1}^{n}x_{i2}y_i \\
\vdots \\
a_1\sum_{i=1}^{n}x_{im}x_{i1} + a_2\sum_{i=1}^{n}x_{im}x_{i2} + \cdots + a_m\sum_{i=1}^{n}x_{im}x_{im} &= \sum_{i=1}^{n}x_{im}y_i
\end{align}
```

できました。ここから単回帰と手順が異なります。
この式を、**行列に置き換える**作業を行います。まずは、これらを単純な$m$行1列の行列にしてみましょう。

```math
\begin{pmatrix}
a_1\sum_{i=1}^{n}x_{i1}x_{i1} + a_2\sum_{i=1}^{n}x_{i1}x_{i2} + \cdots + a_m\sum_{i=1}^{n}x_{i1}x_{im} \\
a_1\sum_{i=1}^{n}x_{i2}x_{i1} + a_2\sum_{i=1}^{n}x_{i2}x_{i2} + \cdots + a_m\sum_{i=1}^{n}x_{i2}x_{im} \\
\vdots \\
a_1\sum_{i=1}^{n}x_{im}x_{i1} + a_2\sum_{i=1}^{n}x_{im}x_{i2} + \cdots + a_m\sum_{i=1}^{n}x_{im}x_{im}
\end{pmatrix}
=
\begin{pmatrix}
\sum_{i=1}^{n}x_{i1}y_i \\
\sum_{i=1}^{n}x_{i2}y_i \\
\vdots \\
\sum_{i=1}^{n}x_{im}y_i
\end{pmatrix}
```

次に　$a_1 \ a_2 \ \cdots \ a_m$ を縦に並べた以下の行列$a$を考えます。

```math
a = \begin{pmatrix}
a_1 \\
a_2 \\
\vdots \\
a_m
\end{pmatrix}
```

この行列を使って、左辺を表してみましょう。

```math
\begin{pmatrix}
\sum_{i=1}^{n}x_{i1}x_{i1} & \sum_{i=1}^{n}x_{i1}x_{i2} & \cdots & \sum_{i=1}^{n}x_{i1}x_{im} \\
\sum_{i=1}^{n}x_{i2}x_{i1} & \sum_{i=1}^{n}x_{i2}x_{i2} & \cdots & \sum_{i=1}^{n}x_{i2}x_{im} \\
\vdots & \vdots & \ddots & \vdots \\
\sum_{i=1}^{n}x_{im}x_{i1} & \sum_{i=1}^{n}x_{im}x_{i2} & \cdots & \sum_{i=1}^{n}x_{im}x_{im}
\end{pmatrix}
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots \\
a_m
\end{pmatrix}
=
\begin{pmatrix}
\sum_{i=1}^{n}x_{i1}y_i \\
\sum_{i=1}^{n}x_{i2}y_i \\
\vdots \\
\sum_{i=1}^{n}x_{im}y_i
\end{pmatrix}
```

こうですね。元の行列から$a$を外に出したイメージでしょうか。このまま、他の行列を$x$と$y$の比較的単純な行列で表すことができれば、$a$ を求めることができそうですね。やってみましょう。

一番左の行列は、$X^T$と$X$の積になっています。$X^T$は、$X$の**転置行列**（行と列が逆になった行列）を表します。

```math
\begin{pmatrix}
\sum_{i=1}^{n}x_{i1}x_{i1} & \sum_{i=1}^{n}x_{i1}x_{i2} & \cdots & \sum_{i=1}^{n}x_{i1}x_{im} \\
\sum_{i=1}^{n}x_{i2}x_{i1} & \sum_{i=1}^{n}x_{i2}x_{i2} & \cdots & \sum_{i=1}^{n}x_{i2}x_{im} \\
\vdots & \vdots & \ddots & \vdots \\
\sum_{i=1}^{n}x_{im}x_{i1} & \sum_{i=1}^{n}x_{im}x_{i2} & \cdots & \sum_{i=1}^{n}x_{im}x_{im}
\end{pmatrix}
=
X^TX
```

右辺の行列は、$X^T$と$y$の積で表せますね。

```math
\begin{pmatrix}
\sum_{i=1}^{n}x_{i1}y_i \\
\sum_{i=1}^{n}x_{i2}y_i \\
\vdots \\
\sum_{i=1}^{n}x_{im}y_i
\end{pmatrix}
=
X^Ty
```

まとめてみましょう

```math
X^TXa = X^Ty
```

綺麗にまとまりました。
$a=$ 〜　の形にしたいので、両辺に$(X^TX)^{-1}$を左からかけ、左辺を$a$にしましょう。

```math
a = (X^TX)^{-1}X^Ty
```

<br>
できました。
$(X^TX)^{-1}$は$X^TX$の**逆行列**を表します。逆行列とは、積をとった時に**単位行列**となる行列で、単位行列とは、積をとった時にどの行列に対しても変化を与えない行列を指します。行列において、逆行列は「逆数」、単位行列は「1」の様な働きします。
<br>

### Pythonで実装
では、ここまでをPythonで実装してみましょう。

```python
import numpy as np
from scipy.stats import zscore # 標準化用の関数

class linear_sta:
    def fit(self, X, y):
        # 標準偏差と平均の保存
        self.X_mean = np.array(X.mean(axis=0)) # 平均
        self.X_std = np.array(X.std(axis=0)) # 標準偏差
        self.y_mean = y.mean()
        self.y_std = y.std()
        
        # 標準化して行列に変換
        self.X = np.matrix(zscore(X, axis=0))
        self.XT = np.matrix(zscore(X.T, axis=1))
        self.y = np.matrix(zscore(np.array(y)).reshape(-1,1))
         
        # 回帰係数
        XT_X_inverse = np.linalg.pinv(self.XT * self.X)  # XTとXの積の逆行列
        XT_y = self.XT * self.y  # XTとyの積
        self.a = np.array(XT_X_inverse * XT_y).ravel()　# 積をとって平坦化
    
    # 予測
    def predict(self, X):
        sta_X = np.array((X - self.X_mean) / self.X_std) # fitで入れた説明変数でスケーリング
        return np.dot(sta_X, self.a) * self.y_std + self.y_mean # 回帰係数との内積をとって標準化した分を戻す
```

例によって```fit()```で計算が行われます。一つ一つ見ていきましょう。

```python
# 標準偏差と平均の保存
self.X_mean = np.array(X.mean(axis=0)) # 平均
self.X_std = np.array(X.std(axis=0)) # 標準偏差
self.y_mean = y.mean()
self.y_std = y.std()
```

平均値と標準偏差を変数に格納しています。Xとyが標準化されているという前提のもとで回帰係数($a$)を求めたため、予測を行う際に、入力データをこのデータに基づいてスケーリングする必要があるからです。また、```axis```のデフォルト値は*None*なので、「列ごとに」を意味する```axis=0```を指定しています。

```python
# 標準化して行列に変換
self.X = np.matrix(zscore(X, axis=0))
self.XT = np.matrix(zscore(X.T, axis=1))
self.y = np.matrix(zscore(np.array(y)).reshape(-1,1))
```

標準化したデータを行列として変数に代入しています。```zscore()```で標準化、```np.matrix```で行列に変換しています。```np.array()```でndarrayに変換しても同じ演算ができますが、行列として扱う場合は```matrix```の方が便利です。```.T```は転置行列、```reshape(-1,1)```は 1 $\times$ $n$ の配列を $n$ $\times$ $1$ に変換して返しています。

```python
# 回帰係数
XT_X_inverse = np.linalg.pinv(self.XT * self.X)  # XTとXの積の逆行列
XT_y = self.XT * self.y  # XTとyの積
self.a = np.array(XT_X_inverse * XT_y).ravel() # 積をとって平坦化
```

定義した行列を使って回帰係数を求めています。```np.linalg.pinv()```は逆行列を返します。```np.linalg.inv()```または```**-1```でも逆行列が求められますが、元の行列が正則行列でない（逆行列を持たない）ときにエラーを吐いてしまいます。```np.linalg.pinv()```では擬似的な逆行列を生成してエラーを回避できるため、採用しています。

```predict()```も見ていきましょう

```python
    # 予測
    def predict(self, X):
        sta_X = np.array((X - self.X_mean) / self.X_std) # fitで入れた説明変数でスケーリング
        return np.dot(sta_X, self.a) * self.y_std + self.y_mean # 回帰係数との内積をとって標準化した分を戻す
```

まず、入力データから```fit()```で定義した平均を引き、標準偏差で割ることでスケーリングをしています。そして```np.dot()```で回帰係数との内積をとり、目的変数の標準偏差をかけて平均を足すことで、スケーリングした分を元に戻しています。

### 予測してみる

実際に予測してみましょう。[こちら](https://archive.ics.uci.edu/ml/datasets/Computer+Hardware)のデータセットを使います。コンクリートの成分から圧縮強度ってやつを求めてみましょう。

```python
import pandas as pd
from sklearn.model_selection import train_test_split

data_path = 'https://archive.ics.uci.edu/ml/machine-learning-databases/concrete/compressive/Concrete_Data.xls'
df = pd.read_excel(data_path) # データセット読み込み

y_columns = 'Concrete compressive strength(MPa, megapascals) ' # 圧縮強度
X = df.drop(columns=y_columns) # 説明変数
y = df[y_columns] # 目的変数

X_train, X_test, y_train, y_test = train_test_split(X, y) # 訓練データとテストデータに分割

linear_s = linear_sta() # インスタンス生成
linear_s.fit(X_train, y_train) # 回帰係数の計算

pre_X = X_test.iloc[:3] # テストデータの上三つ
pre_y = np.array(y_test[:3]) # その正解
print(f'予測値：　{linear_s.predict(pre_X)}')
print(f'正解：　{pre_y}')
```

実行するとこうなりました。

```
予測値：　[40.02633476 25.91338266 31.51666727]
正解：　[40.8686899  19.41564416 38.56239268]
```

まあまあですかね？　```train_test_split()```がランダムなので、試行毎に結果は変わります。

### 回帰係数を直す

今求めた回帰係数を見てみましょう。

```python
print(linear_s.a)
```

するとこうなります。

```
[ 0.77131819  0.55916812  0.37295602 -0.15677616  0.10616719  0.09854152
  0.11658664  0.43031912]
```

0付近に固まっていますが、これはデータを標準化したためです。上でやった通りこれでも予測はできますが、回帰係数は直線の傾きを表しており、各説明変数が目的変数に与えている影響として見ることができるため、可能なら元のスケールで表したいですね。

```predict()```の中を数式で表してみましょう。入力データを縦に並べた行列を$k$、$X$の列ごとの標準偏差・平均を縦に並べた行列をそれぞれ $\sigma_x, \bar{X}$ とすると、出力を以下の様に表せます。

```math
(k - \bar{X}) \oslash \sigma_x ・ a \times \sigma_y + \bar{y}
```

$\oslash$はアダマール積の記号です。行列の格要素ごとに除算を行います。乗算の場合は$\otimes%$を使います。

上の式を、$k \ ・A + B$ の形に変形します。

```math
\begin{align}
(k - \bar{X}) \ ・ \frac{a}{\sigma_x}\sigma_y + \bar{y} \\
k \ ・ \frac{a}{\sigma_x}\sigma_y - \bar{X} \ ・ \frac{a}{\sigma_x}\sigma_y + \bar{y}
\end{align}
```

こうなりました。これで、元の回帰係数 $A = \frac{a}{\sigma_x}\sigma_y$、切片 $B = - \bar{X} \ ・ A + \bar{y}$ が求まりました。やったぜ。

### Pythonで実装

以上の内容を、Pythonにも反映させましょう。

```python
class Linear:
    def fit(self, X, y):
        # 標準偏差と平均
        self.X_mean = np.array(X.mean(axis=0))
        self.X_std = np.array(X.std(axis=0))
        
        # 標準化して行列に変換
        self.X = np.matrix(zscore(X, axis=0))
        self.XT = np.matrix(zscore(X.T, axis=1))
        self.y = np.matrix(zscore(np.array(y)).reshape(-1,1))
         
        # 回帰係数と切片
        XT_X_inverse = np.linalg.pinv(self.XT * self.X)  # XTとXの積の逆行列
        XT_y = self.XT * self.y  # XTとyの積
        self.A = np.array(XT_X_inverse * XT_y).ravel() / self.X_std * y.std()
        self.B = -np.dot(self.b, self.X_mean) + y.mean()
    
    # 予測
    def predict(self, X):
        return np.dot(X, self.A) + self.B
```

```predict()``` がシンプルになりましたね。予測もしてみましょう。

```python
linear = Linear()
linear.fit(X_train, y_train)
print(f'予測値：　{linear_s.predict(pre_X)}')
```

```
予測値：　[40.02633476 25.91338266 31.51666727]
```

しっかり同じ値になりました。

## おまけ

*scikit-learn*を使って先ほど求めた回帰係数が正しいか確かめてみましょう。

```python
from sklearn import linear_model

reg = linear_model.LinearRegression()
reg.fit(X_train, y_train)

print('自作')
print(linear.A)
print(linear.B)
print('----------------------')
print('sklearn')
print(reg.coef_)
print(reg.intercept_)
```

結果

```
自作
[ 0.12291481  0.10560904  0.09741413 -0.12366291  0.29375934  0.02143313
  0.02399962  0.11264403]
-35.85204847096638
----------------------
sklearn
[ 0.12291481  0.10560904  0.09741413 -0.12366291  0.29375934  0.02143313
  0.02399962  0.11264403]
-35.852048470969976
```

ちゃんと一致していました

これで終わります。ありがとうございました!

参考
- [【大学数学】最小二乗法(回帰分析)【確率統計】](https://www.youtube.com/watch?v=Zz1sgYxrA-k)
- [最小二乗法による線形重回帰分析](https://datachemeng.com/wp-content/uploads/ordinaryleastsquares.pdf)