０からのニューラルネットワーク：Python実装

# はじめに
　この記事の執筆は、機械学習エンジニアとしてキャリアをスタートさせたばかりの筆者が、ニューラルネットワーク（NN）って一体なんぞや？を独学で勉強していくうちに、「皆が思っているほどNNって複雑ではないのでは」と思うようになったことに端を発している。

　実際に勉強していくうちに、NN（の基本的な概念や構造）は実はそこまで難しいものではない、という実感を得た。

　この記事の目的は、NNに関して専門的で正確で厳密な説明をすることではなく、初心者の初心者による初心者のための、あくまで理解を助けるためのプラグマティックな知識を伝えることだ。

　ニューラルネットワークが実際には内部でどんなことをしているのかを、実際にPythonでゼロから実装しながら学んでいくことが主旨となる。
 
# 1. NNの構成要素：”ニューロン”
　まずは、NNの最も基本的な構成要素である、ニューロンについて。
　誤解を恐れずに極端に単純化すれば、ニューロンとは「**入力を受け取り、なんらかの計算をし、出力をする**」ものである。
 　以下は2つの値を入力値としてとるニューロンの図。
  ![図1.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3989773/e38d3c7e-1f8c-403b-979b-cdd6afbb87dd.jpeg)
ここでは3つの操作がなされている。

1．まず、それぞれの入力に重み（w）がかけられる。
```math
\displaylines{
x_1→x_1*w_1\\
x_2→x_2*w_2}
```
2．次に、重みがかけられた入力値とバイアス（b）を足し合わせる（最後に図の＋の部分）。
```math
x_1*w_1+x_2*w_2+b
```
3． 最後に、それらの和が活性化関数（AF：Activation Function）に渡される。
```math
y=f(x_1*w_1+x_2*w_2+b)
```

　活性化関数はまとまりがない入力を、扱いやすい形の出力に変えるために使われる。よく使われる活性化関数にはシグモイド関数やReLU関数（正規化線形関数）などがある。
 
 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3989773/0caafaff-c4a8-4a71-86af-6f903bda7eb1.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3989773/61da7243-2dc5-409f-9493-cd92dc074808.png)

※シグモイド関数は、(-∞, +∞)の入力を(0, 1)に圧縮するようなイメージ。

### 例 
　2つの入力値をとるニューロンがシグモイド関数を活性化関数として持っており、以下のパラメータを持っているとする。
```math
\displaylines{
\begin{align}
w&=[0,\:1]\\
b&=4
\end{align}
}
```
　w = [0, 1] は、w1 = 0, w2 = 1をベクトルの形で書き直したものである。ここにx = [2, 3] の入力を与える。ここではドット積（内積）を使用して、
```math
\displaylines{
\begin{align}
(w·x)+b&=w_1*x_1+w_2*x_2+b\\
&=0*2+1*3+4\\
&=7\\
\\
f(w·x+b)&=f(7)=0.999
\end{align}
}
```
　入力 x = [2, 3] に対するニューロンの出力は0.999となる。**NNの基本はこれだけ！**
　このような、入力を渡して出力を得るプロセスは、**フィードフォワード**（**feed forward**）と呼ばれる。
 
 ## Python実装①：ニューロン
  以下では、NNの内部でどんな計算が行われているかをより理解するために、Tensorflowなどは用いずNumPyを使って実装する。
```Python
import numpy as np

def sigmoid(x):
  # 活性化関数：f(x) = 1 / (1 + e^(-x))
  return 1 / (1 + np.exp(-x))

class Neuron:
  def __init__(self, weights, bias):
    self.weights = weights
    self.bias = bias

  def feedforward(self, inputs):
    # 入力に重みをかけ、バイアスを加え、活性化関数を使う
    total = np.dot(self.weights, inputs) + self.bias
    return sigmoid(total)

weights = np.array([0, 1]) # w1 = 0, w2 = 1
bias = 4                   # b = 4
n = Neuron(weights, bias)

x = np.array([2, 3])       # x1 = 2, x2 = 3
print(n.feedforward(x))    # 0.9990889488055994
```
最終行のコメントアウトしている部分、0.9990889488055994が、先ほど手計算して得た0.999とほぼ同値であることがわかる。

