前回に引き続いて、ニューラルネットワークの演習問題をPythonで実装します。今回は後方誤差伝播（Back Propagation）も行うため、本格的なニューラルネットワークです。

組み込みで実装する場合は一瞬ですが、全部自力で実装する場合は本当に骨が折れます。覚悟してください。自分はもう二度と自力でやりたくありません。

追記：gistからソースファイルをダウンロードできます
自分で実装：https://gist.github.com/koshian2/e75c4120764ce25a5b2c0ecf6626e490
組み込み：https://gist.github.com/koshian2/4b26f582fa3bbcc0ab28e119e3cbda8b

これまでの目次

##◎sklearn.neural_networkの組み込みでニューラルネットワーク
実は前回でやったのとほとんど変わりません。組み込みの場合はコーディングがすぐ終わるので楽しいです。正則化のパラメーターαを変更して、どのように変わるかを観察します。

自分で実装する場合と同じの、隠れ層が1つだけでそこには26ノードがあるニューラルネットワークを作ります。

import numpy as np
import matplotlib.pyplot as plt

from scipy.io import loadmat
from sklearn.neural_network import MLPClassifier

# データの読み込み
def load_data1():
    data = loadmat("ex4data1")
    return np.array(data['X']), np.ravel(np.array(data['y']))

X_data, y = load_data1()
m = len(X_data[:, 1])

## ネットワークの可視化
def display_data(X, alpha, accuracy):
    fig = plt.figure(figsize = (5, 5))
    fig.subplots_adjust(hspace=0.05, wspace=0.05)
    for i in range(X.shape[0]):
        ax = fig.add_subplot(5, 5, i+1, xticks=[], yticks=[])
        ax.imshow(X[i, :].reshape(20, 20, order="F"), cmap="gray")
    plt.suptitle(f"Alpha = {alpha}, Accuracy = {accuracy}")
    plt.show()

# 正則化の値を変えてみる
alphas = [0.0001, 1, 50, 1e10]
for a in alphas:
    print("Alpha =", a)
    clf = MLPClassifier(hidden_layer_sizes=(26, ), solver="adam", random_state=114514, max_iter=500, alpha=a)
    clf.fit(X_data, y)
    accuracy = clf.score(X_data, y)*100
    print("Training Set Accuracy:" , accuracy)
    print()
    display_data(clf.coefs_[0][:, 1:].T, a, np.round(accuracy, 3))

出力

Alpha = 0.0001
Training Set Accuracy: 100.0

Alpha = 1
Training Set Accuracy: 95.7

Alpha = 50
Training Set Accuracy: 53.24

Alpha = 10000000000.0
Training Set Accuracy: 10.0

長くても20秒ぐらいで出てきます。αを小さくするほどオーバーフィッティングしやすく、大きくするほどアンダーフィッティングしやすくなります。αのデフォルトは0.0001です（1番目）。scikit-learn組み込みのLogisticRegressionの場合、λの逆数であるCを使っていましたが、このαはλと同じように見えます（逆数でないことは確かですが、スケールが完全に同一のものかはわかりません）。

本来のオーバー／アンダーフィッティングの判定は、訓練データではなくテストデータに対して行うべきなので、1番目の100％が本当に完璧なモデルができたのか、オーバーフィッティングなのかはこれだけでは判断できません。

入力層→隠れ層へのパラメーター（バイアスノードの1つは抜いています）を可視化すると次のようになります。

詳細は後ほど説明しますが、白いところは重要視している点（ピクセル）、黒い所はウェイトを軽く置いているところです。αが低いのでかなりべちょっとしていますね。

αを1にしました。こちらのほうが数字の輪郭のポイントを抑えているような感じはします。

αを50にしました。かなりぼやけてきています。どれも0か1のように見えます。

αを限りなく大きくしてみました。精度10%なので、サイコロ振って当てたのと同じ精度です。ほぼ真っ黒か真っ白で、ノードが事実上何もしていないことがわかります。
こうやって可視化してみると、目が見えるプロセスを追体験しているようで、「ニューラルネットワーク＝神経回路」と言われても確かにそうかもね、と思えるようになります。

