社畜エンジニア発掘戦線

駆けだしAIエンジニア

①線形分類問題（２値分類）

第一問

ディープラーニングでトライできる問題として、回帰問題に続き、「分類問題」にチャレンジしていきたいと思います。2次平面上にプロットされ、「0と1」にラベリングされている集合データをディープラーニングで分類してみます。

今回もPYTHONISTA3でコードを書きながら解いていってみましょう。

設問１．入力x1、x2、正解ラベルtが (x1, x2, t) = (2, 2, 0)、(-2, -2, 1) で、それぞれにガウシアンノイズが付加されたデータセットを作成せよ

さて、今回は分類問題にチャレンジです。設問を見ると、入力がx1とx2と増えています。これは２次元平面上にx座標、y座標としてプロットできるということですね。(2, 2)の座標にあたるデータの正解ラベルが「0」、(-2, -2)の方は「1」ということですね。イメージはこんな感じ。
f:id:sutokun:20190204093515j:plain

ではコードで書いていってみます。いつも通りnumpyとmatplotlibをインポートするんですが、今回は2次平面なのでmatplotlibのカラーマップも合わせてインポートします。

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

次にサンプル数（num_of_sam）とガウシアンノイズの広がり具合（std_dv）を定義しておきます。データ数は40、広がり係数は1.8としましょう。

num_of_sam = 40
std_dv = 1.8

では実際に座標を設定します。

group1 = np.array([2,2])
group2 = np.array([-2,-2])

このままではただの点なので、ガウシアンノイズをのせたサンプル数を足します。

group1 = np.array([2,2])+np.random.randn(num_of_sam, 2)*std_dv
group2 = np.array([-2,-2])+np.random.randn(num_of_sam, 2)*std_dv

さらに、group1、group2をひとつのデータセットになるように結合します。このXの1列目がx1、2列目がx2にあたります。

X = np.vstack((group1, group2))

グラフ化してみます、ただノイズの広がりが1.8だと大きすぎて分からないのでグラフは「std_dv=0.5」で描いています。「X[:,0]」と「X[:,1]」がそれぞれx1、x2、グラフのx座標とy座標です。

plt.plot(X[:,0], X[:,1], 'o')
plt.show()

f:id:sutokun:20190204094518p:plain
うまくできています。しかし、今の状態では正解ラベルがないので全て同じ点（同じ色）の集合になっています。それぞれのグループに「0」と「1」を設定します。

t_group1 = np.zeros((num_of_sum, 1))
t_group1 = np.ones((num_of_sum, 1))

Xと同じように結合します。

t = np.vstack((t_group1, t_group2))

データセットの完成です。カラーマップでグラフ化してみましょう。ノイズの広がりは「std_dv=1.8」に戻します。「vmax、vmin」は色を表す値の最大、最小です。正解ラベルの0と1に対応しています。「cmap=cm.bwr」は色の種類を示しています。今回は0が青、1が赤の設定です。

plt.scatter(X[:,0], X[:,1], vmin=0, vmax=1, c=t[:,0], cmap=cm.bwr, marker='o', s=50)
plt.show()

f:id:sutokun:20190206113643p:plain
いい感じです。この問題は「線形分類」なので、直線でこの2グループを分類します。とにかく進めていきましょう。

設問２．入力(x1, x2)に対する重みを(w1, w2)としてシグモイド関数を活性化関数として出力yを求めよ

ディープラーニングの問題解決フローは同じなので、ご察しの通り次はニューラルネットワークの設計です。入力はx1とx2、出力は0か1かをラベリングするyひとつなので、こんな感じ。
f:id:sutokun:20190203061819j:plain

出力には非線形回帰で用いたシグモイド関数を使います。ここで使う関数を活性化関数と呼びます。正解ラベルは0か1なので、この間で推移するシグモイド関数が効果的です。今回は入力が2つになっていますが、非線形回帰問題とかに比べると結構シンプルですね。

まずはシグモイド関数が使えるようにPythonの関数として定義しておきます。

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

次に重みの初期値を設定します。w1とw2をまとめて行列Wで書きましょう。

W = np.random.randn(2,1)

準備が整ったので出力を求めましょう。XとWは線形結合（内積計算）でシグモイド関数を通します。

y = sigmoid(np.dot(X, W))

念のため、print(Y.shape)で出力の型を見ておきましょうか。正解ラベルtと型が一致しているはずです。

print(t) #(80, 1)
print(y) #(80, 1)

これで設問2は終了です。

設問３．クロスエントロピーを損失関数として損失量を求めよ

流れ通り、次は損失関数の設定ですね。回帰問題でスタンダードだった二乗和誤差と異なり、今回はクロスエントロピーを用います。これも場面によっていくつか定義式があるんですが、今回は次式の定義を用いて進めていきたいと思います。

クロスエントロピー：{\displaystyle E = -\sum_{i} \Bigl( t_{i} \log y_{i} +(1-t_{i}) \log (1-y_{i}) \Bigr) }

前回の二乗和誤差関数よりも複雑なので、こちらもシグモイド同様、Pythonの関数で書いておきましょう。

def loss(y, t):
    return -np.sum(np.multiply(t, np.log(y)) + np.multiply((1 - t), np.log(1 - y)))

では設問2で求めたYとデータセットの正解ラベルtからクロスエントロピーを用いて損失量を求めてみます。

E = loss(y, t)

print(E)してみるとちゃんとスカラー値が出てきます、これで設問2は終了です。

設問４．損失関数のW=(w1, w2)による微分式を求めよ

ここの流れももうおなじみです、連鎖率を用いて記述してみます。

{\displaystyle \frac{\partial E}{\partial W}\ = \frac{\partial y }{\partial W}\ \frac{\partial E }{\partial y}\ }
ではクロスエントロピーのy微分から計算してみます。
{\displaystyle \frac{\partial E }{\partial y}\ = \frac{\partial }{\partial y}\ \Bigl(-\sum_{i} \Bigl( t_{i} \log y_{i} +(1-t_{i}) \log (1-y_{i}) \Bigr) \Bigr) = \sum_{i} \Bigl( -\frac{t_{i}}{y_{i}} + \frac{1 - t_{i}}{1 - y_{i}} \Bigr) }
次にyのW微分です。yは入力XとパラメータWの内積をシグモイド関数に通したものです。シグモイド関数の微分を忘れないように。
{\displaystyle \frac{\partial y_{i} }{\partial W}\ = \frac{\partial }{\partial W}\ \Bigl( \sigma (X_{i} \cdot W) \Bigr) = \sum_{i} \Bigl( X_{i}y_{i}(1-y_{i}) \Bigr) }

以上より、

{\displaystyle \frac{\partial E}{\partial W}\ = \sum_{i} \Bigl( X_{i}y_{i}(1-y_{i}) \Bigr) \sum_{i} \Bigl( -\frac{t_{i}}{y_{i}} + \frac{1 - t_{i}}{1 - y_{i}} \Bigr) = \sum_{i} \Bigl( X_{i}(y_{i}-t_{i}) \Bigr) }
気持ちいい鬼約分ですね、かなりシンプルになりました。

