誤差逆伝搬法②

前回

https://qiita.com/H_T_6221/items/f66eed28045357c6672d

## δE/δwjiを求める
偏微分の連鎖律を使ってδE/δwjiを分解し、入力層から1層への重みパラメータwjiの学習則を導出する。

```math
\frac{δE}{δw_{ji}} = \frac{δE}{δb_{j}} \frac{δb_{j}}{δw_{ji}}
```

初めの中間層ニューロンの入力総和bjでの偏微分δE/δbjを、出力層の入力総和の微分δE/δak=δkの類似性から、以下のように定義する。

```math
\frac{δE}{δb_{j}} = δ_{j}^{(1)}
```

この(1)は1層(中間層)であることを意味する。

偏微分δbj/δwjiは以下のようになる。

```math
\frac{δb_{j}}{δw_{ji}} = \frac{δ}{δw_{ji}} \sum_{i=0}^{D} w_{ji}x_{i} = x_{i}
```

よって、wjiの更新規則は以下のようになる。

```math
w_{ji} (t + 1) = w_{ji} (t) - α \frac{δE}{δw_{ji}} = w_{ji} (t) - α δ_{j}^{(1)} x_{i}
```

これはvkjの更新規則と同じ形をしている。つまり、wjiも結合の先で生じた誤差と、結合のもとの入力に比例する形で変更されるとわかる。

偏微分の連鎖律を使ってδjが何かを求めていく。
```math
δ_{j}^{(1)} = \frac{δE}{δb_{j}} = \{ \sum_{k=0}^{K-1} \frac{δE}{δa_{k}} \frac{δa_{k}}{δz_{j}} \} \frac{δz_{j}}{δb_{j}}
```

分解したδE/δakが2層のδkで表すことができる。δak/δzjは以下のようになる。

```math
\frac{δa_{k}}{δz_{j}} = \frac{δ}{δz_{j}} \sum_{j=0}^{M} v_{kj} z_{j} = v_{kj}
```

δzj/δbjは中間層の活性化関数をh()とすると以下のように表す。

```math
\frac{δz_{j}}{δb_{j}} = \frac{δ}{δb_{j}} h(b_{j}) = h'(b_{j})
```

よって、次の式になる。

```math
δ_{j}^{(1)} = h'(b_{j}) \sum_{k=0}^{K-1} v_{kj} δ_{k}^{(2)}
```

h'(bj)は活性化関数の微分であり、その値は常に正の数になる。次に和の中は出力先の誤差であるδkをvkjの重みで集めている形になっている。つまり、δjは結合先で生じた誤差δkを逆方向に伝達させて計算させているとみなせる。

ここまでの計算は、2層のネットワークを想定していたが、この誤差逆伝搬の法則は、ネットワークの層がもっと増えても出現し、この法則を使えば簡単に重みパラメータの学習則が導き出せる。

この性質から、フィードフォワードニューラルネットワークの勾配法が、誤差逆伝搬法と呼ばれている。

```plantuml
@startuml
left to right direction

rectangle 入力層x {
    usecase x0
    usecase x1
    usecase "x2(ダミー=1)" as x2
}
rectangle 中間層(1層)b/z {
  usecase "b0|z0" as a0
  usecase "b1|z1" as a1
  usecase "z2(ダミー)" as a2
}
rectangle 出力層(2層)a/y {
  usecase "a0|y0" as f1
  usecase "a1|y1" as f2
  usecase "a2|y2" as f3
}

x0 --> a0
x1 --> a0
x2 --> a0

x0 --> a1
x1 --> a1
x2 --> a1

a0 --> f1
a1 --> f1
a2 --> f1

a0 --> f2
a1 --> f2
a2 --> f2

a0 --> f3
a1 --> f3
a2 --> f3
@enduml
```

学習の方法として、以下のようなことをする。

①ネットワークにxを入力し、その出力yを得る。この際に計算されるb,z, aを保持する。
②出力yを目標データtと比べて、その誤差を計算する。この誤差は、出力層の各ニューロンに割り当てる。
③出力層の誤差を使って、中間層の誤差を計算する。(誤差逆伝搬法)
④結合元の信号強度と結合先の誤差の情報を使って、重みパラメータを更新する。