##・ニューラルネットワークを自分で実装
1からコードを書いていくマゾはこちら。一度やれば十分ですし、一度やって直感的にわかっても、何やってるんだがはっきりは腑に落ちなかったです。長いので覚悟してください。

###1.データの読み込み

import numpy as np
import copy
import matplotlib.pyplot as plt

from scipy.io import loadmat

## データの読み込み
def load_data1():
    data = loadmat("ex4data1")
    return np.array(data['X']), np.ravel(np.array(data['y']))

X_data, y = load_data1()
m = len(X_data[:, 1])

# Theta1,2の仮の値を読み込む
def load_weights():
    data = loadmat("ex4weights")
    return np.array(data['Theta1']), np.array(data['Theta2'])

Theta1, Theta2 = load_weights()

演習で使うデータが7MB近くあるので、ハードコーディングは諦めました。MNISTのデータを5000個集めただけなので、scikit-learnに同梱のMNISTのデータを使っても同じことができます。

ex4weightsは前回使った、既に計算済みの$\Theta^{(1)},\Theta^{(2)}$です。最終的にこの値は最適化するので、本来はいらないものです。コスト関数がちゃんと実装できているかの確認に使います。

###2.定数と活性化関数の定義

# 定数の定義
INPUT_LAYER_SIZE = 400
HIDDEN_LAYER_SIZE = 25
NUM_LABELS = 10

## 関数の定義
# シグモイド関数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

# シグモイド関数の導関数
def sigmoid_gradient(z):
    g = sigmoid(z)
    return g * (1 - g)

ここは特に問題ないです。シグモイド関数の微分について解析的に証明しておきましょう。シグモイド関数の定義は

$$g(x) = \frac{1}{1+e^{-x}}$$

でした。これを微分します。

\begin{align}
\frac{d}{dx}g(x)&=-\frac{1}{(1+e^{-x})^2}\frac{d}{dx}(1+e^{-x}) \\
&=\frac{1}{1+e^{-x}}\frac{e^{-x}}{1+e^{-x}} \\
&=\frac{1}{1+e^{-x}}\frac{1+e^{-x}-1}{1+e^{-x}} \\
&=\frac{1}{1+e^{-x}}\bigl(1-\frac{1}{1+e^{-x}}\bigr) \\
&=g(x)(1-g(x))
\end{align}

微分してもきれいな積で表されます。

###3.コスト関数の定義
ここが鬼門です。コスト関数の中でForward PropagationとBack Propagationの両方をやらないといけないので分量が多いです。Andrew Ng先生の講義を見たほうが数百倍わかりやすいと思いますが、Forward PropagationとBack Propagationがそれぞれ何をやっているか概念的に説明を。

Forward Propagationは、yについての推定値の$h_\Theta(x)$とコスト関数$J(\Theta)$の値を計算するためのものです。これはロジスティック回帰の分類器を多層かつ並列において、再帰的にコストを計算している（NNはそういうもの）だけです。
最急降下法でコストを最小化する際に、コスト関数の勾配の値が必要になることはこれまでと変わりません。通常の回帰分析やロジスティック回帰では、$J(\Theta)$の形が与えられれば数式展開で勾配が計算できましたが、ニューラルネットワークではコスト関数が非線形になってしまうので、解析的に解くことがまず不可能です。

そこで勾配を計算するために、Back Propagationというアルゴリズムを使います。ネットを見ているといろんな説明がありますが、Courseraの講座でAndrew Ng先生が言っていた「Back Propagationは勾配を計算するために効率の良いアルゴリズム（要約）」という説明が一番しっくりきました。後の勾配チェック（Gradient Checking）でも触れますが、コンピューターで解かせる場合、$J(\Theta)$に近い2つの値を入れて微分係数を近似的に求める方法よりも遥かに計算量は少ないです。
Back Propagationは文字通りネットワークを逆にたどって計算していくことなのですが、「そんなんで勾配が求まるの？」と思うかもしれません。自分も不思議だったのですが、Gradient Checkingをすると確かにBack Propagationで勾配が求まっています。隠れ層なし2ノードで$ax+b$みたいな簡単なロジスティック回帰をNN上でやると、Back Propagationでの値と解析的に解いた勾配の値が同じ――といったことが確認できるはずです。自分は力尽きました。