# 2. ニューロンを組みあわせる：”ネットワーク”
　ニューラルネットワークとは、つまるところたくさんのニューロンがつながりあっているものにすぎない。
 
![図２.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3989773/2ca2e86f-053d-4d15-8f7e-bd06150e10ee.jpeg)
- IL：入力層（Input Layer）
- HL：隠れ層/中間層（Hidden Layer）
- OL：出力層（Output Layer）

　このネットワークには、2つの入力値、2つのニューロンからなる隠れ層（h1, h2）、1つのニューロンからなる出力層（o1）がある。o1への入力値はh1とh2の出力値であり、このようなつながりにより、「ネットワーク」と呼ばれるのである。

　ここで、**隠れ層**は入力層と出力層の間のあらゆる層のことを言う。隠れ層は当然、複数あることもある。
### 例：フィードフォワード
　全てのニューロンが同じ重み w = [0, 1]、同じバイアス b = 0、同じ活性化関数（シグモイド関数）を持つとする。入力値を x = [2, 3] とすると、
```math
\displaylines{
\begin{align}
h_1=h_2&=f(w·x+b)\\
&=f(0*2+1*3+0)\\
&=f(3)\\
&=0.9526\\
\\
o_1&=f(w·[h_1,\:h_2]+b)\\
&=f(0*h_1+1*h_2+0)\\
&=f(0.9526)\\
&=0.7216
\end{align}
}
```
　このNNにおける入力値 x = [2, 3] に対する出力は0.7216となる。

　NNは隠れ層とニューロンをいくつでも持つことができる。しかし、層とニューロンが増えても、「入力値を出力層に向かって餌やり（feed）しながら前進（forward）し、出力を得る」という基本は変わらない。

## Python実装②：ニューラルネットワーク
```Python
import numpy as np

# 前の実装①から続く...

class NeuralNetwork:
  """
  - 2つの入力
  - 2つのニューロンからなる隠れ層
  - 1つのニューロンからなる出力層
  - それぞれのニューロンは同じ重みとバイアスを持つ
    - w = [0, 1]
    - b = 0
  """
  def __init__(self):
    wights = np.array([0, 1])
    bias = 0

    # 実装①のNeuronクラス
    self.h1 = Neuron(weights, bias)
    self.h2 = Neuron(weights, bias)
    self.o1 = Neuron(weights, bias)

  def feedforward(self, x):
    out_h1 = self.h1.feedforward(x)
    out_h2 = self.h2.feedforward(x)

    # o1への入力はh1, h2の出力
    out_o1 = self.o1.feedforward(np.array([out_h1, out_h2]))

    return out_o1
    
network = NeuralNetwork()
x = np.array([2, 3])
print(network.feedforward(x)) # 0.7216325609518421
```
手計算で得た0.7216と同値。

 # 3. ニューラルネットワークの学習Ⅰ
 　下記のデータセットを用意する。

| 名前 | 体重(kg) | 身長(cm) | 性別 |
| ---- | :-----: | :-----: | :-: |
|花子|60       |165      |女     |
|太郎|72       |182      |男     |
|裕太|69       |177      |男     |
|明美|54       |162      |女     |

　ニューラルネットワークに、身長と体重から性別を予測してもらおう。

　まず、データを扱いやすい形にする。身長と体重を標準化し、男性を0、女性を1のバイナリに変換する。
 
| 名前 | 体重(標準化) | 身長(標準化) | 性別 |
| ---- | :--------: | :--------: | :-: |
|花子|-3.75         |-6.5         |1   |
|太郎|8.25          |10.5         |0   |
|裕太|5.25          |5.5          |0   |
|明美|-9.75         |-9.5         |1   |