この一連の手続きはデータ1つ分に対して更新する。実際にN個のデータがあるので、データを1つずつ変えながら①～④の手続きをN回行って学習の1ステップとする。

この導出した学習則は2層のフィードフォワードニューラルネットワークに限ったものではなく、3層や4層などのネットワークでも同じ形で誤差を計算して求めることができる。

## 誤差逆伝搬法の実装
誤差逆伝搬法でδE/δw, δE/δv を求めるプログラムを作る。
関数名はdce_fnnとする。入力する方法はce_fnnとし、δE/δw, δE/δvはdw, dvとする。
関数の出力は、これらを結合したdvwとする。

<details>
<summary>実装</summary>

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
import time

# データ生成
np.random.seed(seed=1)
N = 200
K = 3
T = np.zeros((N, 3))
X = np.zeros((N, 2))
X_range0 = [-3, 3]
X_range1 = [-3, 3]
Mu = np.array([[-.5, -.5], [.5, 1.0], [1, -.5]])
Sig = np.array([[.7, .7], [.8, .3], [.3, .8]])
Pi = np.array([0.4, 0.8, 1])

for n in range(N):
    wk = np.random.rand()
    for k in range(K):
        if wk < Pi[k]:
            T[n, k] = 1
            break
    for k in range(2):
        X[n, k] = (np.random.randn() * Sig[T[n, :] == 1, k] + Mu[T[n, :] == 1, k])

# 2分類のデータをテスト・訓練データに分割
TestRatio = 0.5
X_n_training = int(N * TestRatio)
X_train = X[:X_n_training, :]
X_test = X[X_n_training:, :]
T_train = T[:X_n_training, :]
T_test = T[X_n_training:, :]

# データを'class_data.npz'に保存
#np.savez('class_data.npz', X_train=X_train, T_train=T_train, X_test=X_test, T_test=T_test, X_range0=X_range0, X_range1=X_range1)

# データの図示
def show_data(x, t):
    wk, n = t.shape
    c = [[0, 0, 0], [.5, .5, .5], [1, 1, 1]]
    for i in range(n):
        plt.plot(x[t[:, i] == 1, 0], x[t[:, i] == 1, 1], linestyle='none', marker='o', markeredgecolor='black', color=c[i], alpha=0.8)
    plt.grid(True)

# シグモイド関数
def sigmoid(x):
    y = 1 / (1 + np.exp(-x))
    return y

# ネットワーク
def fnn(wv, M, K, x):
    N, D = x.shape
    w = wv[:M * (D + 1)]
    w = w.reshape(M, (D + 1))
    v = wv[M * (D + 1):]
    v = v.reshape((K, M + 1))
    b = np.zeros((N, M + 1))
    z = np.zeros((N, M + 1))
    a = np.zeros((N, K))
    y = np.zeros((N, K))

    for n in range(N):
        for m in range(M):
            b[n, m] = np.dot(w[m, :], np.r_[x[n, :], 1])
            z[n, m] = sigmoid(b[n, m])

        z[n, M] = 1
        wkz = 0
        for k in range(K):
            a[n, k] = np.dot(v[k, :], z[n, :])
            wkz = wkz + np.exp(a[n, k])
        for k in range(K):
            y[n, k] = np.exp(a[n, k]) / wkz
    return y, a, z, b

def ce_fnn(wv, M, K, x, t):
    N, D = x.shape
    y, a, z, b = fnn(wv, M, K, x)
    ce = -np.dot(np.log(y.reshape(-1)), t.reshape(-1)) / N
    return ce

def dce_fnn_num(wv, M, K, x, t):
    epsilon = 0.001
    dwv = np.zeros_like(wv)
    for iwv in range(len(wv)):
        wv_modified = wv.copy()
        wv_modified[iwv] = wv[iwv] - epsilon
        mse1 = ce_fnn(wv_modified, M, K, x, t)
        wv_modified[iwv] = wv[iwv] + epsilon
        mse2 = ce_fnn(wv_modified, M, K, x, t)
        dwv[iwv] = (mse2 - mse1) / (2 * epsilon)
    return dwv