前置きが長くなってしまいましたが、とりあえずコードを。

# コスト関数
def nn_cost_function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda_):
    # UnrollではなくTupleとして入れる
    Theta1 = nn_params[0]
    Theta2 = nn_params[1]
    m = X.shape[0]
    # yが1～10の値なので分類器ごとのブール値になるように変換
    Y = np.zeros((m, num_labels))
    for i in range(num_labels):
        Y[:, i] = y == (i+1)

    ## Forward Propagation
    # 入力層にバイアスユニットを追加
    A1 = np.c_[np.ones((m, 1)), X]
    # 第2層のパラメーターを計算
    Z2 = np.dot(A1, Theta1.T)
    A2 = sigmoid(Z2)
    # 第2層にバイアスユニットを追加
    A2 = np.c_[np.ones((A2.shape[0], 1)), A2]
    # 推定値
    h_theta = sigmoid(np.dot(A2, Theta2.T))
    # 誤差の計算
    J = np.sum(np.sum(-Y * np.log(h_theta) - (1 - Y) * np.log(1 - h_theta))) / m
    J += lambda_ / 2 / m * (np.sum(np.sum(Theta1[:, 1:] ** 2)) + np.sum(np.sum(Theta2[:, 1:] ** 2)))

    ## Back Propagation
    # δ^(3)＝推定値-訓練データのY
    delta_3 = h_theta - Y
    # δ^(2)の計算
    delta_2 = np.dot(delta_3, Theta2) * np.c_[np.zeros((m, 1)), sigmoid_gradient(Z2)]
    # δ^(2)からバイアスユニットのパラメーターを取り除く
    delta_2 = delta_2[:, 1:]
    # 誤差を集計し、⊿^(1), ⊿^(2)を計算
    Delta2 = np.dot(delta_3.T, A2)
    Delta1 = np.dot(delta_2.T, A1)
    # 正則化を除いた勾配行列を求める
    Theta2_grad = Delta2 / m
    Theta1_grad = Delta1 / m

    ## 正則化の項を勾配に追加
    if lambda_ > 0:
        # Numpyの行列の代入は参照渡しなので、ディープコピーする
        temp = copy.deepcopy(Theta1)
        temp[:, 0] = 0
        Theta1_grad += temp * (lambda_ / m)
        # Theta2についても同様に足す
        temp = copy.deepcopy(Theta2)
        temp[:, 0] = 0
        Theta2_grad += temp * (lambda_ / m)

    # 返り値はJ, (Theta1_grad, Theta2_grad)
    return J, (Theta1_grad, Theta2_grad)

Octaveでは行列をベクトルに変換して（Unroll）使っていましたが、PythonではTupleなどのデータ構造があるので、わざわざ展開する必要はないかと思います。Pythonでも最適化する際にUnrollしたほうが便利ですが、行列に再代入するときにバグが怖いので好みが分かれると思います。行列ではなくテンソルで扱うと明瞭になるかもしれない。

yの入力値は1～10までの値ですが、出力層は$y==i$を推定するロジスティック回帰が10個あることを思い出すと、yはこれらのフラグの行列になります。

Forward PropagationとBack Propagationコメントを追っていけば多分分かるはずです（書いててもだいたいこんな感じだろうで、さっぱり腑に落ちねえ！）。

正則化はこれまでと同様です。Pythonの参照渡しに注意しましょう。置き換えて代入し、元の値を残しておきたいときはディープコピーしておきます。

###4.初期値のランダム化

# ランダムな初期値を設定
def rand_initialize_weights(L_in, L_out):
    # Θに係数に対して、（Forward Propagationで）入力元のレイヤーのユニットの数=L_in
    # 出ていくレイヤーのユニットの数=L_out
    epsilon_init = 0.12
    # ε_initをアドホックに与えたが、ε_init = √6÷√(L_in + L_out)とおくとよいらしい
    # ε_init=0.12 で逆算すると、L_in + L_outが416.7になるので、入力層=401→隠れ層=25なのでだいたいあってる
    return np.random.uniform(-epsilon_init, epsilon_init, (L_out, 1+L_in))

