Pythonを使用したNNの実装


NNの学習の一環として学んだことの備忘録

## 問題設定
画像のような「色のついた部分を1、色のついていない部分を0」として予測ができるようなNNを実装することが目的

関数的に表すと
(x,y)→[NN]→0か1
入力→[NN]→出力
![正解.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/8b61e51d-a6c1-5617-41a0-5ae449d02842.png)

## NN以外で使用する関数
二次元配列を入力としてヒートマップを出力する関数
```python
def draw_heatmap(data,n):#作成したデータを可視化するために2次元配列のヒートマップを出力する関数
    fig, ax = plt.subplots()
    ax.pcolor(data, cmap="Reds")
    ax.set_title("self_made")
    ax.set_xticks(np.array([n*0.25,n*0.75]))
    ax.set_yticks(np.array([n*0.25,n*0.75]))
    ax.invert_yaxis()
    ax.xaxis.tick_top()
    #正解の線の描画
    ax.grid(color="y",linestyle="dotted",linewidth=1)
    plt.show()
```
## 理想的なNN
どのようになれば最小のパラメータで目的の関数と同じ働きをするNNを構築できるかを考えることは実験を行う上で重要である。

活性化関数：ステップ関数（入力xが0以下では0、0以上で1を返す）
中間層：4

アーキテクチャは以下の通りである。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/5e74b13e-088e-1647-d8d0-583e7a62997c.png)

中間層のそれぞれで
xが25％以上で1を
xが75％以下で1を
yが25％以上で1を
yが75％以下で1を
を出力しすべての中間層が1を出力した場合に出力層で1を出力するように実装している。

実装したコードは以下の通りである。
```python
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm#かっこいいから採用

size_num = 1000#要素の数は二乗
persentage = 1.0
data_size = (size_num,size_num)#データの大きさを指定（タプルで指定）
ideal_list = np.zeros(data_size)
np.random.seed(0)

# --入力層--

class InputLayer:    
    def forward(self,input_data):
        self.y = input_data
        
# -- 中間層 --

class MiddleLayer:
    def __init__(self):  # 初期設定
        self.w = np.array([[0.5,0],[-0.5,0],[0,0.5],[0,-0.5]])  # 重み（行列）
        self.b = np.array([-0.5*0.25*size_num,0.5*0.75*size_num,-0.5*0.25*size_num,0.5*0.75*size_num])  # バイアス（ベクトル）

    def forward(self, x):  # 順伝播
        self.x = x
        u = np.dot(self.w, x) + self.b
        rate = 10000#これで固定
        self.y = 1/(1+np.exp(-1*rate*u))  # シグモイド関数

# -- 出力層 --
class OutputLayer:
    def __init__(self):  # 初期設定
        self.w = np.array([0.2,0.2,0.2,0.2])  # 重み（行列）
        self.b = np.array([-0.7])  # バイアス（ベクトル）
        
        self.nu_w = 0#Adamに使用
        self.s_w = 0
        self.nu_b = 0
        self.s_b = 0
    
    def forward(self, x):  # 順伝播
        self.x = x
        u = np.dot(self.x, self.w) + self.b
        rate = 10000#これで固定
        self.y = 1/(1+np.exp(-1*rate*u))  # シグモイド関数
        
        
# -- インスタンスの生成 --
input_layer = InputLayer()
middle_layer = MiddleLayer()
output_layer = OutputLayer()

def predict(x,y):
    input_layer.forward(np.array([x,y]))
    middle_layer.forward(input_layer.y)
    output_layer.forward(middle_layer.y)
    return output_layer.y
predict_data = np.zeros(data_size)
for i in tqdm(range(size_num)):
    for j in range(size_num):
        predict_data[i][j] = predict(i,j)
```
このNNモデルでは実行後predict(x,y)と入力することで予測ができるようになっている。
（ただしx,yの定義域は0以上size_num以下である。）

ヒートマップとして表すと視覚的には同じほど近似できていることがわかる。
![trivial.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/9a8df108-2ecc-621d-38e8-c8a2bc84f819.png)


## sk-learnで実装
### NNの設定
活性化関数：ReLU（入力xが0以下では0、0以上でxを返す）
最適化手法：Adam（Momentum × RMSProp）
中間層：それぞれ出力結果のヒートマップ上部に記載