def show_wv(wv, M):
    N = wv.shape[0]
    plt.bar(range(1, M * 3 + 1), wv[:M * 3], align="center", color='black')
    plt.bar(range(M * 3 + 1, N + 1), wv[M * 3:], align="center", color='cornflowerblue')
    plt.xticks(range(1, N + 1))
    plt.xlim(0, N + 1)

def fit_fnn_num(wv_init, M, K, x_train, t_train, x_test, t_test, n, alpha):
    wvt = wv_init
    err_train = np.zeros(n)
    err_test = np.zeros(n)
    wv_hist = np.zeros((n, len(wv_init)))
    for i in range(n):
        wvt = wvt - alpha * dce_fnn_num(wvt, M, K, x_train, t_train)
        err_train[i] = ce_fnn(wvt, M, K, x_train, t_train)
        err_test[i] = ce_fnn(wvt, M, K, x_test, t_test)
        wv_hist[i, :] = wvt
    return wvt, wv_hist, err_train, err_test

# 解析的微分
def dce_fnn(wv, M, K, x, t):
    N, D = x.shape

    # wv をwとvに戻す
    w = wv[:M * (D + 1)]
    w = w.reshape(M, (D + 1))
    v = wv[M * (D + 1):]
    v = v.reshape((K, M + 1))

    # ①xを入力して yを得る
    y, a, z, b = fnn(wv, M, K, x)
    # 出力変数の準備
    dwv = np.zeros_like(wv)
    dw = np.zeros((M, D + 1))
    dv = np.zeros((K, M + 1))
    delta1 = np.zeros(M)
    delta2 = np.zeros(K)

    for n in range(N):
        # ②出力層の誤差を求める
        for k in range(K):
            delta2[k] = (y[n, k] - t[n, k])
        # ③中間層の誤差を求める
        for j in range(M):
            delta1[j] = z[n, j] * (1 - z[n, j]) * np.dot(v[:, j], delta2)
        # ④v の勾配dvを求める
        for k in range(K):
            dv[k, :] = dv[k, :] + delta2[k] + z[n, :] / N
        # ④w の勾配dwを求める
        for j in range(M):
            dw[j, :] = dw[j, :] + delta1[j] + np.r_[x[n, :], 1] / N
    # dw と dv を合体させて dwv とする
    dwv = np.c_[dw.reshape((1, M * (D + 1))), dv.reshape((1, K * (M + 1)))]
    dwv = dwv.reshape(-1)
    return dwv

# Show VW
def show_dwv(wv, M):
    N = wv.shape[0]
    plt.bar(range(1, M * 3 + 1), wv[:M * 3], align="center", color='black')
    plt.bar(range(M * 3 + 1, N + 1), wv[M * 3:], align="center", color='cornflowerblue')
    plt.xticks(range(1, N + 1))
    plt.xlim(0, N + 1)

# 動作確認
M = 2
K = 3
N = 2
nWM = M * 3 + K * (M + 1)
np.random.seed(1)
WV = np.random.normal(0, 1, nWM)

dWV_ana = dce_fnn(WV, M, K, X_train[:N, :], T_train[:N, :])
print("analytical dWV")
print(dWV_ana)

dWV_num = dce_fnn_num(WV, M, K, X_train[:N, :], T_train[:N, :])
print("numerical dWV")
print(dWV_num)