ではこの内容をコードに落とします。シグマの方向に注意（axis=0）。

dW = np.sum(X*(y-t),axis=0)

以上で設問4は終了です。

設問６．勾配降下法を用いてパラメータを更新し、イテレーションを30回繰り返した後の出力yを求めよ

ここは今までどおり勾配降下法なので、コードに落とすだけです。とりあえず学習率は0.002としましょう。いつも通り損失量を保存するコードを追記しておきます。

勾配降下法：{\displaystyle P_{new} = P_{old} - \mu \frac{\partial E}{\partial P}\ }

learning_rate = 0.002
E_save = []

num_of_itr = 1
for i in range(num_of_itr):
    y = sigmoid(np.dot(X, W))
    E = loss(y, t)
    E_save = np.append(E_save, E)

    dW = np.sum(X*(y-t),axis=0)
    W = W - learning_rate*np.reshape(dW,(2,1))

dWは要素が2つの行列ですが、型が(1,2)に反転しているので、reshapeでもとのWと同じ型にそろえています。これでパラメータWが更新され、出力を求める準備が整いました。

さて、これで求められる出力をどのように表示するかが問題です。普通に数字だけ出しても「それで？」って感じです。実際に2次元グラフ上に描画してみましょう。色々なやり方があると思いますが、いちばん直感的かな〜と思うのが、2次元グラフの格子点を使って出力を表示してみる方法です。入力データセットXの(x1, x2)をそれぞれ格子点として用意し、各格子点に対して更新された重みWを使って出力yを求めます。
f:id:sutokun:20190203055111j:plain

ちょっとなんのこっちゃという感じですが、とりあえず進めます。

まずは格子点を用意します。どの範囲で表示するかをgrid_range、その範囲に何点の格子点を突っ込むかをresolutionで定義します。-5から5までの範囲を60分割するように設定します。

grid_range = 5
resolution = 60
x1_grid = x2_grid = np.linespace(-grid_range, grid_range, resolution)

これで格子点の外側の目盛りを作成できました。次に中の格子点を作っていきます。numpyのmeshgridを使うと楽に作成できます。

xx, yy = np.meshgrid(x1_grid, x2_grid)
X_grid = np.c_[xx.ravel(), yy.ravel()]

これで格子点の作成が完了しました。一度グラフにプロットしてみましょう。

plt.plot(X_grid[:,0], X_grid[:,1], 'o')
plt.show()

f:id:sutokun:20190203061659j:plain
目がチカチカしますが、ちゃんと格子点がグラフを埋めていますね。

ではこの格子点におけるそれぞれの出力を求めてみます。

Y_grid = sigmoid(np.dot(X_grid, W))

出力はsigmoid関数なので、0から1までの値を取ります。正解ラベルは0か1なのでこのどちらかに入るよう、四捨五入しましょう。

Y_predict = np.around(Y_grid)

ではこの出力をカラーマップで描画してみましょう。同じ色調で描くと色がきついので格子点はalphaで透かしています。

plt.scatter(X[:,0], X[:,1], vmin=0, vmax=1, c=T[:,0], cmap=cm.bwr, marker='o', s=50)
plt.scatter(X_grid[:,0], X_grid[:,1], vmin=0, vmax=1, c=Y_predict[:,0], cmap=cm.bwr, marker='o', s=50,alpha=0.3)
plt.show()

f:id:sutokun:20190206113733p:plain
いい具合に青と赤で分類できました。

設問７．損失量の推移をグラフ化せよ

【解決済】
ここで問題発生です。今回の問題では損失関数にクロスエントロピーを使っていますが、logの計算に不具合が生じます（0で割り算しちゃう系）。これはディープラーニングの本質と関係しているのかなんとも分からず、詳細を調べるのはめんどくさ…非効率なので、今回この設問はパスすることにします。回避方法がひらめいたら修正します（すんません）。【追記】少し調べたらイテレーション5回あたりからエラーが発生します。やはり学習が進んでyの値が更新されてくるとlog計算がだんだん難しくなってくる感じですかね。print(y)してみるとすごい綺麗な「1」が出力されていました。クロスエントロピーの式は「log(1-y)」を含むのでマイナス∞に発散したようです。じゃあなんでこんな綺麗な「1」が出力されるんだ？やっぱり問題が簡単すぎるのか…？【追記】嗚呼ぁぁああああ！！！！凡ミスぅうううう！！！パラメータを更新するときに学習率（learning_rate）をかけ忘れていました。つまり、学習率0.002かけるとか言っといて実は学習率1の天才ニューラルネットワークさんだったわけですね。そりゃ綺麗な「1」も出力されるよ。コレは恥ずかしい。このページのコードは修正しますがこの追記は戒めとして残しておきます。

（いろいろトラブルもありましたが）損失量をグラフ化してみます。

plt.plot(E_save)
plt.show()

f:id:sutokun:20190206114157p:plain
よく収束してくれています。損失量はけっこう残っていますが、もう少しイテレーションを回せば良くなりそうかな？

まとめ

最後にイテレーションの推移でどのように分類が進むかGIFでアニメ化します。（iPadは何回か処理落ちしましたが祈れば何とかなりました）
f:id:sutokun:20190206114402g:plain

問題自体は難しくないですが、データセットとカラーグラフの描画に時間がかかってしまったような感じですね。回帰問題と同様、バイアスを付加していないので分類ラインは原点を中心に回転します。次回はコレにバイアスを付加して自由度のある線形分類にトライしようと思います。

全体コード
github.com

次の問題
②線形分類問題（２値分類：バイアス付加） - 社畜エンジニア発掘戦線

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter

①非線形回帰問題（データフィッティング）

第一問


直線をデータにフィットさせる線形問題はいくつか取り組んできました。ここからは直線ではフィットできないデータにカーブをフィットさせる非線形問題にチャレンジしていきたいと思います。

今回もPYTHONISTA3でコードを書きながら解いていってみましょう。

設問１．ガウシアンノイズを付加したy = sin(x)に準ずるデータセット（インプット：x、正解ラベル：t）を作成せよ

さて、今回から非線形問題ということで、データセットする関数が直線ではなくsinカーブとなりました。とは言え、基本的に解決フローは今までと同じです。そしてデータセットは直線：y=oxのところがy=sin(x)に置きかわるだけです。範囲はsinのだいたい1周期にあたる0から6まででトライしましょう。サンプル数は40にします。ノイズの広がりが大きいとフィッティングが乱れやすいので、標準偏差は0.1で極力狭くしておきます。

import numpy as np
import matplotlib.pyplot as plt

num_of_sam = 40
x = np.random.uniform(0, 6, num_of_sam)
noise = np.random.normal(0, 0.1, num_of_sam)
t = np.sin(x)+noise

グラフにプロットしてみます。

plt.plot(x, t, 'o')
plt.show()

f:id:sutokun:20190204062421p:plain
オッケーですね、設問１は終了です。