回帰分析やロジスティック回帰では問題になりませんでしたが、ニューラルネットワーク特有の問題として、Θの初期値を全て同じ値（例えば全部0）にすると、ノードのウェイトが全て同じになってしまい、初期条件の対称性（Symmetry）が反復しても解消されず、良い値に収束しないというものがあります。初期値を全て同じ0と決め打ちすると、Forward Propagationが同一層（この式では2層目だが他の層でも同じ）で全て等しくなる、つまり、

$$a_1^{(2)}=a_2^{(2)}=\dots$$

同様にBack Propagationを実行していくと、同様に等しくなり、

$$\delta_1^{(2)}=\delta_2^{(2)}=\dots$$

コスト関数の偏微分も同じく等しくなってしまいます。

$$\frac{\partial}{\partial\Theta_{01}^{(1)}}J(\Theta)=\frac{\partial}{\partial\Theta_{02}^{(1)}}J(\Theta)=\dots $$

つまり最急降下法を導入すると、Θが直線的にしか変化しません。最急降下法の解が初期値に依存することは、局所解が複数ある非凸関数の最適化ならどれでもおこること（NNは非凸関数）なので、たまたま対称性というワードが出てきたものの、これ自体は特に不思議なことではありません。

この対称性を破るには、$\epsilon_{init}$という小さな値を与えて、$-\epsilon_{init}$から$\epsilon_{init}$の間の一様乱数を初期値を入れるといいそうです。このコードではアドホックに0.12と入れていますが、具体的に良い$\epsilon_{init}$の値は、

$$\epsilon_{init}=\frac{\sqrt6}{\sqrt{s_l+s_{l+1}}}$$

で計算できるとのことです。ここで$s_l$と$s_{l+1}$は$\Theta^{(l)}$が対応する層のノードの数を表します。

###5.コスト関数のテスト（１）
コスト関数が正しく計算できているかテストします。まずはコストの値$J(\Theta)$から。正則化なし（λ＝0）、正則化あり（λ＝1）の両方で確認します。

## 関数のテスト
# 事前に読み込んだTheta1, Theta2で計算, λ=0とする（正則化なし）
lambda_ = 0
J,grad = nn_cost_function((Theta1, Theta2), INPUT_LAYER_SIZE, HIDDEN_LAYER_SIZE, NUM_LABELS, X_data, y, lambda_)
print("Feedforward Using Neural Network ...")
print("Cost at parameters (loaded from ex4weights):", J)
print("(this value should be about 0.287629)\n")
# λ＝1で正則化
lambda_ = 1
J,grad = nn_cost_function((Theta1, Theta2), INPUT_LAYER_SIZE, HIDDEN_LAYER_SIZE, NUM_LABELS, X_data, y, lambda_)
print("Checking Cost Function (w/ Regularization) ...")
print("Cost at parameters (loaded from ex4weights):", J)
print("(this value should be about 0.383770)\n")

大丈夫そうですね。

出力

Cost at parameters (loaded from ex4weights): 0.287629165161
(this value should be about 0.287629)

Checking Cost Function (w/ Regularization) ...
Cost at parameters (loaded from ex4weights): 0.383769859091
(this value should be about 0.383770)

###6.シグモイド関数の微分の値を確認
ちょっと脱線して、シグモイド関数の導関数を定義します。

g = sigmoid_gradient(np.array([-1, -0.5, 0, 0.5, 1]))
print("Sigmoid gradient evaluated at [-1 -0.5 0 0.5 1]:")
print(g)
print()

まあこういう関数ですよと。0を中心に山のできるグラフになります（ガウス関数と似てる）。

出力

Sigmoid gradient evaluated at [-1 -0.5 0 0.5 1]:
[ 0.19661193  0.23500371  0.25        0.23500371  0.19661193]