plt.figure(1, figsize=(8, 3))
plt.subplots_adjust(wspace=0.5)
plt.subplot(1, 2, 1)
show_dwv(dWV_ana, M)
plt.title('analitical')
plt.subplot(1, 2, 2)
show_dwv(dWV_num, M)
plt.title('numerical')
plt.show()
```

</details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/7caf3a9d-1136-db13-0640-bbc099495b15.png)

左図は解析的微分、右図は数値微分になる。
ランダムに生成した重みパラメータWVに対して解析的微分値dWV_anaを出力し、前回に作成した数値微分値sWV_numの値も表示する。
解析的微分の値は、以前に計算した数値微分の値とほぼ一致することがわかる。

数値微分で解いた分類問題を誤差逆伝搬法で解いてみる。

<details>
<summary>実装</summary>

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
import time

# データ生成
np.random.seed(seed=1)
N = 200
K = 3
T = np.zeros((N, 3))
X = np.zeros((N, 2))
X_range0 = [-3, 3]
X_range1 = [-3, 3]
Mu = np.array([[-.5, -.5], [.5, 1.0], [1, -.5]])
Sig = np.array([[.7, .7], [.8, .3], [.3, .8]])
Pi = np.array([0.4, 0.8, 1])

for n in range(N):
    wk = np.random.rand()
    for k in range(K):
        if wk < Pi[k]:
            T[n, k] = 1
            break
    for k in range(2):
        X[n, k] = (np.random.randn() * Sig[T[n, :] == 1, k] + Mu[T[n, :] == 1, k])

# 2分類のデータをテスト・訓練データに分割
TestRatio = 0.5
X_n_training = int(N * TestRatio)
X_train = X[:X_n_training, :]
X_test = X[X_n_training:, :]
T_train = T[:X_n_training, :]
T_test = T[X_n_training:, :]

# データを'class_data.npz'に保存
#np.savez('class_data.npz', X_train=X_train, T_train=T_train, X_test=X_test, T_test=T_test, X_range0=X_range0, X_range1=X_range1)

# データの図示
def show_data(x, t):
    wk, n = t.shape
    c = [[0, 0, 0], [.5, .5, .5], [1, 1, 1]]
    for i in range(n):
        plt.plot(x[t[:, i] == 1, 0], x[t[:, i] == 1, 1], linestyle='none', marker='o', markeredgecolor='black', color=c[i], alpha=0.8)
    plt.grid(True)

# シグモイド関数
def sigmoid(x):
    y = 1 / (1 + np.exp(-x))
    return y

# ネットワーク
def fnn(wv, M, K, x):
    N, D = x.shape
    w = wv[:M * (D + 1)]
    w = w.reshape(M, (D + 1))
    v = wv[M * (D + 1):]
    v = v.reshape((K, M + 1))
    b = np.zeros((N, M + 1))
    z = np.zeros((N, M + 1))
    a = np.zeros((N, K))
    y = np.zeros((N, K))

    for n in range(N):
        for m in range(M):
            b[n, m] = np.dot(w[m, :], np.r_[x[n, :], 1])
            z[n, m] = sigmoid(b[n, m])

        z[n, M] = 1
        wkz = 0
        for k in range(K):
            a[n, k] = np.dot(v[k, :], z[n, :])
            wkz = wkz + np.exp(a[n, k])
        for k in range(K):
            y[n, k] = np.exp(a[n, k]) / wkz
    return y, a, z, b

def ce_fnn(wv, M, K, x, t):
    N, D = x.shape
    y, a, z, b = fnn(wv, M, K, x)
    ce = -np.dot(np.log(y.reshape(-1)), t.reshape(-1)) / N
    return ce

def dce_fnn_num(wv, M, K, x, t):
    epsilon = 0.001
    dwv = np.zeros_like(wv)
    for iwv in range(len(wv)):
        wv_modified = wv.copy()
        wv_modified[iwv] = wv[iwv] - epsilon
        mse1 = ce_fnn(wv_modified, M, K, x, t)
        wv_modified[iwv] = wv[iwv] + epsilon
        mse2 = ce_fnn(wv_modified, M, K, x, t)
        dwv[iwv] = (mse2 - mse1) / (2 * epsilon)
    return dwv

def show_wv(wv, M):
    N = wv.shape[0]
    plt.bar(range(1, M * 3 + 1), wv[:M * 3], align="center", color='black')
    plt.bar(range(M * 3 + 1, N + 1), wv[M * 3:], align="center", color='cornflowerblue')
    plt.xticks(range(1, N + 1))
    plt.xlim(0, N + 1)

def fit_fnn_num(wv_init, M, K, x_train, t_train, x_test, t_test, n, alpha):
    wvt = wv_init
    err_train = np.zeros(n)
    err_test = np.zeros(n)
    wv_hist = np.zeros((n, len(wv_init)))
    for i in range(n):
        wvt = wvt - alpha * dce_fnn_num(wvt, M, K, x_train, t_train)
        err_train[i] = ce_fnn(wvt, M, K, x_train, t_train)
        err_test[i] = ce_fnn(wvt, M, K, x_test, t_test)
        wv_hist[i, :] = wvt
    return wvt, wv_hist, err_train, err_test

def show_fnn(wv, M, K):
    xn = 60
    x0 = np.linspace(X_range0[0], X_range0[1], xn)
    x1 = np.linspace(X_range1[0], X_range1[1], xn)
    xx0, xx1 = np.meshgrid(x0, x1)
    x = np.c_[np.reshape(xx0, xn * xn), np.reshape(xx1, xn * xn)]
    y, a, z, b = fnn(wv, M, K, x)
    plt.figure(1, figsize=(4, 4))
    for ic in range(K):
        f = y[:, ic]
        f = f.reshape(xn, xn)
        f = f.T
        cont = plt.contour(xx0, xx1, f, levels=[0.8, 0.9], colors=['cornflowerblue', 'black'])
        cont.clabel(fmt='%1.1f', fontsize=9)
    plt.xlim(X_range0)
    plt.ylim(X_range1)

# 解析的微分
def dce_fnn(wv, M, K, x, t):
    N, D = x.shape

    # wv をwとvに戻す
    w = wv[:M * (D + 1)]
    w = w.reshape(M, (D + 1))
    v = wv[M * (D + 1):]
    v = v.reshape((K, M + 1))

    # ①xを入力して yを得る
    y, a, z, b = fnn(wv, M, K, x)
    # 出力変数の準備
    dwv = np.zeros_like(wv)
    dw = np.zeros((M, D + 1))
    dv = np.zeros((K, M + 1))
    delta1 = np.zeros(M)
    delta2 = np.zeros(K)

    for n in range(N):
        # ②出力層の誤差を求める
        for k in range(K):
            delta2[k] = (y[n, k] - t[n, k])
        # ③中間層の誤差を求める
        for j in range(M):
            delta1[j] = z[n, j] * (1 - z[n, j]) * np.dot(v[:, j], delta2)
        # ④v の勾配dvを求める
        for k in range(K):
            dv[k, :] = dv[k, :] + delta2[k] + z[n, :] / N
        # ④w の勾配dwを求める
        for j in range(M):
            dw[j, :] = dw[j, :] + delta1[j] + np.r_[x[n, :], 1] / N
    # dw と dv を合体させて dwv とする
    dwv = np.c_[dw.reshape((1, M * (D + 1))), dv.reshape((1, K * (M + 1)))]
    dwv = dwv.reshape(-1)
    return dwv

# Show VW
def show_dwv(wv, M):
    N = wv.shape[0]
    plt.bar(range(1, M * 3 + 1), wv[:M * 3], align="center", color='black')
    plt.bar(range(M * 3 + 1, N + 1), wv[M * 3:], align="center", color='cornflowerblue')
    plt.xticks(range(1, N + 1))
    plt.xlim(0, N + 1)

def fit_fnn(wv_init, M, K, x_train, t_train, x_test, t_test, n, alpha):
    wv = wv_init.copy()
    err_train = np.zeros(n)
    err_test = np.zeros(n)
    wv_hist = np.zeros((n, len(wv_init)))
    for i in range(n):
        wv = wv - alpha * dce_fnn_num(wv, M, K, x_train, t_train)
        err_train[i] = ce_fnn(wv, M, K, x_train, t_train)
        err_test[i] = ce_fnn(wv, M, K, x_test, t_test)
        wv_hist[i, :] = wv
    return wv, wv_hist, err_train, err_test

# 動作確認
startTime = time.time()

M = 2
K = 3
np.random.seed(1)
WV_init = np.random.normal(0, 0.01, M * 3 + K * (M + 1))

N_step = 1000
alpha = 1

WV, WV_hist, Err_train, Err_test = fit_fnn(WV_init, M, K, X_train, T_train, X_test, T_test, N_step, alpha)

calculation_time = time.time() - startTime
print("Calculation time:{0:.3f} sec".format(calculation_time))

plt.figure(1, figsize=(12, 3))
plt.subplots_adjust(wspace=0.5)

# 学習誤差の表示
plt.subplot(1, 3, 1)
plt.plot(Err_train, 'black', label='training')
plt.plot(Err_test, 'cornflowerblue', label='test')
plt.legend()

# 重みの時間発展の表示
plt.subplot(1, 3, 2)
plt.plot(WV_hist[:, :M * 3], 'black')
plt.plot(WV_hist[:, M * 3:], 'cornflowerblue')

# 境界線の表示
plt.subplot(1, 3, 3)
show_data(X_test, T_test)
M = 2
K = 3
show_fnn(WV, M, K)
plt.show()
```