設問２．隠れ層が20のニューロンで設計されるニューラルネットワークのパラメータ（w、b）の行列型を求めよ

ココから線形問題よりも難易度がぐっと上がります。ですが、解決フローに則って進めていくことに変わりはないので、ひとつずつトライしていきましょう。

まずは、ニューラルネットワークの設計です。非線形と名のつく問題には隠れ層が必須です。入力と出力の間に位置する隠れ層が20のニューロンで構築されるので、次のような感じになります。ニューロンとは矢印を受けたり出したりする玉のことです。
f:id:sutokun:20190127151059j:plain
いきなりめっちゃパラメータが増えましたね。入力（x）から入るデータセットは隠れ層（h）を通り、出力（y）として出ていきます。隠れ層で実行される計算は今まで通りxw+bですが、ここに非線形関数σを寄与させなければいけません。つまり、隠れ層の出力はh=σ(wx+b)となります。

なぜ非線形関数（曲線）を寄与させなければいけないのかというと、早い話が「直線をいくら足してもいいフィッティングにはならない」からです。wx+bのままでは直線ですが、σに入れてやることで曲線に変換できます。そんな感じの理由です。

ここに用いる関数も場面によって使い分ける必要がありますが、スタンダードなのはシグモイド関数と呼ばれるものです。

この関数のグラフはこんな形をしています。
f:id:sutokun:20190127151027j:plain

このシグモイド関数を２０個（ニューロン数）足し合わせていい感じのフィッティングカーブを作ろう！というのが大まかな流れです。ちなみに出力yには特になんの関数も必要なく線形結合（ただの足し合わせ）で大丈夫です。

では、シグモイド関数をPythonの関数に落としておきましょう。

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

これで「sigmoid(oo)」とするだけでシグモイド関数の計算ができます。

ひと通り準備も終えたところで、パラメータについて考えてみます。隠れ層のニューロン数は20なので重み（w）とバイアス（b）はそれぞれ20個必要です。次に、隠れ層からの出力20個が、ニューラルネットワーク全体の出力（y）の入力となります。それぞれに重み（w）が掛けられて出力（y）へ向かいますが、出力に作用するバイアスはひとつだけです。
f:id:sutokun:20190127155152j:plain

整理すると、最初のwとbは(20, 1)のデータとなり、その次はwが(1, 20)、bは(1, 1)、つまりスカラー値となります。wは前と後ろで型が反転していることに注意です。出力は線形結合なので反転させて内積計算にしなければいけません。

コードに落としましょう。W1とW2が前側のパラメータなので(20, 1)となり、W2とB2が後ろ側のパラメータなのでW2は(1, 20)、B2は(1, 1)となっています。いちいち全部の初期値を設定していると面倒なので、ランダムで値が入るようにしています。

dimention = 20
W1 = np.random.randn(dimention, 1)
B1 = np.random.randn(dimention, 1)

W2 = np.random.randn(1, dimention)
B2 = np.random.randn(1, 1)

これで設問２が終了です。

設問３．ニューラルネットワークの出力（y）を計算せよ

パラメータの初期値も設定できたので、実際に計算を進めていきたいと思います。まずは隠れ層の出力です。

H = sigmoid(x*W1+B1)

型を確認しておきましょう。入力（x）は40個のデータなので40列のデータがニューロン数の20行に収まっています。

print(H.shape) #(20,40)

次は、計算したHを用いて出力yを求めます。線形結合なので内積にしなければいけません。dotは内積の意味です。

y = np.dot(W2, H)+B2

こちらも念のため、型が正しいか確認しておきましょう。入力（x）と一致するはずです。が、この計算だとyは行列として出力されますが、もとのxは配列です。そのためshapeで出力される型は異なります。これでも数値計算は進みますが、最後グラフ化するときに型をそろえる必要があります。

print(x.shape) #(40, )
print(y.shape) #(1, 40)
設問４．二乗和誤差関数から損失量を求め、各パラメータによる損失関数の微分式を導出せよ

次は損失関数です、今回も二乗和誤差関数を用いましょう。

E = np.sum((t-y)**2)

さて、ここからは計算です。損失関数を４つのパラメータ（W1、B1、W2、B2）で微分します。これは逆伝搬計算にあたるのでW2、B2から計算していきましょう。その前に、Eをyで微分した結果はよく使うので先にこれだけ計算して変数δで置いておきます。

{\displaystyle \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial y}\ \sum_{i} (t_{i} - y_{i})^2 = 2 \sum_{i} (y_{i}-t_{i}) = \delta_{i} }

まずはW2、B2から計算です、これまで同様「連鎖率」を使います。

{\displaystyle \frac{\partial E}{\partial W_2}\ = \frac{\partial y}{\partial W_2}\ \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial W_2}\ \Bigl(\sum_{i} (H_{i}W_2)+B_2\Bigr) \cdot \delta_{i} = 2\sum_{i} H_{i} (y_{i}-t_{i}) }
{\displaystyle \frac{\partial E}{\partial B_2}\ = \frac{\partial y}{\partial B_2}\ \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial B_2}\ \Bigl(\sum_{i} (H_{i}W_2)+B_2\Bigr) \cdot \delta_{i} = 2\sum_{i} (y_{i}-t_{i}) }

次にW1、B2の計算です。しかし、その前にシグモイド関数の微分を確認しておきましょう。

{\displaystyle \frac{\partial}{\partial x}\ \sigma(x) = \sigma(x)(1-\sigma(x)) }

ややこしい中身の計算は紙とペンでトライしてもらうとして、結果だけ書くと非常にシンプルな式にまとまります。

このシグモイド関数の微分結果を用いて1層目のパラメータ微分を進めます。こちらも同様プロセスですが、2層分戻らなければいけないのでダブルで「連鎖律」を使います。

{\displaystyle \frac{\partial E}{\partial W_1}\ = \frac{\partial H}{\partial W_1}\ \frac{\partial y}{\partial H}\ \frac{\partial E}{\partial y}\ = \frac{\partial H}{\partial W_1}\ \frac{\partial y}{\partial H}\ \cdot \delta_{i} }

ボリュームがあるのでひとつずつ計算します。

{\displaystyle \frac{\partial H}{\partial W_1}\ = \frac{\partial }{\partial W_1}\ \sigma (x_{i}W_1+B_1) = \sum_{i} x_{i}H_{i}(1-H_{i}) }
{\displaystyle \frac{\partial y}{\partial H}\ = \frac{\partial }{\partial H}\ (H_{i}W_2+B_2) = W_2 }

以上をまとめると、

{\displaystyle \frac{\partial E}{\partial W_1}\ = W_2 \sum_{i} x_{i}H_{i}(1-H_{i}) \cdot \delta_{i} =2 W_2 \sum_{i} x_{i}H_{i}(1-H_{i})(y_{i}-t_{i}) }

こりゃ大変ですね、B1については結果だけ。

{\displaystyle \frac{\partial E}{\partial B_1}\ = W_2 \sum_{i} H_{i}(1-H_{i}) \cdot \delta_{i} =2 W_2 \sum_{i} H_{i}(1-H_{i})(y_{i}-t_{i}) }