### 損失関数
　NNに学習させる前に、まずはどれくらい上手く学習できたかを計測できる必要がある。すなわち、あまり上手でない学習による「損失」を計測する関数を考える。
　ここでは、**平均二乗誤差**（MSE：Mean Square Error）を用いて損失を計算する。
```math
MSE=\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat y_i)^2
```
- n はサンプル数。ここでは4
- y は予測しようとしている変数。ここでは性別
- yi は実際の y の値。たとえばy1は花子で、1（女性）
- ^yi は推測された y の値（ネットワークの出力）

　(yi - ^yi)^2は**二乗誤差**と呼ばれる。損失関数は二乗の誤差の平均を計算する（平均二乗誤差）。予測が上手くいけばいくほど、損失関数は小さくなる。

> 予測が上手 = 損失が小さい
> **ネットワークによる学習 = 損失の最小化**

### 例：損失の計算
　たとえば、ネットワークが常に1（女性）を返すとする。損失は、
| 名前 | 予想 | 実際 | (実際 - 予想)^2 |
| ---- | :-: | :-: | :------------: |
|花子  |1    |1     |0              |
|太郎  |1    |0     |1              |
|裕太  |1    |0     |1              |
|明美  |1    |1     |0              |
```math
MSE=\frac{1}{4}(0+1+1+0)=0.5
```

## Python実装③：損失関数（MSE）
```Python
import numpy as np

def mse_loss(y_actual, y_pred):
  """
  - y_actual : 実際のyの値
  - y_pred   : 予測されたyの値
  y_actual, y_predは同じ長さを持つ配列
  """
  return((y_actual - y_pred)**2).mean()

y_actual = np.array([1, 0, 0, 1])
y_pred = np.array([1, 1, 1, 1])

print(mse_loss(y_actual, y_pred)) # 0.5
```

# 4. ニューラルネットワークの学習Ⅱ
　学習の目標が**損失の最小化**ということがわかれば、NNにおける重みとバイアスが予測に影響を与えているので、それらをいじればいいことがわかる。
　しかし、損失を減らせるような重みとバイアスをどのように設定すればいいのだろうか。
 
:::note warn
以下では結構複雑な他変数の微積分が出てきます。微積分が苦手な場合は、数学の部分は飛ばしてかまいません。
:::
　まずは単純化のために、花子のみのデータセットを想定する。

| 名前 | 体重(標準化) | 身長(標準化) | 性別 |
| ---- | :--------: | :--------: | :-: |
|花子|-3.75         |-6.5         |1   |

すると、平均二乗誤差は花子の二乗誤差となる
```math
\displaylines{
\begin{align}
MSE&=\frac{1}{1}\sum_{i=1}^{1}(y_i-\hat y_i)^2\\
&=(y_1-\hat y_1)^2\\
&=(1-\hat y_1)\\
\end{align}
}
```
　損失は重みとバイアスによる関数であるから、ネットワーク内のそれぞれの重みとバイアスにラベルをつける。
 
![図3.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3989773/a01a1a36-c839-461b-9542-210ea0e6b2a5.jpeg)
損失関数Lは、以下のような多変数関数になる。
```math
L(w_1,w_2,w_3,w_4, w_5,w_6,b_1,b_2,b_3)
```
w1の値を変えたとき、損失はどのように変化するだろうか。