</details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/ad28d1e7-5620-35f7-a517-07c250a558a6.png)

## 学習後のニューロンの特性
何故、2層のネットワークは曲線の境界線を作り出すことができたのか。学習の結果、各ニューロンはどういう性質を獲得したのか。bj, zj, ak, ykの特性を図示する。

<details>
<summary>実装</summary>

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
import time

# データ生成
np.random.seed(seed=1)
N = 200
K = 3
T = np.zeros((N, 3))
X = np.zeros((N, 2))
X_range0 = [-3, 3]
X_range1 = [-3, 3]
Mu = np.array([[-.5, -.5], [.5, 1.0], [1, -.5]])
Sig = np.array([[.7, .7], [.8, .3], [.3, .8]])
Pi = np.array([0.4, 0.8, 1])

for n in range(N):
    wk = np.random.rand()
    for k in range(K):
        if wk < Pi[k]:
            T[n, k] = 1
            break
    for k in range(2):
        X[n, k] = (np.random.randn() * Sig[T[n, :] == 1, k] + Mu[T[n, :] == 1, k])

# 2分類のデータをテスト・訓練データに分割
TestRatio = 0.5
X_n_training = int(N * TestRatio)
X_train = X[:X_n_training, :]
X_test = X[X_n_training:, :]
T_train = T[:X_n_training, :]
T_test = T[X_n_training:, :]

