1次元入力の直線モデル

# 1次元入力の直線モデル
機械学習における区分けは「教師あり」と「教師なし」がある。
その中でも「教師あり」は、さらに、**回帰**と**分類**に問題を分けることができる。

**回帰**は、入力に対して連続した数値を対応付ける問題。
**分類**は、入力に対して順番のないラベルを対応付ける問題。

1次元入力の直線モデルでは回帰について取り扱う。

年齢x, 身長tについてデータをセットする。
```python
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(seed=1)              # 乱数を固定
X_min = 4                           # Xの下限値
X_max = 30                          # Xの上限値
X_n = 16                            # データ個数
X = 5 + 25 * np.random.rand(X_n)
Prm_c = [170, 108, 0.2]             # 生成パラメータ

T = Prm_c[0] - Prm_c[1] * np.exp(-Prm_c[2] * X) + 4 * np.random.randn(X_n)
np.savez('ch5_data.npz', X=X, X_min=X_min, X_max=X_max, X_n=X_n, T=T)

print(np.round(X, 2)) # 小数点以下を四捨五入
print(np.round(T, 2))

plt.figure(figsize=(4, 4))
plt.plot(X, T, marker='o', linestyle='None', markeredgecolor='black', color='cornflowerblue')
plt.xlim(X_min, X_max)
plt.grid(True)
plt.show()
```
### 実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/f01d27c3-ecba-0f8f-8326-a54d6a9bd74c.png)
```
[15.43 23.01  5.   12.56  8.67  7.31  9.66 13.64 14.92 18.47 15.48 22.13
 10.11 26.95  5.68 21.76]
[170.91 160.68 129.   159.7  155.46 140.56 153.65 159.43 164.7  169.65
 160.71 173.29 159.31 171.52 138.96 165.87]
 ```

このとき、xを入力変数、tを目標変数と呼ぶ。
目的としては、データベースにない人の年齢xに対して、その人の身長tを予測すること。

## 直線モデル
データのばらつきに対し、ある程度の誤差を許せば、データ上に直線を引くことができる。
その直線を使って、入力値xから目標変数tを予測することができる。

直線の式は
```
y(x) = w0 * x + w1
```
で表せ、直線モデルと呼ぶ。

データに合う直線を求めるためのw0, w1を決める方法について求めていく。

## 二乗誤差関数
「データに合う」目安として、データ数Nとして誤差Jを定義する。
```
J = 1/N * Σ(yn - tn)^2
```

このときにynに対しての直線モデルにxnを入れた場合は
```
yn = y(xn) = w0 * xn + w1
```
となる。

Jは平均二乗誤差と呼ばれ、直線とデータ点の差の平均になる。

w0とw1が決まると、直線ynが求められ、データ点tnと比較することができ、平均二乗誤差Jを計算することができる。

w0とw1の碁盤目の点でJの値を計算し、グラフを出力する。

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(seed=1)              # 乱数を固定
X_min = 4                           # Xの下限値
X_max = 30                          # Xの上限値
X_n = 16                            # データ個数
X = 5 + 25 * np.random.rand(X_n)
Prm_c = [170, 108, 0.2]             # 生成パラメータ

T = Prm_c[0] - Prm_c[1] * np.exp(-Prm_c[2] * X) + 4 * np.random.randn(X_n)

# 平均誤差関数
def mse_line(x, t, w):
    y = w[0] * x + w[1]
    mse = np.mean((y - t)**2)
    return mse

# 計算
xn = 100    # 等高線表示の解像度
w0_range = [-25, 25]
w1_range = [120, 170]
x0 = np.linspace(w0_range[0], w0_range[1], xn)
x1 = np.linspace(w1_range[0], w1_range[1], xn)

xx0, xx1 = np.meshgrid(x0, x1)

J = np.zeros((len(x0), len(x1)))

for i0 in range(xn):
    for i1 in range(xn):
        J[i1, i0] = mse_line(X, T, (x0[i0], x1[i1]))

plt.figure(figsize=(9.5, 4))
plt.subplots_adjust(wspace=0.5)

ax = plt.subplot(1, 2, 1, projection='3d')
ax.plot_surface(xx0, xx1, J, rstride=10, cstride=10, alpha=0.3, color='blue', edgecolor='black')
ax.set_xticks([-20, 0, 20])
ax.set_yticks([120, 140, 160])
ax.view_init(20, -60)