　この問いを数式化すると以下のようになる。
```math
\frac{\partial L}{\partial w_1}
```
上記の偏導関数を次のように書きかえてみよう。
```math
\frac{\partial L}{\partial w_1}=\frac{\partial L}{\partial \hat y}*\frac{\partial \hat y}{\partial w_1}
```
右辺の掛け算の左側については、すでに L = (1 - ^yi)^2 で算出しているため、以下のようになる。
```math
\displaylines{
\begin{align}
\frac{\partial L}{\partial \hat y}&=\frac{\partial (1-\hat y)^2}{\partial \hat y}\\
&=2(1-\hat y)(-1)\\
&=-2(1-\hat y)
\end{align}
}
```
ここで、右辺の掛け算の右側についても考えてみよう。h1, h2, o1はニューロンの出力であるから、
```math
\displaylines{
\hat y = o1 = f(w_5h_1+w_6h_2+b_3)\\
fはシグモイド関数
}
```
w1はh1にしか影響しないため、式は以下のようになる。
```math
\displaylines{
\frac{\partial \hat y}{\partial w_1}
=\frac{\partial \hat y}{\partial h_1}*\frac{\partial h_1}{\partial w_1}\\
\frac{\partial \hat y}{\partial h_1}
=w_5*f'(w_5h_1+w_6h_2+b_3)\\
}
```
また、
```math
\displaylines{
h_1=f(w_1x_1+w_2x_2+b_1)\\
\frac{\partial h_1}{\partial w_1}=x_1*f'(w_1x_1+ w_2x_2 + b_1)
}
```
ここでは、x1は体重、x2を身長とする。シグモイド関数f(x)を微分すると、
```math
\displaylines{
f(x)=\frac{1}{1+e^{-x}}\\
f'(x)=\frac{e^{-x}}{(1+e^{-x})}=f(x)*(1-f(x))
}
```
この式には後で立ち戻る。

　さて、当初の問題「w1の値を変えたとき、損失はどのように変化するか」について、計算可能な式に置き換えた。
```math
\frac{\partial L}{\partial w_1}=\frac{\partial L}{\partial \hat y_1} * \frac{\partial \hat y_1}{\partial h_1}*\frac{\partial h_1}{\partial w_1}
```
　このように、偏微分計算を後ろから前にしてく方法を「**誤差逆伝播法」（backpropagation**）と呼び、NNにおいて非常に重要な考え方である。

### 例：偏微分計算
　とはいえ、記号だらけでは一体何をしているのかわかりづらいだろう。ということで、花子の例に再び戻る。
 
| 名前 | 体重(標準化) | 身長(標準化) | 性別 |
| ---- | :--------: | :--------: | :-: |
|花子|-3.75         |-6.5         |1   |

　まずは全ての重みを1、バイアスを0とし、このデータについてfeedforwardで学習させよう。
 ```math
\displaylines{
\begin{align}
h_1&=f(w_1x_1+w_2x_2+b_1)\\
&=f(-2+-1+0)\\
&=0.0474
\end{align}
}
```
```math
\displaylines{
\begin{align}
h_2&=f(w_3x_1+w_4x_2+b_2)=0.0474
\end{align}
}
```
```math
\displaylines{
\begin{align}
o_1&=f(w_5h_1+w_6h_2+b_3)\\
&=f(0.0474+0.0474+0)\\
&=0.524
\end{align}
}
```
このNNの出力は0.524と、男（0）でも女（1）でもなさそうな数値だった。次に、w1の変化に対するLの変化を見てみよう。
```math
\displaylines{
\frac{\partial L}{\partial w_1}=\frac{\partial L}{\partial \hat y_1} * \frac{\partial \hat y_1}{\partial h_1}*\frac{\partial h_1}{\partial w_1}
}
```
```math
\displaylines{
\begin{align}
\frac{\partial L}{\partial \hat y_1}&=-2(1-\hat y_1)\\
&=-2(1-0.524)\\
&=-0.952
\end{align}
}
```
```math
\displaylines{
\begin{align}
\frac{\partial \hat y_1}{\partial h_1}&=w_5*f'(w_5h_1+w_6h_2+b_3)\\
&=1*f'(0.0474+0.0474+0)\\
&=f(0.0948)*(1-f(0.0948))\\
&=0.249
\end{align}
}
```
```math
\displaylines{
\begin{align}
\frac{\partial h_1}{\partial w_1}&=x_1*f'(w_1x_1+ w_2x_2 + b_1)\\
&=1*f'(0.0474+0.0474+0)\\
&=-2*f'(-2+-1+0)\\
&=-2*f(-3)*(1-f(-3))\\
&=-0.0904
\end{align}
}
```
```math
\displaylines{
\begin{align}
\frac{\partial L}{\partial w_1}&=-0.952*0.249*-0.0904\\
&=0.0214
\end{align}
}
```
この式からわかるのは、w1を増やすと、Lは「ごくわずかに増加する」ことがわかった。