# データを'class_data.npz'に保存
#np.savez('class_data.npz', X_train=X_train, T_train=T_train, X_test=X_test, T_test=T_test, X_range0=X_range0, X_range1=X_range1)

# データの図示
def show_data(x, t):
    wk, n = t.shape
    c = [[0, 0, 0], [.5, .5, .5], [1, 1, 1]]
    for i in range(n):
        plt.plot(x[t[:, i] == 1, 0], x[t[:, i] == 1, 1], linestyle='none', marker='o', markeredgecolor='black', color=c[i], alpha=0.8)
    plt.grid(True)

# シグモイド関数
def sigmoid(x):
    y = 1 / (1 + np.exp(-x))
    return y

# ネットワーク
def fnn(wv, M, K, x):
    N, D = x.shape
    w = wv[:M * (D + 1)]
    w = w.reshape(M, (D + 1))
    v = wv[M * (D + 1):]
    v = v.reshape((K, M + 1))
    b = np.zeros((N, M + 1))
    z = np.zeros((N, M + 1))
    a = np.zeros((N, K))
    y = np.zeros((N, K))

    for n in range(N):
        for m in range(M):
            b[n, m] = np.dot(w[m, :], np.r_[x[n, :], 1])
            z[n, m] = sigmoid(b[n, m])

        z[n, M] = 1
        wkz = 0
        for k in range(K):
            a[n, k] = np.dot(v[k, :], z[n, :])
            wkz = wkz + np.exp(a[n, k])
        for k in range(K):
            y[n, k] = np.exp(a[n, k]) / wkz
    return y, a, z, b