plt.subplot(1, 2, 2)
cont = plt.contour(xx0, xx1, J, 30, color='black', levels=[100, 1000, 10000, 100000])
cont.clabel(fmt='%1.0d', fontsize=8)
plt.grid(True)
plt.show()
```
### 実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/6ec49004-a041-d9cf-9b64-b2acd7c6a461.png)

w0(横の線), w1(奥の線)の面に対して平均二乗誤差J(高さの線)は谷の形状をしている。
w0の値に対してJが大きく変化しているのは直線の傾きの方がJに及ぼす影響が大きいため。

また、右図の等高線グラフではJの値が100を下回る箇所がある。
これからw0, w1に対してもJは最小値を取ることがわかる。

## 勾配法
Jが最も小さくなるようなw0, w1の求め方として、勾配法がある。
w0, w1面に対して点Jがあり、Jが最も減少する方向へw0, w1を進めることを繰り返し、"お椀状の底"となるw0, w1が求められる。

Jを最小にするには、w0, w1で偏微分する。w0, w1でしたベクトルを(δJ/δw0, δJ/δw1)で表され、Jの勾配と呼び、ΔwJと表す。

w0, w1の値である**w**の更新方法(学習則)を行列表記で表すと、
```
w(t + 1) = w(t) - α * ΔwJ|w(t)
```

t時点のwに対して、学習率αをかけたt時点のJの勾配で差を取ることでt+1時点のwが求めることができる。

平均二乗誤差Jの等高線上に**w**の更新の様子を青い線で図示する。
```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(seed=1)              # 乱数を固定
X_min = 4                           # Xの下限値
X_max = 30                          # Xの上限値
X_n = 16                            # データ個数
X = 5 + 25 * np.random.rand(X_n)
Prm_c = [170, 108, 0.2]             # 生成パラメータ

T = Prm_c[0] - Prm_c[1] * np.exp(-Prm_c[2] * X) + 4 * np.random.randn(X_n)

# 平均誤差関数
def mse_line(x, t, w):
    y = w[0] * x + w[1]
    mse = np.mean((y - t)**2)
    return mse

# 平均二乗誤差の勾配
def dmse_line(x, t, w):
    y = w[0] * x + w[1]
    d_w0 = 2 * np.mean((y - t) * x)
    d_w1 = 2 * np.mean((y - t))
    return d_w0, d_w1

dw_w = dmse_line(X, T, [10, 165])
print(np.round(dw_w))

def fit_line_num(x, t):
    w_init = [10.0, 165.0]  # 初期パラメータ
    alpha = 0.001           # 学習率
    i_max = 100000          # 繰り返しの最大数
    eps = 0.1               # 繰り返しをやめる勾配の絶対値の閾値

    w_i = np.zeros([i_max, 2])
    w_i[0, :] = w_init
    for i in range(1, i_max):
        dmse = dmse_line(x, t, w_i[i - 1])
        w_i[i, 0] = w_i[i - 1, 0] - alpha * dmse[0]
        w_i[i, 1] = w_i[i - 1, 1] - alpha * dmse[1]
        if max(np.absolute(dmse)) < eps:
            break
    w0 = w_i[i, 0]
    w1 = w_i[i, 1]
    w_i = w_i[:i, :]
    return w0, w1, dmse, w_i

plt.figure(figsize=(4, 4))
# MSEの等高線表示
xn = 100    # 等高線表示の解像度
w0_range = [-25, 25]
w1_range = [120, 170]
x0 = np.linspace(w0_range[0], w0_range[1], xn)
x1 = np.linspace(w1_range[0], w1_range[1], xn)
xx0, xx1 = np.meshgrid(x0, x1)

J = np.zeros((len(x0), len(x1)))

for i0 in range(xn):
    for i1 in range(xn):
        J[i1, i0] = mse_line(X, T, (x0[i0], x1[i1]))

cont = plt.contour(xx0, xx1, J, 30, colors='black', levels=[100, 1000, 10000, 100000])
cont.clabel(fmt='%1.0d', fontsize=8)
plt.grid(True)

# 勾配法呼び出し
W0, W1, dMSE, W_history = fit_line_num(X, T)

# 結果表示
print('繰り返し回数 {0}'.format(W_history.shape[0]))
print('W=[{0:.6f}, {1:.6f}'.format(W0, W1))
print('dMSE=[{0:.6f}, {1:.6f}]'.format(dMSE[0], dMSE[1]))
print('MSE={0:.6f}'.format(mse_line(X, T, [W0, W1])))