なんとか微分式の導出までできたので、コードに落としておきます。

dW2 = 2*np.sum(H*(y-t),axis=1)
dB2 = 2*np.sum(y-t)

dW1 = 2*W2*np.sum(x*H*(1-H)*(y-t),axis=1)
dB1 = 2*W2*np.sum(H*(1-H)*(y-t))

dW2、dW1、dB1はHの項が絡むのでニューロンの数だけ個別にパラメータの更新が必要です。そのため、axis=1を加えて行要素のみにsumが効くように調整しています。…なんですが、バイアスの変化はとてもシビアに出力へ影響します。バイアスを個別に更新すると損失関数がうまく収束に向かわない傾向があるのでdB1からはaxis=1を除いて、全ての要素のsumで更新されるようにしています。

長くなりましたが、これで設問４が終了です。

設問５．降下法を用いてパラメータの更新式を導出せよ

微分式が求められたのでココはもう簡単ですね。Pはパラメータ（w、b）の意味です。

降下法：{\displaystyle P_{new} = P_{old} - \mu \frac{\partial E}{\partial P}\ }

この式に先の微分式を代入するだけ。

{\displaystyle W_{2new} = W_{2old} - \mu \frac{\partial E}{\partial W_{2}}\ = W_{2old} - \mu \Bigl(2 W_2 \sum_{i} x_{i}H_{i}(1-H_{i})(y_{i}-t_{i})\Bigr) }
{\displaystyle B_{2new} = B_{2old} - \mu \frac{\partial E}{\partial B_{2}}\ = B_{2old} - \mu \Bigl(2 W_2 \sum_{i} H_{i}(1-H_{i})(y_{i}-t_{i})\Bigr) }

{\displaystyle W_{1new} = W_{1old} - \mu \frac{\partial E}{\partial W_{1}}\ = W_{1old} - \mu \Bigl( 2\sum_{i} H_{i} (y_{i}-t_{i})\Bigr) }
{\displaystyle B_{1new} = B_{1old} - \mu \frac{\partial E}{\partial B_{1}}\ = B_{1old} - \mu \Bigl( 2\sum_{i} (y_{i}-t_{i})\Bigr) }

ではコードに落とします。学習率（learning_rate）はこの後定義するので、とりあえずコードを書いてしまいましょう。

W1 = W1-learning_rate*dW1.T
W2 = W2-learning_rate*dW2
B1 = B1-learning_rate*dB1.T
B2 = B2-learning_rate*dB2

W1とB1は計算される出力は同じなのですが、型が反転します。そのためTで転置させています。（この辺は本質ではないのでコーディングが上手な方はうまくアレンジして下さい）

設問５が終了です。

設問６．学習率を0.002として、イテレーションを3000回繰り返した後の出力を求めよ

ようやくゴールが見えてきました。イテレーションの数がハンパなく増えましたが、非線形の収束にはこれくらい必要なようです。ではいつもも通りforを使って繰り返しのコードを書いてみましょう。どうせこの後にコスト関数のグラフ化をするので、先に損失量を保存するコードも入れておきます。

learning_rate = 0.002
E_save = []

num_of_itr = 3000
for i in range(num_of_itr):
    #forward propagation
    H = sigmoid(x*W1+B1)
    y = np.dot(W2,H)+B2
    E = np.sum((t-y)**2)
    E_save = np.append(E_save,E)

    #differential
    dW2 = 2*np.sum(H*(y-t),axis=1)
    dB2 = 2*np.sum(y-t)
    dW1 = 2*W2*np.sum(x*H*(1-H)*(y-t),axis=1)
    dB1 = 2*W2*np.sum(H*(1-H)*(y-t))

    #back propagation
    W1 = W1-learning_rate*dW1.T
    W2 = W2-learning_rate*dW2
    B1 = B1-learning_rate*dB1.T
    B2 = B2-learning_rate*dB2

では、ご待望のグラフ化です。描画するのはカーブなのでx軸は200点と細かめに刻んでおきます。最終的に出力されたW1、B1、W2、B2を使って描画用の入力（X_line）を順伝搬計算します。が、このままではグラフが描けません。上の方でちょろっと書いたんですが、この順伝搬計算の出力は行列で出力されますが、入力は配列です。そのためそれらの型をそろえないとグラフが描けません。出力はravelを使って配列にしてしまいます。（ここも本質ではないので、コーディングが上手な方はアレンジして下さい）

X_line = np.linspace(0, 6, 200)
H_line = sigmoid(X_line*W1+B1)
Y_line = np.ravel(np.dot(W2,H_line)+B2)

plt.plot(x, t, 'o')
plt.plot(X_line, Y_line, 'o')
plt.show()

f:id:sutokun:20190204062717p:plain
尻尾のあたりが微妙ですが頑張ってフィッテイングしてくれました感じは伝わります。

設問７．損失量のイテレーション回数による推移をグラフ化せよ

いつも通り、損失関数で求められる損失量を確認してみましょう。既に損失量を保存するコードは書いているので、それをプロットするだけです。

plt.plot(E_save)
plt.show()

f:id:sutokun:20190204062803p:plain
大きな問題なく、収束に向かってくれています。尻尾のあたりがうまくフィットしないのに損失量の推移は落ち着いてしまってこれ以上フィッティングが進みません。んん、これはどうしたもんかな〜という感じです。

まとめ

非線形問題の第一問が完了しました。同じくGIFでアニメ化します。イテレーションは3000回でそんなに画像出力はできないので、100回ごとに画像出力させて全部で30枚をアニメ化しています。
f:id:sutokun:20190204063008g:plain

大きな問題解決フローは今までと変わりませんが、なにせパラメータの数が一気に増えたので頭が爆発しそうになります。ニューラルネットワークを使った非線形フィッティングはあくまでシグモイド関数の足し合わせで表現されるので、完ぺきにきれいな曲線になることはありません。じゃあシグモイド関数以外の関数を使ったらどうなるんだ、とか面白そうです。また問題にして考えてみたいと思います。隠れ層の出力をグラフ化してシグモイド関数たちがどう動くのか、とかも面白そう。

全体コード
github.com

次の問題
①線形分類問題（２値分類） - 社畜エンジニア発掘戦線

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter

③線形回帰問題（損失関数の可視化）

第三問


今回も線形回帰問題ですが、イテレーションごとの損失量の推移に着目したいと思います。前回の最後にちょろっと書きましたが、まだパラメータの数が多くないのでいろいろとグラフ化して視覚的に変化を追うことができます。最初にトライしたバイアスなしの問題をベースに可視化問題にチャレンジしてみます。

参照問題 -> ①線形回帰問題（データフィッティング） - 社畜エンジニア発掘戦線

今回もPYTHONISTA3でコードを書きながら解いていってみましょう。

設問１．ガウシアンノイズを付加したy = 3xに準ずるデータセットを作成し、バイアスなしでフィッティング直線を出力せよ