def ce_fnn(wv, M, K, x, t):
    N, D = x.shape
    y, a, z, b = fnn(wv, M, K, x)
    ce = -np.dot(np.log(y.reshape(-1)), t.reshape(-1)) / N
    return ce

def dce_fnn_num(wv, M, K, x, t):
    epsilon = 0.001
    dwv = np.zeros_like(wv)
    for iwv in range(len(wv)):
        wv_modified = wv.copy()
        wv_modified[iwv] = wv[iwv] - epsilon
        mse1 = ce_fnn(wv_modified, M, K, x, t)
        wv_modified[iwv] = wv[iwv] + epsilon
        mse2 = ce_fnn(wv_modified, M, K, x, t)
        dwv[iwv] = (mse2 - mse1) / (2 * epsilon)
    return dwv

def show_wv(wv, M):
    N = wv.shape[0]
    plt.bar(range(1, M * 3 + 1), wv[:M * 3], align="center", color='black')
    plt.bar(range(M * 3 + 1, N + 1), wv[M * 3:], align="center", color='cornflowerblue')
    plt.xticks(range(1, N + 1))
    plt.xlim(0, N + 1)

def fit_fnn_num(wv_init, M, K, x_train, t_train, x_test, t_test, n, alpha):
    wvt = wv_init
    err_train = np.zeros(n)
    err_test = np.zeros(n)
    wv_hist = np.zeros((n, len(wv_init)))
    for i in range(n):
        wvt = wvt - alpha * dce_fnn_num(wvt, M, K, x_train, t_train)
        err_train[i] = ce_fnn(wvt, M, K, x_train, t_train)
        err_test[i] = ce_fnn(wvt, M, K, x_test, t_test)
        wv_hist[i, :] = wvt
    return wvt, wv_hist, err_train, err_test

def show_fnn(wv, M, K):
    xn = 60
    x0 = np.linspace(X_range0[0], X_range0[1], xn)
    x1 = np.linspace(X_range1[0], X_range1[1], xn)
    xx0, xx1 = np.meshgrid(x0, x1)
    x = np.c_[np.reshape(xx0, xn * xn), np.reshape(xx1, xn * xn)]
    y, a, z, b = fnn(wv, M, K, x)
    plt.figure(1, figsize=(4, 4))
    for ic in range(K):
        f = y[:, ic]
        f = f.reshape(xn, xn)
        f = f.T
        cont = plt.contour(xx0, xx1, f, levels=[0.8, 0.9], colors=['cornflowerblue', 'black'])
        cont.clabel(fmt='%1.1f', fontsize=9)
    plt.xlim(X_range0)
    plt.ylim(X_range1)

# 解析的微分
def dce_fnn(wv, M, K, x, t):
    N, D = x.shape

    # wv をwとvに戻す
    w = wv[:M * (D + 1)]
    w = w.reshape(M, (D + 1))
    v = wv[M * (D + 1):]
    v = v.reshape((K, M + 1))

    # ①xを入力して yを得る
    y, a, z, b = fnn(wv, M, K, x)
    # 出力変数の準備
    dwv = np.zeros_like(wv)
    dw = np.zeros((M, D + 1))
    dv = np.zeros((K, M + 1))
    delta1 = np.zeros(M)
    delta2 = np.zeros(K)

    for n in range(N):
        # ②出力層の誤差を求める
        for k in range(K):
            delta2[k] = (y[n, k] - t[n, k])
        # ③中間層の誤差を求める
        for j in range(M):
            delta1[j] = z[n, j] * (1 - z[n, j]) * np.dot(v[:, j], delta2)
        # ④v の勾配dvを求める
        for k in range(K):
            dv[k, :] = dv[k, :] + delta2[k] + z[n, :] / N
        # ④w の勾配dwを求める
        for j in range(M):
            dw[j, :] = dw[j, :] + delta1[j] + np.r_[x[n, :], 1] / N
    # dw と dv を合体させて dwv とする
    dwv = np.c_[dw.reshape((1, M * (D + 1))), dv.reshape((1, K * (M + 1)))]
    dwv = dwv.reshape(-1)
    return dwv