###7.Gradient Checking（1）～数値計算で微分を計算～
次にコスト関数の勾配がちゃんと機能しているか確認するのですが、デバッグするためにJの差分を取って微分を取るという数値計算で（Back Propagationも数値計算ですが、ここで説明するのははありきたりで貪欲な微分の計算方法で）求めます。つまり、今作ろうとしているのはデバッグ用の関数です。

何をやっているのかというと、$\epsilon$をとても小さい値で取ったときに、1変数の場合

$$J'(\theta)\approx\frac{J(\theta+\epsilon)-J(\theta-\epsilon)}{2\epsilon}$$

を計算します（εをどんどん小さくすれば最終的に微分と一致します）。片側微分ではないのがポイント。ごく普通の方法ですが、偏微分の計算では$J(\Theta)$の計算に2×(Θの変数の個数)回Forward Propagationを行い、さらにいろいろ行列計算もするので遅いです。裏を返せばいかにBack Propagationが効率的なアルゴリズムかがわかります。

## Gradient Checking（Back Propagationのデバッグ）
# Gradient Checking用にJ(θ±ε)を計算する（とても遅い）
def compute_numerical_gradient(cost_wrap_func, theta):
    Theta1 = theta[0]
    Theta2 = theta[1]
    pertub1 = np.zeros(Theta1.shape)
    pertub2 = np.zeros(Theta2.shape)
    numgrad = (np.zeros(Theta1.shape), np.zeros(Theta2.shape))
    e = 1e-4
    # 両側微分を求める関数
    calc_cost = lambda: (cost_wrap_func((Theta1 + pertub1, Theta2 + pertub2)) - cost_wrap_func((Theta1 - pertub1, Theta2 - pertub2))) / 2 / e
    for p in range(np.size(Theta1)):
        index = np.unravel_index(p, Theta1.shape, order="F")
        pertub1[index] = e
        numgrad[0][index] = calc_cost()
        pertub1[index] = 0
    for p in range(np.size(Theta2)):
        index = np.unravel_index(p, Theta2.shape, order="F")
        pertub2[index] = e
        numgrad[1][index] = calc_cost()
        pertub2[index] = 0
    return numgrad

この数値計算で求めた偏微分の値と、Back Propagationで求めた偏微分が一致する（一致しなくてもいいが、相対誤差が極めて小さくなること）を確認します。だからGradient Checkingはデバッグ用の関数なのです。

###8.Gradient Checking（2）～Back Propagationとの比較～
compute_numerical_gradientを定義したことで、ようやくGradient Checkingができます。Gradient Checking用の初期値を与える関数と、Gradient Checkingの本体の関数を定義します。

# デバッグ用の初期値設定
def debug_initialize_weights(fan_out, fan_in):
    W = np.zeros((fan_out, 1+fan_in))
    W = np.sin(np.arange(1, np.size(W)+1)).reshape(W.shape, order='F') / 10
    return W

# Gradient Checking用の関数
def check_nn_gradients(lambda_):
    # デバッグ用パラメーター
    input_layer_size, hidden_layer_size, num_labels, m = 3, 5, 3, 5
    # Theta1,Theta2にランダムな初期値を与える
    Theta1 = debug_initialize_weights(hidden_layer_size, input_layer_size)
    Theta2 = debug_initialize_weights(num_labels, hidden_layer_size)

    # Xの値もデバッグ用のランダムな値にする
    X = debug_initialize_weights(m, input_layer_size-1)
    # yはmod関数で適当に作る
    y = 1 + np.mod(np.arange(1, m+1), num_labels)

    # コスト関数を使ってコストと勾配の計算
    cost, grad = nn_cost_function((Theta1, Theta2), input_layer_size, hidden_layer_size, num_labels, X, y, lambda_)
    cost_warp = lambda Thetas: nn_cost_function(Thetas, input_layer_size, hidden_layer_size, num_labels, X, y, lambda_)[0]
    numgrad = compute_numerical_gradient(cost_warp, (Theta1, Theta2))

    #BackPropagationと数値計算の値の比較
    # 結果をベクトルに変換
    result_grad= np.append(np.ravel(grad[0], order="F"), np.ravel(grad[1], order="F"))
    result_numgrad = np.append(np.ravel(numgrad[0], order="F"), np.ravel(numgrad[1], order="F"))
    for g, n in zip(result_grad, result_numgrad):
        print(n, g)
    print("The above two columns you get should be very similar.")
    print("(Left-Your Numerical Gradient, Right-Analytical Gradient)\n")
    ##相対誤差の計算
    diff = np.linalg.norm(result_numgrad - result_grad) / np.linalg.norm(result_numgrad + result_grad)
    print("If your backpropagation implementation is correct, then")
    print("the relative difference will be small (less than 1e-9).")
    print("Relative Difference: ", diff)