1回目の問題と同じです、ガバっとコピペしちゃいましょう。

import numpy as np
import matplotlib.pyplot as plt

#dataset
num_of_sam = 20
x = np.random.uniform(0, 2, num_of_sam)
noise = np.random.normal(0, 0.2, num_of_sam)
t = 3*x+noise

#initial_setting
w = 1
learning_rate = 0.004
E_save = []

#iteration
num_of_itr = 25
for i in range(num_of_itr):
    y = w*x
    E = np.sum((y-t)**2)
    E_save = np.append(E_save, E)
    w = w - learning_rate*2*np.sum(x*(y-t))

#plot
X_line = np.linspace(0, 2, 5)
Y_line = w*X_line

plt.plot(x, t, 'o')
plt.plot(X_line, Y_line)
plt.show()

f:id:sutokun:20190204033844p:plain
既にやっているので大丈夫ですね。

設問２．パラメータ：重み（w）に対する損失量の推移をグラフにプロットし、最小値となるwを求めよ

重みはイテレーションごとに都度更新されますが、点ではなく広い範囲でwを変化させて損失量がどのように推移するかを可視化します。損失関数にはシグマが含まれているのでwの関数としてさらっと書くのは困難です。wを繰り返しひとつずつ変化させてそれぞれの損失量を求めていきましょう。繰り返しと言えばfor文ですね。

まずは変化させるw（weight：x軸）とそれに対応する損失量（E_cost_weight：y軸）を保存するリストを設定しておきましょう。

weight = []
E_loss_weight = []

次にweightを何点（num_of_w）刻むか、どれぐらい細かく刻むか（resolution）を設定します。今回は0から10まで0.1ずつ45点刻めるように設定しましょうか。

num_of_w = 45
resolution = 0.1

ではforの中身を書いていきましょう。forが1回計算するごとに損失量を計算し、その時のweightとE_loss_weightを保存していきます。

for j in range(num_of_w):
    y_weight = (j*resolution)*x
    E_weight = np.sum((y_weight-t)**2)

    #strage data
    weight = np.append(weight, j*resolution)
    E_loss_weight = np.append(E_loss_weight, E_weight)

ではグラフ化してみます。

plt.plot(weight, E_loss_weight)
plt.show()

w=3.0で最小値をとる2次関数となりました。確かにy=3xをベースにデータセットを用意したので、最小値がw=3.0になるのは納得ですね。
f:id:sutokun:20190204033919p:plain
これで設問２は終了です。

設問３．設問２のグラフにイテレーションごとの損失量をプロットせよ

設問２で損失量をプロットしたグラフに、実際にイテレーションを回して得られる損失量がどう推移するのか同じグラフにプロットしてみましょう。

まず、同じグラフにプロットするために、イテレーションごとにパラメータ（w）も保存しなければいけません。設問１のコードに追記します、と言っても2行だけですが。

#initial value
w = 1
learning_rate = 0.02
E_loss = [] 
w_itr = [] #空のリストを用意

#iteration
num_of_itr = 10
for i in range(num_of_itr):
    y = w*x
    E = np.sum((y-t)**2)
    E_loss = np.append(E_cost, E) 
    w_itr = np.append(w_itr, w) #イテレーションごとのwを保存
    w = w - learning_rate*2*np.sum(x*(y-t))

では準備が整ったのでプロットしてみましょう。

num_of_w = 45
resolution = 0.1

plt.plot(w_itr, E_save, 'o')
plt.plot(weight, E_loss_weight)
plt.show()

f:id:sutokun:20190204033959p:plain
これを見ると、赤丸（イテレーションごとの損失量）が損失関数（二次曲線）上を転がって最小値（w=3.0）へ向かっていることが分かります。

これで設問３は終了です。

まとめ

これで第三問が完了しました。今回も最後にイテレーションごとに損失量が推移する様子をGIFアニメ化しました。（これはまだiPadでも耐えれる処理量）
f:id:sutokun:20190204034916g:plain
重みパラメータと損失量がどのような関係で推移していくのかを視覚的に捉えることは重要です。今後、損失量が収束せずに頭が爆発することもあると思いますが、そんなときに一度基本に立ち返って「どのように損失関数が収束するのか」を捉える手がかりになればと思います。

全体コード
github.com

次の問題
①非線形回帰問題（データフィッティング） - 社畜エンジニア発掘戦線

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter

②線形回帰問題（バイアス付加）

第二問


前回のに引き続き線形回帰問題にチャレンジします。
前回の問題 -> ①線形回帰問題（データフィッティング） - 社畜エンジニア発掘戦線

今回もPYTHONISTA3でコードを書きながら解いていってみましょう。

設問１．ガウシアンノイズを付加したy = 5x+2に準ずるデータセット（インプット：x、正解ラベル：t）を作成せよ

前回とほとんど同じですが、直線に切片が追加されています。まず、numpyとmatplotlibをインポートします。

import numpy as np
import matplotlib.pyplot as plt

次に、y = 2x+5を描くためのxをセットします。今回はランダムに0から4の範囲でサンプル数は30にしましょう。（前回と数字を変えたことにあんまり意味はありません、最後グラフ化したとき極力分かりやすいように〜ぐらいです）

num_of_sam = 30
x = np.random.uniform(0, 4, num_of_sam)

最後に、y = 2x+5を書いて終わりです。今回も正解ラベルをtとしているのでt = 2x+5としましょう。

t = 2*x+5

前回と同じく、この正解ラベルにガウシアンノイズを乗せます。正規（ガウス）分布に沿った「平均0、標準偏差0.5」のノイズをサンプル数だけ生成します。このノイズをt = 2x+5に足しましょう。

noise = np.random.normal(0, 0.5, num_of_sam) 
t = 2*x+5+noise

これで設問１は終了です、グラフ化してみましょう。

plt.plot(x, t, 'o')
plt.show()

f:id:sutokun:20190204015905p:plain
ちょっと分かりにかもくいですが、切片がプラス5されています。

設問２．初期パラメータ：重み(w)を0.2、バイアス（b）を0として、出力を計算せよ

ではニューラルネットワークを設計します。今回も線形問題であることに変わりはないので「入力層」と「出力層」のみで構築できます。しかし、ここにバイアスを付加しなければいけません。なので設計するニューラルネットワークはこんな感じ。
f:id:sutokun:20190127090406j:plain

切片としてバイアス（b）が付加されています、大きく変わってないですね。ではコードに落としてみます。

w = 0.2
b = 0
y = w*x+b

設問２は完了です。

設問３．二乗和誤差関数を損失関数（E）として、損失量の値を求めよ

今回も損失関数に二乗和誤差関数を用いて損失量を計算します。前回も登場しましたが、二乗和誤差関数は次のような式で表されます。

二乗和誤差関数： \displaystyle E = \sum_{i} (y_{i} -t_{i} )^2 

E = np.sum((y-t)**2)

ここの記述は前回と全く同じです、バイアスが付加されても損失関数の式は変わりません。

設問４．勾配降下法を用いて、パラメータ更新式を導出せよ