# Show VW
def show_dwv(wv, M):
    N = wv.shape[0]
    plt.bar(range(1, M * 3 + 1), wv[:M * 3], align="center", color='black')
    plt.bar(range(M * 3 + 1, N + 1), wv[M * 3:], align="center", color='cornflowerblue')
    plt.xticks(range(1, N + 1))
    plt.xlim(0, N + 1)

def fit_fnn(wv_init, M, K, x_train, t_train, x_test, t_test, n, alpha):
    wv = wv_init.copy()
    err_train = np.zeros(n)
    err_test = np.zeros(n)
    wv_hist = np.zeros((n, len(wv_init)))
    for i in range(n):
        wv = wv - alpha * dce_fnn_num(wv, M, K, x_train, t_train)
        err_train[i] = ce_fnn(wv, M, K, x_train, t_train)
        err_test[i] = ce_fnn(wv, M, K, x_test, t_test)
        wv_hist[i, :] = wv
    return wv, wv_hist, err_train, err_test

def show_activation3d(ax, v, v_ticks, title_str):
    f = v.copy()
    f = f.reshape(xn, xn)
    f = f.T
    ax.plot_surface(xx0, xx1, f, color='blue', edgecolor='black', rstride=1, cstride=1, alpha=0.5)
    ax.view_init(70, -110)
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    ax.set_title(title_str, fontsize=18)


# 動作確認
startTime = time.time()

M = 2
K = 3
np.random.seed(1)
WV_init = np.random.normal(0, 0.01, M * 3 + K * (M + 1))

N_step = 1000
alpha = 1

WV, WV_hist, Err_train, Err_test = fit_fnn(WV_init, M, K, X_train, T_train, X_test, T_test, N_step, alpha)

calculation_time = time.time() - startTime
print("Calculation time:{0:.3f} sec".format(calculation_time))

xn = 15
x0 = np.linspace(X_range0[0], X_range0[1], xn)
x1 = np.linspace(X_range1[0], X_range1[1], xn)
xx0, xx1 = np.meshgrid(x0, x1)
x = np.c_[np.reshape(xx0, xn * xn), np.reshape(xx1, xn * xn)]
y, a, z, b = fnn(WV, M, K, x)

fig = plt.figure(1, figsize=(12, 9))
plt.subplots_adjust(left=0.075, bottom=0.05, right=0.95, top=0.95, wspace=0.4, hspace=0.4)

for m in range(M):
    ax = fig.add_subplot(3, 4, 1 + m * 4, projection='3d')
    show_activation3d(ax, b[:, m], [-10, 10], '$b_{0:d}$'.format(m))
    ax = fig.add_subplot(3, 4, 2 + m * 4, projection='3d')
    show_activation3d(ax, y[:, k], [0, 1], '$z_{0:d}$'.format(m))

for k in range(K):
    ax = fig.add_subplot(3, 4, 3 + k * 4, projection='3d')
    show_activation3d(ax, a[:, k], [-5, 5], '$a_{0:d}$'.format(k))
    ax = fig.add_subplot(3, 4, 4 + k * 4, projection='3d')
    show_activation3d(ax, y[:, k], [0, 1], '$y_{0:d}$'.format(k))

plt.show()
```

</details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3850581/9056a434-f61f-2aa1-8ca9-3846e42373ef.png)

入力総和bjの入出力マップは、シグモイド関数δ()をくぐると、低い部分と高い部分がそれぞれ0と1の範囲に押し込められて、出力zjとなる。
出力層の入力総和akの入出力マップは、z0, z1の2つの入出力マップの線形和で作られる。
aはソフトマックス関数を通り、0～1の範囲に押しつぶされ、ykが作られる。y0, y1, y2の盛り上がっている部分は、それぞれクラスに分類される範囲に対応している。
ykは、ソフトマックス関数で加工された結果なので、各面を全て足せば、高さが1の平面になる。

出力層は、中間層のニューロンを基底関数を線形ソフトマックスモデルとみなすことができる。しかし、ニューラルネットワークモデルの場合は、基底関数の特性も学習によって自動で最適化される特別な線形ソフトマックスモデルなのである。