### 確率的勾配降下法
　NNを学習させるツールが全て揃ったところで、**確率的勾配降下法**（**stochastic gradient descent**）と呼ばれる最適化のためのアルゴリズムを用いて、損失が最小となる重みとバイアスを求めよう。確率的勾配降下法は、以下のように示される。
 ```math
 w_1←w_1-\eta \frac{\partial L}{\partial w_1}
 ```
　ηは**学習率**という定数で、学習の速さを決める。ここでは、w1からη倍したw1に対するLの変動を引いているだけである。
- もし変動が正であるなら、w1は減少する。そのため、Lは減少する
- もし変動が負であるなら、w1は増加する。そのため、Lは増加する

　NN中の全ての重みとバイアスに関してこれらの操作を行えば、損失は徐々に減少し、NNの学習は改善していく、というわけだ。

 　ここまでの学習プロセスをまとめると、
  1. データセットから1つのサンプルを取り出す（これにより確率的勾配降下法が可能になる）
  1. 損失に関する全ての偏微分を計算する
  1. 重みとバイアスを少しずつ変動させる
  1. 1に戻る
  
  これを繰り返していくことで、NNはより「良く」学習する。

## Python実装④：ニューラルネットワーク完成版
| 名前 | 体重(標準化) | 身長(標準化) | 性別 |
| ---- | :--------: | :--------: | :-: |
|花子|-3.75         |-6.5         |1   |
|太郎|8.25          |10.5         |0   |
|裕太|5.25          |5.5          |0   |
|明美|-9.75         |-9.5         |1   |

![図3.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3989773/a01a1a36-c839-461b-9542-210ea0e6b2a5.jpeg)