plt.plot(W_history[:, 0], W_history[:, 1], '.-', color='gray', markersize=10, markeredgecolor='cornflowerblue')
plt.show()
```
### 実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/af6fa667-8063-d5d5-c9f7-cff35666cdf1.png)

```
[5046.  302.]
繰り返し回数 13820
W=[1.539947, 136.176160
dMSE=[-0.005794, 0.099991]
MSE=49.027452
```

ピンの時点から開始し、線の終点で止まることがわかる。
開始地点から徐々にJの小さくなる方向へ向かっていき、終了地点で止まる。

勾配法によって、最小のJがわかったので、データ分布の上に直線を描いてみる。

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(seed=1)              # 乱数を固定
X_min = 4                           # Xの下限値
X_max = 30                          # Xの上限値
X_n = 16                            # データ個数
X = 5 + 25 * np.random.rand(X_n)
Prm_c = [170, 108, 0.2]             # 生成パラメータ

T = Prm_c[0] - Prm_c[1] * np.exp(-Prm_c[2] * X) + 4 * np.random.randn(X_n)

# 平均誤差関数
def mse_line(x, t, w):
    y = w[0] * x + w[1]
    mse = np.mean((y - t)**2)
    return mse

# 平均二乗誤差の勾配
def dmse_line(x, t, w):
    y = w[0] * x + w[1]
    d_w0 = 2 * np.mean((y - t) * x)
    d_w1 = 2 * np.mean((y - t))
    return d_w0, d_w1

# 勾配法
def fit_line_num(x, t):
    w_init = [10.0, 165.0]  # 初期パラメータ
    alpha = 0.001           # 学習率
    i_max = 100000          # 繰り返しの最大数
    eps = 0.1               # 繰り返しをやめる勾配の絶対値の閾値

    w_i = np.zeros([i_max, 2])
    w_i[0, :] = w_init
    for i in range(1, i_max):
        dmse = dmse_line(x, t, w_i[i - 1])
        w_i[i, 0] = w_i[i - 1, 0] - alpha * dmse[0]
        w_i[i, 1] = w_i[i - 1, 1] - alpha * dmse[1]
        if max(np.absolute(dmse)) < eps:
            break
    w0 = w_i[i, 0]
    w1 = w_i[i, 1]
    w_i = w_i[:i, :]
    return w0, w1, dmse, w_i

def show_line(w):
    xb = np.linspace(X_min, X_max, 100)
    y = w[0] * xb + w[1]
    plt.plot(xb, y, color=(.5, .5, .5), linewidth=4)

# 勾配法呼び出し
W0, W1, dMSE, W_history = fit_line_num(X, T)

plt.figure(figsize=(4, 4))
W = np.array([W0, W1])
mse = mse_line(X, T, W)
print("w0={0:.3f}, w1={1:.3f}".format(W0, W1))
print("SD={0:.3f} cm".format(np.sqrt(mse)))
show_line(W)
plt.plot(X, T, marker='o', linestyle='None', color='cornflowerblue', markeredgecolor='black')
plt.xlim(X_min, X_max)
plt.grid(True)
plt.show()
```

### 実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/1ccaf400-3e9c-cc7b-d3c4-1c41da2a68d1.png)
```
w0=1.540, w1=136.176
SD=7.002 cm
```

データと直線は完全に一致しない。
平均二乗誤差の値は49.03cm^2になるが、誤差を二乗しているため、平方根に入れると、7.00cmとなり、「誤差がだいたい7.00cm」になるという。
平均二乗誤差の平方根を**標準偏差**という。
ただし、勾配法で求まる階は極小値であるため、Jが複雑な凹凸がある形をしている場合は近くの凹みで収束してしまうので、一番深い凹み(最小値)を求めることが難しい。

## 直線モデルパラメータの解析解
勾配法のように繰り返しの計算により近似的な値を**数値解**と呼ぶ。
しかし、直線モデルの場合は、近似的な解ではなく、方程式を解くことにより厳密な解を求めることができる。これを**解析解**と呼ぶ。

解析解を使い、データ分布の上に直線を描いてみる。

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(seed=1)              # 乱数を固定
X_min = 4                           # Xの下限値
X_max = 30                          # Xの上限値
X_n = 16                            # データ個数
X = 5 + 25 * np.random.rand(X_n)
Prm_c = [170, 108, 0.2]             # 生成パラメータ

T = Prm_c[0] - Prm_c[1] * np.exp(-Prm_c[2] * X) + 4 * np.random.randn(X_n)

# 平均二乗誤差関数
def mse_line(x, t, w):
    y = w[0] * x + w[1]
    mse = np.mean((y - t)**2)
    return mse

def show_line(w):
    xb = np.linspace(X_min, X_max, 100)
    y = w[0] * xb + w[1]
    plt.plot(xb, y, color=(.5, .5, .5), linewidth=4)

# 解析解
def fit_line(x, t):
    mx = np.mean(x)
    mt = np.mean(t)
    mtx = np.mean(t * x)
    mxx = np.mean(x * x)
    w0 = (mtx - mt * mx) / (mxx - mx**2)
    w1 = mt - w0 * mx
    return np.array([w0, w1])

W = fit_line(X, T)
print("w0={0:.3f}".format(W[0], W[1]))

mse = mse_line(X, T, W)
print("SD={0:.3f}".format(np.sqrt(mse)))

plt.figure(figsize=(4, 4))
show_line(W)
plt.plot(X, T, marker='o', linestyle='None', color='cornflowerblue', markeredgecolor='black')
plt.xlim(X_min, X_max)
plt.grid(True)
plt.show()
```
### 実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/b8957079-4556-3f47-cbdf-b4fcf23d94ff.png)
```
w0=1.558
SD=7.001
```

データ点が直線でフィッティングする場合は解析解が導入できるため、勾配法は使う必要がないということになる。