debug_initialize_weightsはデバッグに使う、適当に当てはまるためのΘをsinを使って計算しています。これは特に意味はなさそうなので深く考えなくていいです。

check_nn_gradientsが本体で、正則化のλの値を引数にして、Back Propagationによる勾配が正しいかどうかを見ます。Xもyも同様に適当に入れているのでここも深い意味はありません。

このアルゴリズムのポイントは、（コード見ればわかるのですが）、数値計算で求めた微分の値と、Back Propagationの値をそれぞれベクトルに格納し、差のベクトルのノルムを和のベクトルで割った相対誤差が一定以下（1e-9など）となることを目指します。Andrew Ng先生いわく「これに何度も助かった」とのことなので、何やっているんだかわからなくなりがちなニューラルネットワークを自分で実装しなければいけないときは、確実にやったほうがいいでしょう。組み込みの場合はこういう機能があるか調べてもよくわかりませんでしたが、組み込みを信用する場合は特にいらないかも？

###9.コスト関数のテスト（２）
デバッグ用の関数ができたので、ようやく勾配のテストができます。正直ヘトヘトです。

#λ=0でGradient Checking
print("Checking Backpropagation...")
check_nn_gradients(0)
print()

#λ=3でGradient Checking
print("Checking Backpropagation (w/ Regularization) ... ")
lambda_ = 3
check_nn_gradients(lambda_)
print()

# 同様にでコスト関数もデバッグ
debug_J, debug_grad = nn_cost_function((Theta1, Theta2), INPUT_LAYER_SIZE, HIDDEN_LAYER_SIZE, NUM_LABELS, X_data, y, lambda_)
print("Cost at (fixed) debugging parameters (w/ lambda =", lambda_, ") :", debug_J)
print("(for lambda = 3, this value should be about 0.576051)\n")

出力

Checking Backpropagation...
-0.00927825234864 -0.00927825235799
0.0088991195879 0.00889911959567
（中略）
0.0531542052462 0.0531542052426
0.0465597186272 0.0465597186259
The above two columns you get should be very similar.
(Left-Your Numerical Gradient, Right-Analytical Gradient)

If your backpropagation implementation is correct, then
the relative difference will be small (less than 1e-9).
Relative Difference:  2.57955061714e-11

Checking Backpropagation (w/ Regularization) ...
-0.00927825234864 -0.00927825235799
0.0088991195879 0.00889911959567
（中略）
-0.00452964426634 -0.00452964427015
0.00150048381942 0.00150048381964
The above two columns you get should be very similar.
(Left-Your Numerical Gradient, Right-Analytical Gradient)

If your backpropagation implementation is correct, then
the relative difference will be small (less than 1e-9).
Relative Difference:  2.48150088452e-11

Cost at (fixed) debugging parameters (w/ lambda = 3 ) : 0.57605124695
(for lambda = 3, this value should be about 0.576051)

出力の「Relative Difference」のところに注目してください。λ＝0とλ＝3の2回Gradient Checkingを行っており、それぞれ2箇所ありますが、数値計算での偏微分とBack Propagationでの偏微分の相対誤差がここにあります。どちらもオーダーがe-11（想定はe-9）なので、この実装は正しいことになります。
個々の偏微分の値は、数値計算のものが左、Back Propagationのものが右となります。長いので省略しました。ほぼ一緒ですね。

###10.訓練データのトレーニング
これまではデバッグで、ようやくここからが本番。