sk-learnで実装で実装したコードは以下の通りである。
```python
import numpy as np
import matplotlib.pyplot as plt
import tqdm#かっこいいから採用

size_num = 500#要素の数は二乗
persentage = 0.3
data_size = (size_num,size_num)#データの大きさを指定（タプルで指定）
ideal_list = np.zeros(data_size)
np.random.seed(0)

def test_data(n,x,y):
    data = 0
    if x>=0.25*n and x<0.75*n and y>=0.25*n and y<0.75*n:
     data =1
    return data
    

sample_list = []#listの中にタプル（x,y,data）で格納
for _ in tqdm.tqdm(range(int(size_num*persentage))):
    randx = int(np.random.rand(1)*size_num)
    randy = int(np.random.rand(1)*size_num)
    sample_list.append(np.array([test_data(size_num,randx,randy),randx,randy]))


from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn import neural_network
from sklearn.metrics import accuracy_score, precision_score, recall_score

data = np.array(sample_list)
labels = data[:, 0:1] # 目的変数を取り出す
features = preprocessing.minmax_scale(data[:, 1:]) # 説明変数を取り出した上でスケーリング
x_train, x_test, y_train, y_test = train_test_split(features, labels.ravel(), test_size=0.3) # トレーニングデータとテストデータに分割



hidden_layers=(40,40)
clf = neural_network.MLPClassifier(activation="relu",
                                   solver='adam',
                                   alpha=0.001,
                                   hidden_layer_sizes=hidden_layers,
                                   max_iter=10000,
                                   ) 
clf.fit(x_train, y_train)

predict = clf.predict(x_test)
loss = clf.loss_curve_
loss_x = list(range(len(loss)))
plt.plot(loss_x,loss)
print()  
print(accuracy_score(y_test, predict),#正解率（accuracy）
      precision_score(y_test, predict, zero_division = 'warn'),#適合率（precision）
      recall_score(y_test, predict,zero_division = 'warn'))#再現率（recall）


print()
#print(clf.get_params)#条件の表示
#print(clf.coefs_)#重みの表示
#print(clf.intercepts_)#バイアスの表示

def predict(x,y):
    x_ = x/size_num
    y_ = y/size_num
    return clf.predict(np.array([[x_,y_]]))

predict_data = np.zeros(data_size)
for i in range(size_num):
    for j in range(size_num):
        predict_data[i][j] = predict(i,j)
        
def draw_heatmap(data,n):#作成したデータを可視化するために2次元配列のヒートマップを出力する関数
    fig, ax = plt.subplots()
    ax.pcolor(data, cmap="Reds")
    ax.set_title(f"hidden_layer={hidden_layers}")
    ax.set_xticks(np.array([n*0.25,n*0.75]))
    ax.set_yticks(np.array([n*0.25,n*0.75]))
    ax.invert_yaxis()
    ax.xaxis.tick_top()
    #サンプル点の描画
    xtrain = x_train * size_num 
    xtest = x_test * size_num
    ax.scatter(xtrain[:,0],xtrain[:,1],c="black")
    ax.scatter(xtest[:,0],xtest[:,1],c="blue")
    #正解の線の描画
    ax.grid(color="y",linestyle="dotted",linewidth=1)
    plt.show()

draw_heatmap(predict_data,size_num)
```
与えるデータ数を150とした場合の結果は以下の通りである。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/ced0d24f-f8c9-6f0d-41aa-4a1163fe2a8e.png)
黒い点が学習に用いたデータで、青い点がテストに使用したデータである。
この場合、テストでの正解数は100％であったが可視化すると正しく予測できていないことがわかる。

さらに与えるデータを増やし、50000点のデータを与えた場合の結果が以下である。
![predict_40_40_50000data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/81f18d81-daef-0cf6-0090-0e9312a3644f.png)
こちらもテストでの正解数は100％だが、先ほどのデータ数150とは異なり、可視化しても正しく予測できていることがわかる。

・このモデルの問題点は先ほど実装したNNモデルに比べて隠れ層の数が多く、計算に時間がかかってしまうこと。
隠れ層が多いために学習に必要なデータの数が多いこと
がある。

## 自作のNNで実装
この問題を解決するためには理想的なNNモデルと同様にステップ関数を使用し、0か1の出力にしなければならない。

しかしsk-learnにはステップ関数は微分不可能な為、活性化関数として用意されていない。

そこで自分でＮＮモデルを実装することにした。
自分で実装したからといってステップ関数が微分不可能であることは解決していないから微分可能なステップ関数を用意しなければならない。

そこでシグモイド関数をカスタマイズしてステップ関数に近似する。
使用するカスタマイズしたシグモイド関数のグラフは以下の通りである。
![sigmoid_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/65774d5a-e215-9303-8fb3-fdd58011e81f.png)
シグモイド関数1/(1+exp(-rate*x))としている。
(rate=1が通常のシグモイド関数)

NNモデルの設定
活性化関数：カスタマイズしたシグモイド関数
最適化手法：Adam
中間層：4