ここから少し計算が必要です。前回は損失関数をパラメータwで微分してwの更新式を求めました。バイアスも同様のプロセスで更新式を求めにいきます。まずは損失関数のパラメータ微分から。

{\displaystyle \frac{\partial E}{\partial b}\ = \frac{\partial }{\partial b}\ \sum_{i} (t_{i} - y_{i})^2 }
「連鎖律」を使います。
{\displaystyle \frac{\partial E}{\partial b}\ = \frac{\partial y}{\partial b}\ \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial b}(x_{i}w+b)\ \frac{\partial }{\partial y}\ \sum_{i} (t_{i} - y_{i})^2 }
微分を進めましょう。
{\displaystyle \frac{\partial E}{\partial b}\ = 1 \cdot \Bigl(-2 \sum_{i} (t_{i} - y_{i})\Bigr) = 2 \sum_{i} (y_{i} -t_{i} ) }
wx+bのbでの微分は1なので、wみたいにややこしいこと考えなくていいですね。

では、パラメータの更新式を求めます。更新方法は今回も降下法です。

{\displaystyle b_{new} = b_{old} - \mu \frac{\partial E}{\partial b}\ }

損失関数をbで微分した式を代入します。

{\displaystyle b_{new} = b_{old} - \mu \Bigl(2 \sum_{i} (y_{i} -t_{i} )\Bigr) }

ちなみに、前回導出したｗの更新式はこちら。

{\displaystyle w_{new} = w_{old} - \mu \Bigl(2 \sum_{i} x_{i}(y_{i} -t_{i} )\Bigr) }

これで設問４は終了です。

設問５．学習率0.0018でイテレーションを30回転させ、その後の出力を計算せよ

今回の学習率は0.0018です。バイアスを付加するとなかなか厄介で、損失関数が発散（過学習）しやすくなります。そのため学習率を極力下げる分、ちょっと時間がかかりそうなのでイテレーション回数を増やしました。wも更新しなければいけないのでお忘れなく。

learning_rate = 0.0018
w = w - learning_rate*2*np.sum(x*(y-t))
b = b - learning_rate*2*np.sum(y-t)

これでパラメータの更新が完了しました。順伝搬から逆伝搬までの部分をインテンドさせて、forの中に入れてしまいます。learning_rateの定義は外に出しましょう。

w = 0.2
b = 0
learning_rate = 0.0018

num_of_itr = 30
for i in range(num_of_itr):
    y = w*x+b
    E = np.sum((y-t)**2)
    w = w - learning_rate*2*np.sum(x*(y-t))
    b = b - learning_rate*2*np.sum(y-t)

これでニューラルネットワークの構築が完成しました。うまくいっていれば最初に設定したプロットデータに沿っていい感じに直線がフィッティングされているはずです。

X_line = np.linspace(0, 2, 5)
Y_line = w*X_line+b

plt.plot(x, t, 'o')
plt.plot(X_line, Y_line)
plt.show()

f:id:sutokun:20190204015955p:plain
概ね合ってそうですがもうちょっと頑張る余地はあるかな、イテレーション回数を増やせば微修正が進んでいきます。

これで設問６は終了です。

設問６．損失量のイテレーション回数による推移をグラフ化せよ

前回と同様、各イテレーションで損失量を保存するコードを追記します。

w = 0.2
b= 0
learning_rate = 0.0018
E_save = [] #空のリストを用意

num_of_itr = 30
for i in range(num_of_itr):
    y = w*x+b
    E = np.sum((y-t)**2)
    E_save = np.append(E_save, E) #E_saveにEの値を保存
    w = w - learning_rate*2*np.sum(x*(y-t))
    b = b - learning_rate*2*np.sum(y-t)

では、最後に損失量を保存したリストをグラフ化します。

plt.plot(E_save)
plt.show()

f:id:sutokun:20190204020026p:plain
バイアスを追加したことで1回目の損失量は大きな値になっています。最初の値がでかすぎて30回目がどうなってるか見えにくいですね、拡大してみます。
f:id:sutokun:20190204020049p:plain
イテレーションをもう少し回すといいフィッティングになりそうですね。こんな感じでイテレーション回数が不足しているかどうかもこのグラフから確認できます。とは言え、入力数も増やしてノイズ幅も広げたので収束してもそこそこの損失量は残りそうですが。

まとめ

これで第二問が完了です。直線がプロットにフィットしていく様子をアニメ化しました。（1枚目の画像を出力し忘れてダイナミックな動きがなくなってしまいました…）
f:id:sutokun:20190204020826g:plain

前回と同じような問題をベースにバイアスを付加しました。バイアスを付加することで直線に自由度が生まれ、原点束縛から開放されました。まだパラメータも少なく、十分グラフ化できるレベルの次元ですが、ニューラルネットワークが複雑化すると重みやバイアスも行列化され、可視化できなくなります。そうなるとこれらパラメータのもともとの役割が何だったのか見えにくくなるので、簡単なモデルできちんと確認しておくのは良いことかと思います。

[ note ]
ここまでGIFアニメ用の画像出力もPYTHONISTA3でやってたんですが、処理落ちが激しくなってきました。これ以上iPadにムチを打つのも非効率なんで、次から画像出力はPC（Mac）でやろうと思います、ちょっとグラフの雰囲気が変わります。画像出力以外の計算はPYTHONISTA3でまだ問題なくできそうです。

全体コード
github.com

次の問題
③線形回帰問題（損失関数の可視化） - 社畜エンジニア発掘戦線

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter

PYTHONISTA3のクセ（グラフの描画）

おつかれまです。

PYTHONISTA3はmatprotlibをインポートすることでグラフ描写が可能です。しかし、このアプリ上でグラフが描画されるときのクセを知ってるとイライラせずに済むこともたくさんあります。

ひとつのプラグラムで複数回グラフ描写するときはfigureで分ける：2019-01-27

ふだんはPCでPythonを書いているのであんまり気にしてなかったんですが、PCだと特にfigureを設定しなくても、plt.show()するたびにその直前の設定で毎度グラフを出力してくれるんですね。ありがたく（気づかず）この仕様にあやかっていました。figureとはグラフを描画するためのレイヤー？絵のキャンパスみたいなものです。


ただ、PYTHONISTA3だとそうもいかなくて、きちんと出力したいグラフごとにfigureを設定してあげないとplt.plot()したデータが同じグラフに延々と重ねられてしまいます。
f:id:sutokun:20190128062829j:plain

なので、グラフを描きたいときは都度figureを入れてあげましょう。

#graph 1
plt.figure()
plt.plot(xxxxxx)
plt.show()

#graph 2
plt.figure()
plt.plot(yyyyyy)
plt.show()
グラフ画像の保存ができる：2019-01-27

PYTHONISTA3は機能が多くないし、そこまで期待してなかったんですが、普通にできました。

plt.savefig('NAME')

保存先はカレントフォルダ（走らせているプログラムがあるフォルダ）です。iCloud上にこのフォルダを作れば出力した画像をPCで編集してアニメーションにしたり、と結構できることが広がります。