## ニューラルネットワークの訓練
import sys, time
# 最急降下法の定義
def gradient_descent(initial_all_theta, cost_warp, eta, maxiter = 100):
    theta_before = initial_all_theta
    for i in range(maxiter):
        J, grad = cost_warp(theta_before)
        theta = copy.deepcopy(theta_before)
        for j in range(len(initial_all_theta)):
            theta[j] = theta[j] - eta * grad[j]
        mes = f"iter = {i}, J = {J}"
        theta_before = theta
        # コンソールを埋め尽くさないように上書きしながら表示
        sys.stdout.write("\r%s" % mes)
        sys.stdout.flush()
        time.sleep(0.01)
    return theta, J

# Θの初期値
initial_Theta1 = rand_initialize_weights(INPUT_LAYER_SIZE, HIDDEN_LAYER_SIZE)
initial_Theta2 = rand_initialize_weights(HIDDEN_LAYER_SIZE, NUM_LABELS)
# 正規化＝1と設定（ここをいろいろ変える）
lambda_ = 1
# ニューラルネットワークの訓練
print("Training Neural Network...")
nn_cost_wrap = lambda Theta: nn_cost_function(Theta, INPUT_LAYER_SIZE, HIDDEN_LAYER_SIZE, NUM_LABELS, X_data, y, lambda_)
nn_Thetas, cost = gradient_descent([initial_Theta1, initial_Theta2], nn_cost_wrap, 2, 500)
print()

前回と同様に最急降下法は自分で実装しました。本当ただの最急降下法なので、収束までにそこそこ時間かかります。この例では500反復（だいたい1分ぐらい）やっています。学習率は大きめ（2ぐらい）がちょうどよかったです。コンソールに表示されるJの値の推移を見てきれいに落ちるように調整すればいいので。

###11.モデルの精度
モデルの精度を評価します。訓練データの一致率です。ここは多クラス分類と同様です。

## 予測と精度
def predict(Theta1, Theta2, X):
    m = X.shape[0]
    num_labes = Theta2.shape[0]
    h1 = sigmoid(np.dot(np.c_[np.ones((m, 1)), X], Theta1.T))
    h2 = sigmoid(np.dot(np.c_[np.ones((m, 1)), h1], Theta2.T))
    p = np.argmax(h2, axis=1)+1
    return p

# 精度
pred = predict(nn_Thetas[0], nn_Thetas[1], X_data)
print("Training Set Accuracy:", np.mean(pred == y) * 100)

出力

Training Neural Network...
iter = 499, J = 0.49221426662002384
Training Set Accuracy: 95.1

精度は95.1%となりました。ニューラルネットワークを1から自分で実装できました！（既に疲労困憊）

###12.ネットワークの可視化
これだけでは何やっているのかわからないので、入力層の係数を画像にして表示してみます。

## ネットワークの可視化
def display_data(X):
    fig = plt.figure(figsize = (5, 5))
    fig.subplots_adjust(hspace=0.05, wspace=0.05)
    for i in range(X.shape[0]):
        ax = fig.add_subplot(5, 5, i+1, xticks=[], yticks=[])
        ax.imshow(X[i, :-1].reshape(20, 20, order="F"), cmap="gray")
    plt.show()

display_data(nn_Thetas[0])#もう二度とニューラルネットワークを自力で書きたくねえ

SAN値が下がりそうな画像がでてきましたが、左上から右下にいたるまで、入力層から隠れ層の25ノードに対する、ピクセルごとの係数を画像で表したものです。画像の黒いところは係数の低いところ、白いところは係数が高いところです。
白い所が集中している箇所は、どのノードでも重要視されているところ、と読み取れます。画像の中心や数字の輪郭が多く出そうな所に集中しているので、感覚的には結構しっくりきます。

###13.オーバーフィッティング
このデータでは、正則化の値を小さし、反復回数を増やすとニューラルネットワークでもオーバーフィッティングすることが確認できるそうです。ただ疲れてしまったので、どうしても暇だという方はやってみてください（組み込みでやったからいいよね）。

以上です。もう二度とニューラルネットワークの実装を自力で書きたくないです。面倒な実装は組み込みにまかせて楽したい！

次回に進む
Coursera Machine LearningをPythonで実装 - [Week6]正則化、Bias vs Variance