```Python
import numpy as np

def sigmoid(x):
  # 活性化関数（シグモイド関数）：f(x) = 1 / (1 + e^(-x))
  return 1 / (1 + np.exp(-x))

def deriv_sigmoid(x):
  # シグモイド関数の微分：f'(x) = f(x) * (1 - f(x))
  fx = sigmoid(x)
  return fx * (1 - fx)

def mse_loss(y_actual, y_pred):
  # 実装③参照
  return ((y_actual - y_pred)**2).mean()

class NeuralNetwork:
  """
  以下のコードは学習目的のため、あくまで単純化したものです。
  実際のNNはもっと複雑であろうことは想像に難くありません。
  """
  def __init__(self):
    # 重み
    self.w1 = np.random.normal()
    self.w2 = np.random.normal()
    self.w3 = np.random.normal()
    self.w4 = np.random.normal()
    self.w5 = np.random.normal()
    self.w6 = np.random.normal()

    # バイアス
    self.b1 = np.random.normal()
    self.b2 = np.random.normal()
    self.b3 = np.random.normal()

  def feedforward(self, x):
    h1 = sigmoid(self.w1 * x[0] + self.w2 * x[1] + self.b1)
    h2 = sigmoid(self.w3 * x[0] + self.w4 * x[1] + self.b2)
    o1 = sigmoid(self.w5 * h1 + self.w6 * h2 + self.b3)
    return o1

  def train(self, data, all_y_actuals):
    """
    - (n x 2)のNumPy配列で、nはデータセットのサンプル数
    - all_y_actualsはn個の要素を持つNumPy配列
    all_y_actualsの要素はデータセットの要素に対応する。
    """
    learn_rate = 0.1
    epochs = 1000 # データセット内をループする回数

    for epoch in range(epochs):
      for x, y_actual in zip(data, all_y_actuals):
        # feedforward：後でこれらの値を利用
        sum_h1 = self.w1 * x[0] + self.w2 * x[1] + self.b1
        h1 = sigmoid(sum_h1)

        sum_h2 = self.w3 * x[0] + self.w4 * x[1] + self.b2
        h2 = sigmoid(sum_h2)

        sum_o1 = self.w5 * h1 + self.w6 * h2 + self.b3 
        o1 = sigmoid(sum_o1)
        y_pred = o1

        # 偏微分を計算
        # 命名規則：d_L_d_w1はw1の変動に対するLに変動を示す
        d_L_d_ypred = -2 * (y_actual - y_pred)

        # ニューロンo1
        d_ypred_d_w5 = h1 * deriv_sigmoid(sum_o1)
        d_ypred_d_w6 = h2 * deriv_sigmoid(sum_o1)
        d_ypred_d_b3 = deriv_sigmoid(sum_o1)

        d_ypred_d_h1 = self.w5 * deriv_sigmoid(sum_o1)
        d_ypred_d_h2 = self.w6 * deriv_sigmoid(sum_o1)

        # ニューロンh1
        d_h1_d_w1 = x[0] * deriv_sigmoid(sum_h1)
        d_h1_d_w2 = x[1] * deriv_sigmoid(sum_h1)
        d_h1_d_b1 = deriv_sigmoid(sum_h1)

        # ニューロンh2
        d_h2_d_w3 = x[0] * deriv_sigmoid(sum_h2)
        d_h2_d_w4 = x[1] * deriv_sigmoid(sum_h2)
        d_h2_d_b2 = deriv_sigmoid(sum_h2)

        # 重みとバイアスを更新
        # ニューロンh1
        self.w1 -= learn_rate * d_L_d_ypred * d_ypred_d_h1 * d_h1_d_w1
        self.w2 -= learn_rate * d_L_d_ypred * d_ypred_d_h1 * d_h1_d_w2
        self.b1 -= learn_rate * d_L_d_ypred * d_ypred_d_h1 * d_h1_d_b1

        # ニューロンh2
        self.w3 -= learn_rate * d_L_d_ypred * d_ypred_d_h2 * d_h2_d_w3
        self.w4 -= learn_rate * d_L_d_ypred * d_ypred_d_h2 * d_h2_d_w4
        self.b2 -= learn_rate * d_L_d_ypred * d_ypred_d_h2 * d_h2_d_b2

        # ニューロンo1
        self.w5 -= learn_rate * d_L_d_ypred * d_ypred_d_w5
        self.w6 -= learn_rate * d_L_d_ypred * d_ypred_d_w6
        self.b3 -= learn_rate * d_L_d_ypred * d_ypred_d_b3
        
      # 各ループ（epoch）ごとに損失を計算する
      if epoch % 10 == 0:
        y_preds = np.apply_along_axis(self.feedforward, 1, data)
        loss = mse_loss(all_y_actuals, y_preds)
        print("Epoch %d loss: %.3f" % (epoch, loss))

# データセットを定義
data = np.array([
  [-3.75, -6.5], # 花子
  [8.25, 10.5],  # 太郎
  [5.25, 5.5],   # 裕太
  [-9.75, -6],   # 明美
])
all_y_actuals = np.array([
  1, # 花子
  0, # 太郎
  0, # 裕太
  1, # 明美
])

# NNに学習させる
network = NeuralNetwork()
network.train(data, all_y_actuals)
```

出力された数値を見ると、損失が徐々に減少しているのがわかる。

　このモデルを利用して、身長と体重から性別を推測することができる。
```Python
akira = np.array([-12.75, -4.5]) # 51kg, 167cm
yuki = np.array([12.25, 3.5]) # 76kg, 175cm
print("アキラ: %.3f" % network.feedwork(akira)) # 0.962
print("ユウキ: %.3f" % network.feedwork(yuki))  # 0.038
```
# まとめ
これまでの学習内容を振り返る。
- ニューロンはニューラルネットワークの構成要素である
- 活性化関数（シグモイド関数）
- ニューラルネットワークとは、単にニューロンがつながりあっているだけ
- 体重と身長という特徴量を入力とし、性別を出力としたデータセットを作成
- 損失関数と平均二乗誤差
- 「学習」とは「損失の最小化」
- 誤差逆伝播法
- 確率的勾配降下法

拙文で（おそらく）最も基本的なニューラルネットワーク内のアルゴリズムについて、少しでも理解の助けになれば幸いです。