iPad（PYTHONISTA3）のコード
f:id:sutokun:20190128064500j:plain

保存された画像はPC（Mac）と同期している
f:id:sutokun:20190128064627p:plain

まあ、アニメーション作成だとiPad上でもできますが。すごい、そうなるとiPadだけで結構なことできるな。

処理落ちする：2019-02-03

「なんでもできるやん！」と意気揚々でしたが、iPadの処理限界はそこまで高くありませんでした。特にグラフの出力を絡めて計算するときにアプリが落ちます。グラフ画像の出力もPYTHONISTA3でやっていましたが、無理そうならPC使います。


元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter


続く（追記予定）

PYTHONISTA3のクセ（全般）

おつかれまです。

PYTHONISTA3はiOS上で走るPythonのプログラミングアプリということで重宝しています。ただ、やはりPCと比較して「できること、できないこと」が多々あります。PYTHONISTA3を使ってみる中でそういう「クセ」をまとめておきたいと思います。

※グラフの描画に関しては別でまとめています。

標準機能のまま使う

ネットとかでPYTHONISTA3の記事をみてみると、「Stash」というシェルを導入して機能を拡張させよう！という話題が多かったんですが、かく言う私もトライしました。しかし、コードを書いていて「なんでそこで詰まるねん」というようなバグが出たり、そもそもそんなにStashの恩恵を受ける場面がなかったので、一度削除、再インストールしてそれ以来何も触らずに使っています。

numpyとmatplotlibは標準で使えるので、ディープラーニングを自作するのには問題ありません。

iCloudと同期している

これすごい便利だな〜と思います。使える機能の差は置いといて、PCのキーボードとエディタソフトのおかげで「コードを書く効率さ」はやはりPCの方が上です。ただ、外出先、電車、飛行機の中とか限られた環境で使うにはiPadがとても便利です。PCとiPadが同期していれば、途中までPCで書いたコードをiPadで引き継げます。

この機能はホントvery goodですね。

謎「NameError」が出る

発生トリガーがよく分からないんですが、「random_sampleが定義されてない」って文句を言われます。
f:id:sutokun:20190130092210j:plain

いや、定義も何も、numpyしかimportしてないし。しょうがないのでアプリを再起動するとなんか直ります。よく分からん。

MacのターミナルでiCloudのPYTHONISTAディレクトリへいく：2019-02-03

PYTHONISTA3がiCloudと良い連携をしてくれることは分かったので、PC環境でも動作確認をしたい。そのためにはターミナルでiCloudのディレクトリまで向かわなければいけません。その行き方をメモしておきます。

/Users/FUTOSHI/Library/Mobile\ Documents/iCloud~com~omz-software~Pythonista3/Documents/neuralnetwork

最後の「neuralnetwork」は自分で作ったディレクトリです。気を付けないと行けないのは途中の「Mobile Documents」、半角スペースをそのまま打つと認識してくれないので「Mobile\ Documents」とする必要があります。

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter


続く（追記予定）

①線形回帰問題（データフィッティング）

第一問


とっつきやすいように高校数学のような形式の問題にしてみました。PYTHONISTA3でコードを書きながら解いていってみましょう。

設問１．ガウシアンノイズを付加したy = 3xに準ずるデータセット（インプット：x、正解ラベル：t）を作成せよ

ガウシアンノイズの乗せ方はいったん置いといて、y=5xのデータを作ってみましょう。まず、numpyとmatplotlibをインポートします。

import numpy as np
import matplotlib.pyplot as plt

次に、y = 3xを描くためのインプット：xをセットします。ランダムに0から2の範囲でサンプル数は20にしましょう。

num_of_sam = 20
x = np.random.uniform(0, 2, num_of_sam)

最後に、y = 3xを書いて終わりです。ただし、今回は正解ラベルをtとしているのでt = 3xとしましょう。

t = 3*x

では、この正解ラベルにガウシアンノイズを乗せます。正規（ガウス）分布に沿った「平均0、標準偏差0.2」のノイズをサンプル数だけ生成します。このノイズを先ほどのt = 3xに足しましょう。

noise = np.random.normal(0, 0.2, num_of_sam) 
t = 3*x+noise

これで設問１は終了です、グラフ化してみましょう。「'o'」はおっきい丸をプロット、という意味です。

plt.plot(x, t, 'o')
plt.show()

ノイズの平均や標準偏差を調整するとデータのばらつき具合が変わるので確認してみてください。ちなみに、xの値はプログラムを走らせるごとに与えられる乱数なので、毎回同じプロットにはなりません。
f:id:sutokun:20190203125702p:plain

設問２．初期パラメータ：重み(w)を1として、出力を計算せよ

さて、本題です。出力を求めるにはニューラルネットワークを設計しなければいけません。線形問題に関しては「入力層」と「出力層」のみで構築できます。これさえ分かっていれば、今回の対象となる関数はy=3x（線形）なのでこの設計像が浮かび上がります。
f:id:sutokun:20190126143346j:plain

非常にシンプルですね、今回の重み（w）は1です。この重みとは線形回帰では直線の「傾き」を表します。つまり、出力はこの重みと入力の積：y = wxとなります。この計算が「順伝搬」にあたります。

w = 1
y = w*x

では、x、t、yのデータ形状を確認してみましょう。

print(y.shape) #(20,)
print(x.shape) #(20,)
print(t.shape) #(20,)

すべて一致しているはずです。今回の問題はシンプルなので「あたりまえやん」となりますが、今後ニューラルネットワークが複雑化するとこれらのデータ形状が行列化します。きちんと入力、正解ラベル、出力の型を整理しておかないと、途中で型が合わずに計算ができなくなることが多々あります、こうなったら頭が爆発します。気にかけておいてください。

設問３．二乗和誤差関数を損失関数（E）として、Eの値を求めよ

損失関数とは、得られた出力yが正解ラベルtとどれくらい離れているか（誤差：損失量）を計算するための関数です。よく用いるものとして「二乗和誤差関数」、「クロスエントロピー」の2種類があります。今回は二乗和誤差関数を用いて損失量を計算します。ちなみに、広義で損失関数をコスト関数とも呼びます。コスト関数とは損失関数に過学習を防ぐ対策が施された関数を指します。過学習とは字の通り、勉強し過ぎてゴールを見失うことです、人間でも勉強しすぎてぶっ飛んじゃう方いますよね。。


二乗和誤差関数： \displaystyle E = \sum_{i} (y_{i} -t_{i} )^2

シグマの前に1/2がつくこともありますが、これは微分したときに係数が消えてくれるようにするテクニックです。この関数の大事なところは極小となる位置なので、「前に係数が付いても位置は変わんないよね」って発想のテクニックです。この数式をコードで書いてみます。

E = np.sum((y-t)**2)

ちなみに、数式の通りですがEはスカラー値です。printで確認してみてください。

print(E)

確かに、ひとつの値のみが出力されます。

設問４．勾配降下法を用いて、パラメータ更新式を導出せよ