実装したコードが以下の通りである。
```python
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm#かっこいいから採用
import warnings#カスタマイズシグモイド関数の発散の警告を非表示にする


# -- 各設定値 --

wb_width = 0.1  # 重みとバイアスの広がり具合
eta = 0.001  # 学習係数
size_num = 500#要素の数は二乗
persentage = 1000
data_size = (size_num,size_num)#データの大きさを指定（タプルで指定）
ideal_list = np.zeros(data_size)
np.random.seed(0)
warnings.simplefilter('ignore')#警告を非表示

"""
def test_data(n,x,y):#ひし形
    data = 0
    boo1 = y+x >= 0.5*size_num
    boo2 = y-x >= -0.5*size_num 
    boo3 = y-x <= 0.5*size_num
    boo4 = y+x <= 1.5*size_num
    if boo1 and boo2 and boo3 and boo4:
     data =1
    return data
"""
def test_data(n,x,y):#正方形
    data = 0
    if x>=0.25*n and x<0.75*n and y>=0.25*n and y<0.75*n:
     data = 1
    return data
#"""
    
sample_list = []#listの中にタプル（x,y,data）で格納
for _ in tqdm(range(int(size_num*persentage))):
    randx = int(np.random.rand(1)*size_num)
    randy = int(np.random.rand(1)*size_num)
    sample_list.append(np.array([test_data(size_num,randx,randy),randx,randy]))



from sklearn import preprocessing
from sklearn.model_selection import train_test_split

data = np.array(sample_list)
labels = data[:, 0:1] # 目的変数を取り出す
features = preprocessing.minmax_scale(data[:, 1:]) # 説明変数を取り出した上でスケーリング
x_train, x_test, y_train, y_test = train_test_split(features, labels.ravel(), test_size=0.3) # トレーニングデータとテストデータに分割


# -- 入力層 --
class InputLayer:    
    def forward(self,input_data):
        self.y = input_data
        
# -- 中間層 --

class MiddleLayer:
    def __init__(self):  # 初期設定
        self.w = wb_width * np.random.randn(4, 2)  # 重み（行列）
        self.b = wb_width * np.random.randn(4)  # バイアス（ベクトル）

        #Adamに使用
        self.nu_w = 0
        self.s_w = 0
        self.nu_b = 0
        self.s_b = 0
        self.q = 0#更新回数

    def forward(self, x):  # 順伝播
        self.x = x
        u = np.dot(self.w, x) + self.b
        rate = 10000#これで固定
        self.y = 1/(1+np.exp(-1*rate*u))  # シグモイド関数

    def backward(self, grad_y):  # 逆伝播
        beta1 = 0.9
        beta2 = 0.999
        rate = 10000
        self.q += 1
        delta = rate * grad_y * (1-self.y)*self.y  # シグモイド関数の微分
        delta = delta.reshape(1,-1)
        self.x = self.x.reshape(1,-1)
        
        self.grad_w = np.dot(delta.T,self.x)
        self.grad_b = np.sum(delta, axis=0)
        self.grad_x = np.dot(delta, self.w)
        
        self.nu_w = beta1*self.nu_w + (1-beta1)*self.grad_w
        self.nu_w = (1/(1-(beta1**self.q))) * self.nu_w
        self.s_w = beta2*self.s_w + (1-beta2)*((self.grad_w)**2)
        self.s_w = (1/(1-(beta2**self.q))) * self.s_w
        
        self.nu_b = beta1*self.nu_b + (1-beta1)*self.grad_b
        self.nu_b = (1/(1-(beta1**self.q))) * self.nu_b
        self.s_b = beta2*self.s_b + (1-beta2)*((self.grad_b)**2)
        self.s_b = (1/(1-(beta2**self.q))) * self.s_b

    def update(self, eta):  # 重みとバイアスの更新
        ep = 1e-8
        adam_w = self.nu_w/np.sqrt(self.s_w + ep)
        adam_b = self.nu_b/np.sqrt(self.s_b + ep)  
        self.w -= eta * adam_w
        self.b -= eta * adam_b


# -- 出力層 --
class OutputLayer:
    def __init__(self):  # 初期設定
        self.w = wb_width * np.random.randn(4)  # 重み（行列）
        self.b = wb_width * np.random.randn(1)  # バイアス（ベクトル）

        #Adamに使用        
        self.nu_w = 0
        self.s_w = 0
        self.nu_b = 0
        self.s_b = 0
        self.q = 0
    
    def forward(self, x):  # 順伝播
        self.x = x
        u = np.dot(self.x, self.w) + self.b
        rate = 10000#これで固定
        self.y = 1/(1+np.exp(-1*rate*u))  # シグモイド関数
    
    def backward(self, t):  # 逆伝播
        beta1 = 0.9
        beta2 = 0.999
        rate = 10000
        self.q += 1
        self.grad_w = rate * self.x * (self.y - t)
        self.grad_x = rate * self.w * (self.y - t)
        self.grad_b = rate * (self.y - t)
        
        self.nu_w = beta1*self.nu_w + (1-beta1)*self.grad_w
        self.nu_w = (1/(1-(beta1**self.q))) * self.nu_w
        self.s_w = beta2*self.s_w + (1-beta2)*((self.grad_w)**2)
        self.s_w = (1/(1-(beta2**self.q))) * self.s_w
        
        self.nu_b = beta1*self.nu_b + (1-beta1)*self.grad_b
        self.nu_b = (1/(1-(beta1**self.q))) * self.nu_b
        self.s_b = beta2*self.s_b + (1-beta2)*((self.grad_b)**2)
        self.s_b = (1/(1-(beta2**self.q))) * self.s_b

    def update(self, eta):  # 重みとバイアスの更新
        ep = 1e-8
        adam_w = self.nu_w/np.sqrt(self.s_w + ep)
        adam_b = self.nu_b/np.sqrt(self.s_b + ep)  
        self.w -= eta * adam_w
        self.b -= eta * adam_b
        
        
# -- インスタンスの生成 --
input_layer = InputLayer()
middle_layer = MiddleLayer()
output_layer = OutputLayer()

for i in tqdm(range(len(x_train))):
    # -- 順伝播 -- 
    input_layer.forward(x_train[i])
    middle_layer.forward(input_layer.y)
    output_layer.forward(middle_layer.y)

    """#伝播の様子の確認
    print(input_layer.y)#入力
    print(middle_layer.y)#中間層
    print(output_layer.y)#出力
    """
    # -- 逆伝播 --
    output_layer.backward(y_train[i])
    output_layer.update(eta)
    middle_layer.backward(output_layer.grad_x)
    middle_layer.update(eta)


correct_num = 0
accurate_list = []
for i in range(len(x_test)):
    input_layer.forward(x_test[i])
    middle_layer.forward(input_layer.y)
    output_layer.forward(middle_layer.y)
    Y = int(output_layer.y)
    if Y == y_test[i]:
        correct_num += 1
    accurate = abs(Y - y_test[i])/(y_test[i]+1e-4)
    accurate_list.append(accurate)


def predict(x,y):
    x_ = x/size_num
    y_ = y/size_num
    input_ = np.array([x_,y_])
    input_layer.forward(input_)
    middle_layer.forward(input_layer.y)
    output_layer.forward(middle_layer.y)
    return output_layer.y

predict_data = np.zeros(data_size)
for i in tqdm(range(size_num)):
    for j in range(size_num):
        predict_data[i][j] = predict(i,j)

        
def draw_heatmap(data,n):#作成したデータを可視化するために2次元配列のヒートマップを出力する関数
    fig, ax = plt.subplots()
    ax.pcolor(data, cmap="Reds")
    ax.set_title("self_made")
    ax.set_xticks(np.array([n*0.25,n*0.75]))
    ax.set_yticks(np.array([n*0.25,n*0.75]))
    ax.invert_yaxis()
    ax.xaxis.tick_top()
    """
    #サンプル点の描画
    xtrain = x_train * size_num 
    xtest = x_test * size_num
    ax.scatter(xtrain[:,0],xtrain[:,1],c="black")
    ax.scatter(xtest[:,0],xtest[:,1],c="blue")
    """
    #正解の線の描画
    ax.grid(color="y",linestyle="dotted",linewidth=1)
    plt.show()

draw_heatmap(predict_data,size_num)
```

結果が以下の通りである。
![self_made4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2780803/cca337e3-9cda-6632-39a9-0c7b659677fb.png)

視覚的に評価するとうまく学習できたとはいえないだろう。

## 総括
中間層を増やすことでデータ数が5000程度あれば有効な近似をすることができた。
しかし、中間層が少ないとデータ数が10000あっても有効な近似が出来なかった。

なぜこのようなことになったのかとして考察すると、
・更新量に微分項があり、ステップ関数に近似したシグモイド関数を使用するとx=0の部分での変化量が大きすぎるため微分可能であっても学習ができない
・人間が150個のデータ点を見て予測できるのかと考えると縁の部分は正確に予測をすることは難しい
・更新量を確認すると最大で1e-118の更新量しかなく求める値にたどり着くまでに0か1かしか値をとらない問題には膨大なデータが必要
ということが考えられる。

特に2つ目の理由に関して、重みの初期値はランダムな0～1の値をとっており、理想的なNNでは半分の重みの値が0である。
学習をし、重みの値を0にすることは重み＝誤差となる必要があり、データ数不足やオーバシュートの観点から完全に一致することは難しい。

学習をしてうまく出来た部分ではデータ数と中間層を増やし、表現力を向上させることで有効な近似ができたのではないかと考える。