Python
DeepLearning
bioinformatics
Chainer
morriswatermazetest

Chainerを使って素人が水迷路試験を4日で終わらせる方法

More than 1 year has passed since last update.

モリス水迷路試験とは

 モリス水迷路試験は1981年にRichard G. Morrisらが考案した行動試験で、主に齧歯類の空間学習能力や記憶力を評価するために用いられます。具体的には、水を貯めた円形のプールにマウスを放ち、固定された足場までの到達時間を測定します。足場は水面下に隠れているのでマウスからは見えませんが、施行を繰り返すことで周囲の景色を手掛かりとしてゴールの場所を特定できるようになります。すなわち、学習が進めばゴールへの到達時間が短縮していくことになるため、認知機能の指標になるのです。トレーニングの日数は厳密には決まっていませんが、Tamara Timicらの論文(Behav Brain Res. 2013 15;241:198-205)にあるように、5日間が標準的なようです。
mwm_ex.png

人生で最も退屈な5日間

 実験操作をしていると分かるのですが、4日間順調に学習が進んだマウスは、当然5日目にも良好な成績を収めます。逆に、毎回100秒以上時間が掛かっていたマウスが、5日目に突然短時間でゴールするという事態は稀です。つまり、ほとんど予想通りの結果を確認するために、我々は単純作業を繰り返さなくてはいけないのです。そのことに気付いた時、私のゴーストが囁きました。
「退屈なことはPythonにやらせよう」

Chainerで構築する多層パーセプトロン(MLP)

 退屈な工程を機械に肩代わりさせるため、自分の脳ミソではなくArtificial Neural Network(ANN)にデータを放り込むことに決めました。
 この記事では、深層学習初心者の私が「Chainerで水迷路結果予測モデルを構築」する過程を振り返りながら、知識の整理をしていくことを目的としています。Chainerをpip installする方法はこちらの記事が非常に分かりやすかったです。

 さて今回は「4日分の入力から最終日の値を計算する」という数理モデルを得たいので、4つの入力ノードと1つの出力ノードからなる多層パーセプトロンを構築することにしました。
ANN_figure.png
 Neural Networkは生物の神経回路を模して作られたものであり、上図のノードはニューロン、リンクはシナプスに相当します。x1からx4の各ノードに水迷路試験のday1からday4の到達時間を入力して出力(y)を得るわけですが、これがday5の到達時間に近い値になれば成功です。
 ちなみに今回のように隠れ層(Hidden Layer)を2層以上使って機械学習をする場合、深層学習(Deep Learning)と呼んで良いようです。MLPのクラス定義は下記のようになりました。

import numpy as np
import chainer.functions as F
import chainer.links as L
from chainer import Variable, Chain, optimizers

class MLP(Chain):
    def __init__(self, n_units, n_output):
        super(MLP, self).__init__(
            l1 =L.Linear(None, n_units),
            l2 =L.Linear(None, n_units),
            l3 =L.Linear(None, n_output)
        )

    def __call__(self, x):
        h1 = F.leaky_relu(self.l1(x))
        h2 = self.l2(h1)
        return self.l3(h2)

 まどろっこしい書き方ですが、慣例に則りChainクラスを継承した上で、スーパークラスのコンストラクタを借りて__init__を定義します。仮引数をn_units, n_outputとすることで、インスタンス生成時に入出力ノードの数を指定できるようにしました。
 入力層から出力層までの4層は3つの線形結合(links)で繋がれます。具体的には、l1,l2,l3のそれぞれにおいて「入力値に重み(w)をかけてバイアス(b)を足したもの」が次のノードに渡されるわけです。なおchainer.links.Linearの第一引数にNoneを与えると自動的に入力値の個数に対応してくれるので便利です。

 続いて__call__メソッドを定義します。ここでは引数(x)を与えてMLPインスタンスを呼び出した時に、最終的に出力される値を定義する、すなわち順伝播(forward propagation)時の処理を書きます。入力層にxが渡されると、リンク1(l1)で線形に処理されて隠れ層1(h1)に到達します。h1では次のノードに伝える前に活性化関数(F.leaky_relu())で処理され、リンク2(l2)で再び線形に処理されて隠れ層2(h2)に辿り着きます。h2では活性化関数を使わず、そのままリンク3(l3)で線形に処理して返り値となります(上図のy)。

 隠れ層1で活性化関数を使うのは非線型のモデルを構築するためです。線形結合だけのモデル(線形分類器)では学習精度が上がらないことが分かっており、入出力処理にパンチを効かせるために活性化関数を使います。今回は負数の入力のみ0.01倍に変換するLeaky ReLUを使ってみました(普通のReLUだと負数が全て0になってしまいますからね)。

model = MLP(4,1) #4入力1出力のMLPインスタンスを作成
optimizer = optimizers.Adam() #Adaptive Moment Estimationを採用
optimizer.setup(model)

 クラス定義ができたら早速インスタンスを作りましょう。4入力1出力のMLPインスタンスを作成し、modelと名付けます。modelは生まれたばかりのヘナチョコMLPなので、各線形結合の持つ重み(w)はランダムな初期値になっています。今の時点でy = model(x)とcallしても、謎の数値が返ってくるばかりです。学習を繰り返しながらlinksのparameterをうまく調整して、正解との誤差を少なくしていく最適化作業が必要になります。

 出力層における誤差を微分すると一つ上の層の誤差が分かるので、層を遡りながら勾配を計算していくのが誤差逆伝播法(backward propagation, backpropagation)です。こうして得られた勾配から、最終的な誤差が最小になるようなparameterを求める手法が確率的勾配降下法(SGD)であり、optimizersモジュールに定義されています。今回はSGDの一種であるAdamを採用することとし、setup()メソッドを使ってmodelをセットしました。Backpropagationの処理は後ほどトレーニングを行うブロックに書きます。

学習用データを準備する

 深層学習のモデルが出来上がったので、今度は学習データ(教師データ)を準備しましょう。準備といっても簡単で、以下の3つの作業を順番に行うだけです。
 ①インプット用とアウトプット用のリストに分ける
 ②それらをNumpyの多次元配列に変換する
 ③最後にChainerで扱えるVariable変数に代入する

xls_ex.png
 まずは上図のエクセルデータから通常のPythonリストにデータを振り分けます。エクセルからの読み込みには巷で評判の良いxlrdモジュールを使用しました。

import xlrd

book = xlrd.open_workbook("data_set.xlsx")
sheet = book.sheet_by_index(0) #sheet1はindex(0)
N = sheet.nrows -1 #ヘッダを除いた行数がマウスの総数

raw_4ds = [] #day1-day4をinput用リストに入れる
raw_5th = [] #day5をoutputリストに入れる
for row in range(1,N+1):
    _4day = [(sheet.cell(row, col).value)/120 for col in range(2,6)]
    raw_4ds.append(_4day)
    raw_5th.append([(sheet.cell(row, 6).value/120)])

 入力値リスト(raw_4ds)にはday1からday4の結果をマウス毎にまとめて二次元配列として格納します。出力値リスト(raw_5th)はday5の値だけなので本来は一次元(ベクトル)なのですが、Chainerでの処理のため要素数1の二次元配列の形にします。
 ここで数値を120で割っているのは、入力値を0-1の範囲に収めるためです。水迷路試験はマウスを最大120秒間泳がせるのですが、入力値が0-1の方が正確に計算できるという都市伝説に準拠しました。

x = Variable(np.array(raw_4ds).astype(np.float32))
t = Variable(np.array(raw_5th).astype(np.float32).reshape(len(raw_5th),1))

 入力値リストと出力値リストをそれぞれnp.array()に引数として渡して、Numpyのndarray型に変換します。この際にデータ型をnp.float32に変えておきましょう(デフォルトの64bitでは処理できない)。ndarrayに変換したあとは更にChainerで扱うVariableというオブジェクトに変換します。
 これで準備はできたのですが、最後にもう一工夫しておくと便利です。

division_ratio = 0.7 #学習データと検証データの分割比を決める
ind = int(N * division_ratio) #分割境界のインデックス
x_train, x_test = x[:ind], x[ind:]
t_train, t_test = t[:ind], t[ind:]

 読み込んだデータを全て学習データとして使うのではなく、何割かを検証用データとして分ける処理です。上の例では入出力の各リストから70%を学習データとし、残りの30%をテストデータに振り分けています。さて、ここまでくれば、あとは猛烈に学習させるのみです。

トレーニング開始

for epoch in range(3000):
    model.zerograds()
    y = model(x_train)
    loss = F.mean_squared_error(y, t_train)
    loss.backward()
    optimizer.update()
    print("LOSS FUNCTION:",loss.data) #損失関数が減っていく様を表示します

 それでは始めましょう。まずモデルが有する勾配情報を0に初期化した上で、学習データ(x_train)から得られた出力値をyに代入します。出力値yと実測値であるtとの誤差を算出するため、平均二乗誤差(Mean Squared Error: MSE)を返す関数に引数として渡しています。この「誤差を算出するための関数」が損失関数(Loss Function)と呼ばれるもので、作例ではlossという変数名を付けています。
 次に損失関数の勾配を逆伝播を用いて計算します。backward()メソッドを呼ぶことでBackpropagationが開始され、計算結果が各ノードの情報(grad)として記録されます。その後、optimizerのupdate()メソッドが呼ばれると、Adamアルゴリズムに則って誤差が最小になるようにlinksの重みが調整され、gradがdataに反映されて更新処理が完了となります。

 勾配を初期化してからupdateするまでを1 epochと言います。今回はfor文を使って3000 epochの繰り返し学習を記述していますが、損失関数の値を条件指定すればwhile文でも書けそうですね。ちなみに作例は学習データ全体を使ってparameter更新をしているのでバッチ処理ということになります(次回はミニバッチも試してみたいと思います)。

予測性を検証

 人間の代わりに3000回も退屈な数値を勉強させられたANNですが、果たして予測性能はいかほどになっているでしょうか。

from scipy import stats

real, pred = t_test.data, model(x_test).data #実測値と予測値を.dataで入手
for tup in zip(real*120, pred*120): #120倍して秒表示に戻す
    print(tup)

print("real_data",np.mean(real*120),"±",np.std(real*120)/len(t_test)**0.5)
print("pred_data",np.mean(pred*120),"±",np.std(pred*120)/len(x_test)**0.5)

R_val, p_val = stats.pearsonr(real*120, pred*120)
print("PEARSON's R_val:", R_val, "p_val:", p_val)

 学習を終えたmodelに検証用の入力データ(x_test)を渡し、予測値(pred)を出力結果として得ます。実測値と予測値のタプルを表示した上で、Scipyを使って各群の平均値と標準誤差を表示させました。また実測値と予測値の相関分析を行い、ピアソンの相関係数も求めています。
 実行結果は以下のようになりました。

LOSS FUNCTION: 0.0030334105249494314
LOSS FUNCTION: 0.003033286426216364
LOSS FUNCTION: 0.0030331630259752274
LOSS FUNCTION: 0.00303303892724216
(array([ 90.], dtype=float32), array([ 88.0375061], dtype=float32))
(array([ 30.], dtype=float32), array([ 22.96634293], dtype=float32))
(array([ 20.], dtype=float32), array([ 24.61615181], dtype=float32))
(array([ 50.], dtype=float32), array([ 51.244133], dtype=float32))
(array([ 40.], dtype=float32), array([ 47.53696823], dtype=float32))
(array([ 29.], dtype=float32), array([ 29.54800034], dtype=float32))
(array([ 40.], dtype=float32), array([ 47.9940033], dtype=float32))
(array([ 30.], dtype=float32), array([ 31.04858589], dtype=float32))
(array([ 19.], dtype=float32), array([ 24.61615181], dtype=float32))
(array([ 40.], dtype=float32), array([ 48.63508606], dtype=float32))
(array([ 45.], dtype=float32), array([ 46.89588547], dtype=float32))
(array([ 38.], dtype=float32), array([ 34.40187836], dtype=float32))
(array([ 90.], dtype=float32), array([ 79.77006531], dtype=float32))
(array([ 18.], dtype=float32), array([ 9.70597363], dtype=float32))
(array([ 33.], dtype=float32), array([ 34.94280243], dtype=float32))
(array([ 50.], dtype=float32), array([ 48.31797791], dtype=float32))
(array([ 76.], dtype=float32), array([ 58.85263824], dtype=float32))
real_data 43.4118 ± 5.26860451724
prediction 42.89 ± 4.73104620938
PEARSON's R_val: [ 0.95093602] p_val: [  4.82376006e-09]

 相関係数 R=0.95と出ました。なんと驚異的な予測精度でしょう。作戦は大成功と言えそうです。これからは4日間だけまじめに実験をして、5日目はPythonに任せることにしましょう。空いた時間でBiSHの動画でも見ながらビールが飲めたら最高ですね。

終わりに

 最後までお付き合い頂きありがとうございました。今回使用した学習データ・検証データは全て私自身がでっち上げた数値を使っているので、結果を鵜呑みにしないようにお気をつけ下さい。

 また今回はあくまでOne of the etudes for programmingなので、ANNモデルの構想自体が不適切な可能性もあります(例えばday1-day5までのデータは時系列に沿った関連のある数値ですので、LSTMを用いたRecurrent Neural Networkを適用するべきなのかも知れません)。ミニバッチやランダム化など学習方法の工夫も必要でしょうし、Trainerを使った方が汎用性が高まるかも知れません。

 ということで、非理工系出身の素人ではありますが、これからも地道に機械学習の勉強を続けていく予定ですので、暖かく見守って頂ければ嬉しいです。