徹底理解：単回帰解析

# はじめに
Pythonには、scikit-learnなど便利なライブラリが用意されており、比較的簡単に解析プログラムを組むことができます。ライブラリの中身はブラックボックスになりがちです。
本稿では、単回帰解析に焦点を当て、そのブラックボックスの中を覗いてみたいと思います。
使用するライブラリは、おもにnumpyとTensorFlowとなります。
本稿を書くに当たり、下記書籍を参考にしました。

* [TensorFlow機械学習クックブック Pythonベースの活用レシピ60+](https://book.impress.co.jp/books/1116101172)
* [統計的機械学習の数理100問 with Python](https://www.kyoritsu-pub.co.jp/bookdetail/9784320125070)

TensorFlow機械学習クックブック Pythonベースの活用レシピ60+は好著なのですが、TensorFlow1系でプログラムが書かれています。せっかくなので、ここではTensorFlow2系でプログラミングしてみたいと思います。

# 単回帰解析
説明変数 $(x_i, y_i)$ にできるだけ適合する回帰直線 $y = Ax + b$ を求めます。
![FireShot Capture 008 - 2-1copied.ipynb - Colaboratory - colab.research.google.com.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/c20e8901-fc27-5abe-203d-06c51a417075.png)

## 最小二乗法(least squares method)
説明変数 $x_1, \cdots, x_N, y_1, \cdots, y_N \in \mathbb{R}$ の直線からの距離 $|y_i - b - Ax_i|$ の二乗の和 $L:=\sum_{i=1}^{N}(y_i - b - Ax_i)^2$ を最小にする回帰直線の切片 $b$ と傾き $A$ を求めます。

$L$ を $A, b$ で偏微分して0とおくと、下記2式が得られます。

```math
\begin{eqnarray}
\frac{{\partial}L}{{\partial}b} &=& -2 \sum_{i=1}^{N}(y_i - b - Ax_i) = 0 \tag{1} \\
\frac{{\partial}L}{{\partial}A} &=& -2 \sum_{i=1}^{N}x_i(y_i - b - Ax_i) = 0 \tag{2}
\end{eqnarray}
```

この連立方程式を解きます。  
(1)(2)を解くと、$\sum_{i=1}^{N}(x_i-\bar{x})^w \neq 0$, すなわち　$x_1 = \cdots = x_N ではない$　とき、

```math
\begin{eqnarray}
\hat{A} &=& \frac{\sum_{i=1}^{N}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{N}(x_i - \bar{x})^2} \tag{3} \\
\hat{b} &=& \bar{y} - \hat{A}\bar{x} \tag{4}
\end{eqnarray}
```

が得られると、[統計的機械学習の数理100問 with Python](https://www.kyoritsu-pub.co.jp/bookdetail/9784320125070)に解説が載っています。  
この連立方程式から解が得る過程をもう少し詳しく見てみます。こちらの[ページ（単回帰分析・最小二乗法の公式はどうすれば求められるのか。統計上の誤差と残差の違い）](https://atarimae.biz/archives/12371)を参考にしました。

```math
\begin{eqnarray}
0 &=& \sum_{i=1}^{n}y_i - \hat{A}\sum_{i=1}^{n}x_i - \hat{b} \cdot n \\
\hat{b} &=& \frac{1}{n}\sum_{i=1}^{n}y_i - \frac{\hat{A}}{n}\sum_{i = 1}^{n}x_i \\
&=& \bar{y} - \hat{A} \cdot \bar{x} \\
\\
0 &=& \sum_{i=1}^{n}x_{i}y_{i} - \hat{A}\sum_{i=1}^{n}x_i^2 - \hat{b}\sum_{i=1}^{n}x_i \\
&=& \sum_{i=1}^{n}{x_i}{y_i} - \bar{y}\sum_{i=1}^{n}x_i - \hat{A} \left( \sum_{i=1}^{n}x_i^2 - \bar{x}\sum_{i=1}^{n}x_i \right) \\
&=& \sum_{i=1}^{n}{x_i}{y_i} - n \bar{x} \bar{y} - \hat{A} \left( \sum_{i=1}^{n}x_i^2 - n \bar{x}^2 \right) \\
\hat{A} &=& \frac{\sum_{i=1}^{n}{x_i}{y_i} - n \bar{x} \bar{y}}{\sum_{i=1}^{n}x_i^2 - n\bar{x^2}} \\
&=& \frac{\sum_{i=1}^{n}{(x_iy_i - \bar{x}\bar{y})}}{\sum_{i=1}^{n}{(x_i^2 - \bar{x}^2)}}
\end{eqnarray}
```

先ほどの解 $\hat{A} = \frac{\sum_{i=1}^{N}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{N}(x_i - \bar{x})^2}$ と $\hat{A} = \frac{\sum_{i=1}^{n}{(x_iy_i - \bar{x}\bar{y})}}{\sum_{i=1}^{n}{(x_i^2 - \bar{x}^2)}}$ を確認してみます。分母に着目すると、$\sum_{i=1}^{N}(x_i - \bar{x})^2 =\sum_{i=1}^{n}{(x_i^2 - \bar{x}^2)}$ ということになります。$(x_i - \bar{x})^2$ を展開すると、$x_i^2 - 2x_i\bar{x} + \bar{x}^2$ となるので、腑に落ちません。$\bar{x}$ が平均値 $\frac{1}{n}{\sum_{i=1}^{n}{x_i}}$ だからだろうと推測しましたが、せっかくなので数学的に確認してみました。

まず、$x_0, x_1, \cdots, x_n$ の合計を $T = \sum_{i=1}^{n}{x_i}$ とおいておきます。

```math
\begin{eqnarray}
\sum_{i=1}^{n}{(x_i^2 - \bar{x}^2)} &=& \frac{\sum_{i=1}^{n}{(n^2x_i^2 - T^2)}}{n^2} \tag{5}\\
\sum_{i=1}^{n}(x_i - \bar{x})^2 &=& \frac{\sum_{i=1}^{n}{(nx_i - T)^2}}{n^2} \tag{6}
\end{eqnarray}
```

(5), (6)より

```math
\begin{eqnarray}
\sum_{i=1}^{n}{(n^2x^2_i - T^2)} &=& \sum_{i=1}^{n}{(nx_i - T)^2} \\
&=& \sum_{i=1}^{n}{(n^2x_i^2 - 2nx_iT + T^2)}
\end{eqnarray}
```
が得られます。

```math
\begin{eqnarray}
-\sum_{i=1}^{n}{T^2} &=& \sum_{i=1}^{n}{(T^2 - 2nx_iT)} \\
-nT^2 &=& nT^2 - \sum_{i=1}^{n}{2nx_iT} \\
-2nT^2 &=& -2nT\sum_{i=1}^{n}{x_i} \\
T &=& \sum_{i=1}^{n}{x_i}
\end{eqnarray}
```

$T=\sum_{i=1}^{n}{x_i}$であるから、両辺が等価であることが分かります。  
少し回り道をしましたが、早速(3)、(4)の式を使って、最小二乗法を実装してみます。

## プログラム
```Python
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf

# データの作成
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100)
y_vals = x_vals + np.random.normal(0, 1, 100)
```

```Python
def min_eq(x, y):
    x_bar, y_bar = np.mean(x), np.mean(y)
    beta_1 = tf.divide(np.dot(tf.subtract(x, x_bar), tf.subtract(y, y_bar)), tf.square(tf.linalg.norm(tf.subtract(x, x_bar))))
    beta_0 = tf.subtract(y_bar, tf.multiply(beta_1, x_bar))
    return beta_1, beta_0
```

```Python
A, b = min_eq(x_vals, y_vals)
print('slope: {}'.format(A.numpy()))
print('y_intercept: {}'.format(b.numpy()))
```

> slope: 0.9931105078364767
y_intercept: 0.06955974394305287

```Python
# 最も適合する直線を取得
best_fit = []
for i in x_vals:
    best_fit.append(A * i + b)

# 結果をプロット
plt.plot(x_vals, y_vals, 'o', label='Data')
plt.plot(x_vals, best_fit, 'r-', label='Best fit line', linewidth=3)
plt.legend(loc='upper left')
plt.show()
```

![3-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/b252734c-b362-45a8-539c-4a809cc82dfb.png)


# 逆行列法
線形回帰(linear regression)は、$A \cdot x = b$ といった一連の行列方程式として表すことができます。行列 $x$ の係数を求めます。$x$ を解くための式は、$x = (A^T \cdot A)^{-1} \cdot A^T \cdot b$ として表すことがでます。

もう少し詳しく見てみます。

直線の方程式は次のように表される。

```math
y = ax + b
```

この直線が点 $(x_1, y_1)$ と点 $(x_2, y_2)$ を通過するならば、

```math
y_1 = ax_1 + b \\
y_2 = ax_2 + b
```

が成立する。これを行列で書くと

```math
\begin{bmatrix}
y_1 \\
y_2
\end{bmatrix}
=
\begin{bmatrix}
x_1 & 1 \\
x_2 & 1
\end{bmatrix}
\begin{bmatrix}
a \\
b
\end{bmatrix}
```

となる。各行列を簡略化して表すと

$$
Y = XA
$$

3点以上が与えられた場合、各項は下記のようになります。

```math
Y =
\begin{bmatrix}
y_1 \\
y_2 \\
y_3
\end{bmatrix} \tag{1}
```

```math
A =
\begin{bmatrix}
a \\
b
\end{bmatrix} \tag{2}
```

```math
X =
\begin{bmatrix}
x_1 && 1\\
x_2 && 1\\
x_3 && 1
\end{bmatrix} \tag{3}
```

$X$ は正方行列ではなく、逆行列を持てません。  
そこで、$A$ の前にある行列項を、行列演算して、正方行列に変換します。  
$X$ の転置行列 $X^T$ を乗じます。

```math
X^T Y = X^T X A \tag{4}
```

$X^T X$は正方行列になるので、これが正則行列であるならば、逆行列 $(X^T X)^{-1}$ が存在します。行列式を式(4)の左から掛けると

```math
(X^T X)^{-1} X^T Y = (X^T X)^{-1} X^T X A \tag{5}
```

右辺Aの前の行列演算結果は単位行列 $E$ となるので、次式が得られます。

```math
(X^T X)^{-1} X^T Y = A \tag{6}
```

このようにして一組のAが得られます。

## プログラム
```Python
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf

# データの作成
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100)
y_vals = x_vals + np.random.normal(0, 1, 100)
```

```Python
# 計画行列Aを作成
x_vals_column = np.transpose(np.matrix(x_vals))
ones_column = np.transpose(np.matrix(np.repeat(1, 100)))
A = np.column_stack((x_vals_column, ones_column))

# 行列bを作成
b = np.transpose(np.matrix(y_vals))
```

```Python
# テンソルを作成
A_tensor = tf.constant(A)
b_tensor = tf.constant(b)
```

```Python
# 逆行列法を適用
tA_A = tf.matmul(tf.transpose(A_tensor), A_tensor)
tA_A_inv = tf.linalg.inv(tA_A)
product = tf.matmul(tA_A_inv, tf.transpose(A_tensor))
solution_eval = tf.matmul(product, b_tensor)
```

```Python
# 係数を抽出
slope = solution_eval[0][0]
y_intercept = solution_eval[1][0]
print('slope: ' + str(slope))
print('y_intercept: ' + str(y_intercept))
```

> slope: tf.Tensor(0.9931105078364761, shape=(), dtype=float64)
y_intercept: tf.Tensor(0.06955974394305653, shape=(), dtype=float64)

```Python
# 最も適合する直線を取得
best_fit = []
for i in x_vals:
    best_fit.append(slope * i + y_intercept)

# 結果をプロット
plt.plot(x_vals, y_vals, 'o', label='Data')
plt.plot(x_vals, best_fit, 'r-', label='Best fit line', linewidth=3)
plt.legend(loc='upper left')
plt.show()
```

![3-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/5448764e-e887-d9a2-62fb-e2b75845716d.png)

# 行列分解法
逆行列法の実装は、行列が大きくなると計算効率が落ちる場合があります。行列 $A$ を分解し、分解した行列で行列演算を行う方法を実装します。TensorFlowに実装されているコレスキー分解関数を利用します。ここでは、連立方程式 $A \cdot x = b$ を $L \cdot L' \cdot x = b$ として解いてみよう。まず $L \cdot y = b$ を解き、次に $L' \cdot x = y$ を解くことで、係数行列 $x$ を求めます。

```Python
import matplotlib.pyplot as pyplot
import numpy as np
import tensorflow as tf
from tensorflow.python.framework import ops

ops.reset_default_graph()

# データを作成
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100)
y_vals = x_vals + np.random.normal(0, 1, 100)

# 計画行列Aを作成
x_vals_column = np.transpose(np.matrix(x_vals))
ones_column = np.transpose(np.matrix(np.repeat(1, 100)))
A = np.column_stack((x_vals_column, ones_column))

# 行列bを作成
b = np.transpose(np.matrix(y_vals))

# テンソルを作成
A_tensor = tf.constant(A)
b_tensor = tf.constant(b)
```

正方行列のコレスキー分解を実行します。

```Python
# コレスキー分解
tA_A = tf.matmul(tf.transpose(A_tensor), A_tensor)
L = tf.linalg.cholesky(tA_A)

# L * y = t(A) * b を解く
tA_b = tf.matmul(tf.transpose(A_tensor), b)
sol1 = tf.linalg.solve(L, tA_b)

# L' * y = sol1 を解く
sol2 = tf.linalg.solve(tf.transpose(L), sol1)
solution_eval = sol2
```

```Python
# 係数を抽出
slope = solution_eval[0][0]
y_intercept = solution_eval[1][0]
print('slope: {}'.format(slope))
print('y_intercept: {}'.format(y_intercept))
```

> slope: 0.9931105078364756
y_intercept: 0.06955974394305893

```Python
# 最も適合する直線を取得
best_fit = []
for i in x_vals:
    best_fit.append(slope + i + y_intercept)

# 結果をプロット
plt.plot(x_vals, y_vals, 'o', label='Data')
plt.plot(x_vals, best_fit, 'r-', label=' Best fit line', linewidth=3)
plt.legend(loc='upper left')
plt.show()
```

![3-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/97d81e69-8834-409a-2d06-241da8ef80a5.png)


# TensorFlowでの線形回帰の実装パターン

```Python
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from tensorflow.python.framework import ops

ops.reset_default_graph()

# データを作成
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100).astype(np.float32)
y_vals = x_vals + np.random.normal(0, 1, 100).astype(np.float32)
```

```Python
# 学習率とバッチサイズを設定
learning_rate = 0.01
batch_size = 25

A = tf.Variable(0.)
b = tf.Variable(0.)
```

線形モデルの方程式  y=Ax+b  を設定する。

```Python
def f(x):
    y = A * x + b
    return y
```

```Python
def loss(y_pred, y_true):
    return tf.reduce_mean(tf.square(y_pred - y_true))
```

```Python
loss_vec = []
for i in range(100):
    with tf.GradientTape() as tape:
        y_predicted = f(x_vals)
        tmp_loss = loss(y_predicted, y_vals)
    # get gradients
    gradients = tape.gradient(tmp_loss, [A, b])
    loss_vec.append(tmp_loss)
    # compute and adjust weights
    A.assign_sub(gradients[0] * learning_rate)
    b.assign_sub(gradients[1] * learning_rate)

    if (i + 1) % 25 == 0:
        print('Step #{} A = {} b = {}'.format(i + 1, A.numpy(), b.numpy()))
        print('Loss = {}'.format(tmp_loss))
```

> Step #25 A = 0.9827542901039124 b = 0.13843616843223572
Loss = 0.9922008514404297
Step #50 A = 0.9839657545089722 b = 0.13037890195846558
Loss = 0.9919329285621643
Step #75 A = 0.9850354790687561 b = 0.12326423078775406
Loss = 0.9917240738868713
Step #100 A = 0.9859801530838013 b = 0.11698180437088013
Loss = 0.9915611147880554

```Python
plt.scatter(x_vals, y_vals)
# plot the best fit line
plt.plot(x_vals, f(x_vals), 'r')
plt.show()

plt.plot(loss_vec, 'k-')
plt.title('L2 Loss per Generation')
plt.xlabel('Generation')
plt.ylabel('L2 Loss')
plt.show()
```

![3-4-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/10d99303-6a70-536b-065c-f5f84d1830e9.png)
![3-4-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/8031c89a-da1b-01ab-d0fd-b9835e172d05.png)

# 線形回帰の損失関数

```Python
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from sklearn import datasets
from tensorflow.python.framework import ops

ops.reset_default_graph()

# データを作成
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100).astype(np.float32)
y_vals = x_vals + np.random.normal(0, 1, 100).astype(np.float32)

# バッチサイズ、学習率、繰り返し回数を設定
batch_size = 25
learning_rate = 0.1 # 学習率が0.4を超えると収束しなくなる。
iterations = 50

# 変数、モデルの演算を定義
A = tf.Variable(tf.random.normal(shape=[1, 1]))
b = tf.Variable(tf.random.normal(shape=[1, 1]))

def model_output(x):
    return tf.add(tf.matmul(x, A), b)
```

損失関数をL1損失関数に変更します。

```Python
def loss_l1(model, x, y):
    return tf.reduce_mean(tf.abs(y - model(x)))
```

```Python
# 最適化関数を指定
my_opt_l1 = tf.keras.optimizers.SGD(learning_rate=learning_rate)

# トレーニンググループを開始
loss_vec_l1 = []
for i in range(iterations):
    rand_index = np.random.choice(len(x_vals), size=batch_size)
    rand_x = tf.Variable(np.transpose([x_vals[rand_index]]))
    rand_y = tf.Variable(np.transpose([y_vals[rand_index]]))

    with tf.GradientTape() as tape:
        temp_loss_l1 = loss_l1(model_output, rand_x, rand_y)
    grads = tape.gradient(temp_loss_l1, [A, b])
    loss_vec_l1.append(temp_loss_l1)
    my_opt_l1.apply_gradients(zip(grads, [A, b]))
    
    if (i + 1) % 25 == 0:
        print('Step #{} A = {} b = {}'.format(i + 1, A.numpy(), b.numpy()))
```

> Step #25 A = [[0.96738863]] b = [[0.88311166]]
Step #50 A = [[0.780318]] b = [[0.5991116]]

```Python
plt.scatter(x_vals, y_vals)
# plot the best fit line
plt.plot(x_vals, model_output(tf.Variable(np.transpose([x_vals]))), 'r')
plt.show()
```

![3-5-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/c8ecf7fc-c2b0-6512-b183-ad1107f95958.png)

L2損失関数を使用するコード

```Python
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100).astype(np.float32)
y_vals = x_vals + np.random.normal(0, 1, 100).astype(np.float32)

# バッチサイズ、学習率、繰り返し回数を設定
batch_size = 25
learning_rate = 0.01 # 学習率が0.4を超えると収束しなくなる。
iterations = 50

# 変数、モデルの演算を定義
A = tf.Variable(tf.random.normal(shape=[1, 1]))
b = tf.Variable(tf.random.normal(shape=[1, 1]))

def model_output(x):
    return tf.add(tf.matmul(x, A), b)

def loss_l2(model, x, y):
    return tf.reduce_mean(tf.square(y - model(x)))

# 最適化関数を指定
my_opt_l2 = tf.keras.optimizers.SGD(learning_rate=learning_rate)

# トレーニンググループを開始
loss_vec_l2 = []
for i in range(iterations):
    rand_index = np.random.choice(len(x_vals), size=batch_size)
    rand_x = tf.Variable(np.transpose([x_vals[rand_index]]))
    rand_y = tf.Variable(np.transpose([y_vals[rand_index]]))

    with tf.GradientTape() as tape:
        temp_loss_l2 = loss_l2(model_output, rand_x, rand_y)
    grads = tape.gradient(temp_loss_l2, [A, b])
    loss_vec_l2.append(temp_loss_l2)
    my_opt_l2.apply_gradients(zip(grads, [A, b]))
    
    if (i + 1) % 25 == 0:
        print('Step #{} A = {} b = {}'.format(i + 1, A.numpy(), b.numpy()))
```

> Step #25 A = [[1.061035]] b = [[-0.483181]]
Step #50 A = [[1.0561559]] b = [[-0.41459158]]

```Python
plt.scatter(x_vals, y_vals)
# plot the best fit line
plt.plot(x_vals, model_output(tf.Variable(np.transpose([x_vals]))), 'r')
plt.show()
```

![3-5-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/7ec897e6-7674-4f50-f585-f6c24a10ca20.png)


```Python
plt.plot(loss_vec_l1, 'k-', label='L1 Loss')
plt.plot(loss_vec_l2, 'r--', label='L2 Loss')
plt.title('L1 and L2 Loss per Generattion')
plt.xlabel('Generation')
plt.ylabel('L1 loss')
plt.legend(loc='upper right')
plt.show()
```

![3-5-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/e8faca91-0918-4013-b3a4-1731bae77401.png)

# デミング回帰
最小二乗線形回帰（least squares linear regression）が直線に対する線の距離を最小化するのに対し、デミング回帰は直線に対する全距離を最小化します。

```Python
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from sklearn import datasets
from tensorflow.python.framework import ops

ops.reset_default_graph()

# データを作成
np.random.seed(1234)
x_vals = np.linspace(0, 10, 100).astype(np.float32)
y_vals = x_vals + np.random.normal(0, 1, 100).astype(np.float32)

# バッチサイズ、学習率、繰り返し回数を設定
batch_size = 25
learning_rate = 0.05 # 学習率が0.4を超えると収束しなくなる。
iterations = 50

# 変数、モデルの演算を定義
A = tf.Variable(tf.random.normal(shape=[1, 1]))
b = tf.Variable(tf.random.normal(shape=[1, 1]))

def model_output(x):
    return tf.add(tf.matmul(x, A), b)
```

直線を $y = mx + b$ 、データ点を $(x_0, y_0)$ とすると、これらの点の垂直距離を次のように定義できます。

```math
d = \frac{|y_0 - (mx_0 + b)|}{\sqrt{m^2 +1}}
```

```Python
def loss(model, x, y):
    deming_numerator = tf.abs(tf.subtract(y, tf.add(tf.matmul(x, A), b)))
    deming_denominator = tf.sqrt(tf.add(tf.square(A), 1))
    return tf.reduce_mean(tf.truediv(deming_numerator, deming_denominator))
```

```Python
# 最適化関数を指定
my_opt = tf.keras.optimizers.SGD(learning_rate=0.05)

# トレーニンググループを開始
loss_vec = []
for i in range(1500):
    rand_index = np.random.choice(len(x_vals), size=batch_size)
    rand_x = tf.Variable(np.transpose([x_vals[rand_index]]))
    rand_y = tf.Variable(np.transpose([y_vals[rand_index]]))

    with tf.GradientTape() as tape:
        temp_loss = loss(model_output, rand_x, rand_y)
    grads = tape.gradient(temp_loss, [A, b])
    loss_vec.append(temp_loss)
    my_opt.apply_gradients(zip(grads, [A, b]))
    
    if (i + 1) % 100 == 0:
        print('Step #{} A = {} b = {}'.format(i + 1, A.numpy(), b.numpy()))
```

> Step #100 A = [[1.0606606]] b = [[-0.02573951]]
Step #200 A = [[1.0778327]] b = [[-0.05546956]]
Step #300 A = [[1.0786034]] b = [[-0.11129877]]
Step #400 A = [[1.0735927]] b = [[-0.13460638]]
Step #500 A = [[0.9873692]] b = [[-0.10372327]]
Step #600 A = [[1.0768895]] b = [[-0.12781776]]
Step #700 A = [[0.9659874]] b = [[-0.15313245]]
Step #800 A = [[1.1023445]] b = [[-0.15335608]]
Step #900 A = [[1.0503265]] b = [[-0.17138588]]
Step #1000 A = [[1.0963748]] b = [[-0.20842536]]
Step #1100 A = [[1.060973]] b = [[-0.2865414]]
Step #1200 A = [[1.0580977]] b = [[-0.24569693]]
Step #1300 A = [[1.0354245]] b = [[-0.1955461]]
Step #1400 A = [[1.0163252]] b = [[-0.19011997]]
Step #1500 A = [[1.1161515]] b = [[-0.20182066]]

```Python
# 係数を抽出
slope = A[0][0]
y_intercept = b[0][0]

# 最も適合する直線を取得
best_fit = []
for i in x_vals:
    best_fit.append(slope + i + y_intercept)

# 結果をプロット
plt.plot(x_vals, y_vals, 'o', label='Data')
plt.plot(x_vals, best_fit, 'r-', label=' Best fit line', linewidth=3)
plt.legend(loc='upper left')
plt.show()
```

![3-6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/268140/5529ee9c-7b39-a98b-ee71-6aff4233d8ba.png)

# まとめ
TensorFlowとnumpyで単回帰解析を実装することにより、数学的な定義に基いたロジックを理解することができました。今後は、重回帰解析、サポートベクトルマシン、ニューラルネットワークなどをTensorFlowとnumpyで実装し、数学的に理解を深めていきたいと思います。
最後に、数学は楽しいです。あと、プログラマーの性か、$\sum_{i=1}^{n}$ と書くべきところをついつい $\sum_{i=0}^{n}$ と書いてしまいます。