ここからがディープラーニングの真髄と言えるでしょうか。先ほど計算した損失関数の値が0に近づくように、パラメータ：重み（w）の値を更新していきます。それはつまり、出力yと正解ラベルtの値がどんどん近くなっていく、ということを意味します。では、どのように更新するのかというと、「微分値」を用います。損失関数（E）を更新したいパラメータ（w）で微分します。詳しい内容は難しい専門書に譲りますが、要はパラメータ（w）をちょっと変化させると損失関数（E）がどれくらい変化して、0に近づいてくれるかということを表しています。とりあえず、その通りに書いてみましょう。

{\displaystyle \frac{\partial E}{\partial w}\ = \frac{\partial }{\partial w}\ \sum_{i} (t_{i} - y_{i})^2 }
「うお、Eの中にwがないから微分できねぇ！」となるんですが、ここで「連鎖律」というテクニックを使います。
{\displaystyle \frac{\partial E}{\partial w}\ = \frac{\partial y}{\partial w}\ \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial w}x_{i}w\ \frac{\partial }{\partial y}\ \sum_{i} (t_{i} - y_{i})^2 }
「約分できて元の形に戻るやん！」という感じで覚えてもらえばいいと思います。

微分を進めましょう。

{\displaystyle \frac{\partial E}{\partial w}\ = \Bigl( \sum_{i} x_{i}\Bigr) \Bigl(-2 \sum_{i} (t_{i} - y_{i})\Bigr) = 2 \sum_{i} x_{i} (y_{i} -t_{i} ) }
「おいおい微分したらなんでシグマがでてくんねん！」とかいろいろツッコミがあるんですが、これであってるらしいです。数学的にややこしいところなんですが、こういうもんだと飲み込んでしまいましょう。

では、具体的な微分式を得られたところで、実際にパラメータを更新します。パラメータをを更新する手法も何パターンかあるんですが、今回は「勾配降下法」を利用します、一番スタンダードな手法ですね。他にも「モーメンタム法」や「AdaGrad法」などいろいろあるんですが、これらは問題の複雑さによって使い分けられます。「降下法」がいちばんシンプルですが、欠点もあります。ただし、今回の線形回帰のような問題には適しています。降下法では次のような式で更新式が与えられます。

{\displaystyle w_{new} = w_{old} - \mu \frac{\partial E}{\partial w}\ }

損失関数をwで微分した式は先ほど求めたので、この式に代入します。

{\displaystyle w_{new} = w_{old} - \mu \Bigl(2 \sum_{i} x_{i} (y_{i} -t_{i} )\Bigr) }

突然でてきたμという係数ですが、これは学習率と呼ばれ、0から1の間で調整します。μの値が大きければ「頭のいい（学習が早い）AI」となりますが、学習が早すぎるとおかしなところにゴールしてしまうことがあります（過学習）。逆にμの値が小さければ、着実に学習を進めていきますが、時間がかかる（なんとも人間らしい）。μはトレードオフの関係を作るので、適切な値を探す必要があります。残念ながら、これは実際の人間の経験と勘によるところがあります。

これで設問４は終了です。

設問５．学習率0.004でイテレーションを25回転させ、その後の出力を計算せよ

イテレーションとはパラメータの更新を指します。パラメータの更新式も出せたので、実際にコードに書いて更新してみます。

learning_rate = 0.004
w = w - learning_rate*2*np.sum(x*(y-t))

これでパラメータの更新が完了しました。この計算が「逆伝搬（バックプロップ）」にあたり、まさにプログラムが「学習」します。さて、パラメータの更新は完了しましたが、1回の更新で学習は終わりません。この問題では、ここまでの順伝搬と逆伝搬を25回繰り返さなければいけません。繰り返しと言えば「for文」ですね。

順伝搬から逆伝搬までの部分をインテンドさせて、forの中に入れてしまいます。learning_rateの定義は外に出して置きましょう。

w = 1
learning_rate = 0.004

num_of_itr = 25
for i in range(num_of_itr):
    y = w*x
    E = np.sum((y-t)**2)
    w = w - learning_rate*2*np.sum(x*(y-t))

損失関数だけ浮いてる感じがしますが、次の設問で使うのでそのままforの中に入れておきます。では、25回更新を終えたwを使って最初のグラフに直線を出力してみましょう。うまくいっていれば最初に設定したプロットデータに沿っていい感じに直線がフィッティングされているはずです。

X_line = np.linspace(0, 2, 5)
Y_line = w*X_line

plt.plot(x, t, 'o')
plt.plot(X_line, Y_line)
plt.show()

f:id:sutokun:20190203125837p:plain

いい感じですね、これで設問６は終了です。

設問６．損失関数のイテレーション回数による推移をグラフ化せよ

最終設問です、損失関数は出力と正解ラベルの誤差を表します。イテレーションを回すごとにこの誤差をできるだけ小さくするようにパラメータは更新されていきます。先ほど、forの中にEの計算を残しました。イテレーションが回るごとにEも再計算されていきますが、このままでは計算されっぱなしで上書きされ続けます。イテレーションごとにEをリスト保存できるようにコードを追記しましょう。

w = 1
learning_rate = 0.004
E_save = [] #空のリストを用意

num_of_itr = 25
for i in range(num_of_itr):
    y = w*x
    E = np.sum((y-t)**2)
    E_save = np.append(E_save, E) #E_saveにEの値を保存
    w = w - learning_rate*2*np.sum(x*(y-t))

では、最後に損失量を保存したリストをグラフ化します。

plt.plot(E_save)
plt.show()

f:id:sutokun:20190203125809p:plain


イテレーション25回でそれなりに収束していますね。損失量がどれくらいのイテレーション回数で収束するかを確認するために、よくこのグラフ描画をします。ニューラルネットワークがより複雑になってくると損失量の推移も複雑になります。これはまた追々問題を解きながら見ていきましょう。

まとめ

これで第一問が完了です。最後に直線がプロットにフィットしていく様子をGIFで作成しました。画像の保存までPYTHONISTA3で行い、GIFの作成は別のソフトを使いました。
f:id:sutokun:20190203130023g:plain

この問題がニューラルネットワークの入り口の入り口で、「機械学習がやりたくても何から始めりゃいいのか初学者」はまずはココから始めていけば理解が広がっていくのではないかと思います。今回の出力はy=xwなので切片がありません。つまり、どれだけデータにフィットしようと頑張っても原点に縛られています。このニューラルネットワークを原点から開放してあげるには傾き：重み（w）と同様に切片：バイアス（b）も設定してあげなければいけません。次はバイアスを付加した線形問題にチャレンジします。

[ note ]
グラフ出力のコードは最小限のものしか書いていません。写真のグラフは見やすいようにラベルなど加工しています。詳しくは下リンクのGithubに全体コードを載せています。PYTHONISTA3はグラフの描画にクセがあるかな〜と感じますが、今のところ問題はないですね。

全体コード
github.com

次の問題
②線形回帰問題（バイアス付加） - 社畜エンジニア発掘戦線

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter