2017/09/30

3次元一様ランダム回転

3次元の学習データを増やすため、ランダムかつ一様に回転する関数を作りました。手法は色々あるようなのですが、今回は、Shoemake's methodと呼ばれている方法を利用しました。実装例が
http://www.mech.utah.edu/~brannon/public/rotation.pdf
のサブルーチンRNDC (ページA-9)にFortranのコードで示されていましたので、これをPythonで書き直しました。

コードは、ここか下を参照。3次元のベクトルを入力すると、ランダムに回転させたベクトルが返ります。回転行列がrなので、これを返すようにすれば、色々なベクトルを同じように回転できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import numpy as np
import math

def uniform_random_rotate(v):
    x0 = np.random.random()
    y1 = 2*math.pi*np.random.random()
    y2 = 2*math.pi*np.random.random()
    r1 = math.sqrt(1.0-x0)
    r2 = math.sqrt(x0)
    u0 = math.cos(y2)*r2
    u1 = math.sin(y1)*r1
    u2 = math.cos(y1)*r1
    u3 = math.sin(y2)*r2
    coefi = 2.0*u0*u0-1.0
    coefuu = 2.0
    coefe = 2.0*u0
    r = np.zeros(shape=(3, 3))
    r[0, 0] = coefi+coefuu*u1*u1
    r[1, 1] = coefi+coefuu*u2*u2
    r[2, 2] = coefi+coefuu*u3*u3

    r[1, 2] = coefuu*u2*u3-coefe*u1
    r[2, 0] = coefuu*u3*u1-coefe*u2
    r[0, 1] = coefuu*u1*u2-coefe*u3

    r[2, 1] = coefuu*u3*u2+coefe*u1
    r[0, 2] = coefuu*u1*u3+coefe*u2
    r[1, 0] = coefuu*u2*u1+coefe*u3
    return np.dot(v, r)

利用例は以下。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from random_rotate import uniform_random_rotate

uvec = np.array([1,0,0])
data = []
for v in range(10000):
    data.append(uniform_random_rotate(uvec))
data = np.array(data)

fig = plt.figure(figsize=(15, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter3D(data[:,0:1], data[:,1:2], data[:,2:3], s=0.2)
plt.show()

プロット結果が以下。

とりあえず見た感じ、半径1の球の表面に一様分布しているように見えます。一様分布をランダムに生成しているのではなくて、ベクトル(1, 0, 0)を一様ランダムに回転させて得られた図です。

2017/09/23

加速度センサーで行動分類

人の行動のラベル付き加速度センサーのデータを準備したので、行動認識を行ってみます。最初なので、何も工夫せずに、LSTM3層+全結合1層で認識してみます。 具体的には、

model = Sequential()
model.add(LSTM(256, return_sequences=True, input_shape=(timesteps, data_dim)))
model.add(LSTM(256, return_sequences=True))
model.add(LSTM(256))
model.add(Dense(num_classes, activation='softmax'))
という構成です。ユニット数は256個としました。

実装を簡単にするため、学習時も評価時もデータからランダムに切り出してくるようにしました。そのため、同じモデルを学習できていても、精度は変動します。正しく評価するには、評価時は評価用データを全て使う必要があります。

入力データの長さを揃える必要があるため、データの長さの幅は64サンプルで揃えました。サンプリング周波数は32Hzなので、2秒分のデータで行動を推定することになります。データ長が64サンプル未満の場合はデータの先頭に0をパディングします。

切り出しは、クラスが均等に出現するようにしました。また、同一クラス内では、各ファイルが均等に出現するようにしました。データ量に関して、クラス間に偏りがあり、特定のクラスばかり学習されないようにするためです。これは、評価用データの切り出しに関しても同じです。

学習の経過を下図に示します。 バッチサイズ32、1エポックあたりのステップ数は100で、100エポック学習させました。評価のステップ数も100なので、3200個のデータで評価しています。

100エポック目の学習データでの精度は86.7%、評価データでの精度は64.4%となりました。残念ながら、あまりうまく学習できていません。評価データで90%くらいにはなるだろうと考えていましたが、そこまで簡単ではなさそうです。

コードは以下の通りです。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
# This script trains a neural network model using accelerometer data
import keras, os, datetime
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation
from keras.layers import LSTM

TRAIN_EPOCH = 1

def make_model(timesteps, data_dim, num_classes):
    model = Sequential()
    model.add(LSTM(256, return_sequences=True,
                   input_shape=(timesteps, data_dim)))
    model.add(LSTM(256, return_sequences=True))
    model.add(LSTM(256))
    model.add(Dense(num_classes, activation='softmax'))
    print(model.summary())
    return model

def read_accel(filename):
    with open(filename) as f:
        data = []
        for line in f:
            for v in line.rstrip().split(" "):
                if not v.isdigit():
                    print(v, "is not digit")
            data.append([int(v) for v in line.rstrip().split(" ")])
        data = np.array(data)
        return data

# Data generator
def data_gen(filelist, batch_size, window_width):
    # Make a dict for converting a class ID to files
    num_class = 0
    with open(filelist) as f:
        id2f = dict()
        for line in f:
            id, filename = line.split(" ", maxsplit=2)
            id = int(id)
            id2f.setdefault(id, []).append(read_accel(filename.rstrip()))
            num_class = max(num_class, id)
        num_class += 1
    yield num_class # Return the number of classes at first

    # Generate data eternally
    while True:
        data = []
        labels = []
        for i in range(batch_size):
            while True:
                c = np.random.randint(0, num_class) # Select a class
                if c in id2f:
                    break
            # Select a file in the class
            f = np.random.randint(0, len(id2f[c]))
            # Select a window position
            p = np.random.randint(0, max(0, len(id2f[c][f])-window_width)+1)
            # Make a label for the selected class
            labels.append(keras.utils.to_categorical(c, num_classes=num_class))
            # Add selected data
            data.append(id2f[c][f][p:p+window_width,:])
            if(len(data[-1]) < window_width):
                padding = window_width - len(data[-1])
                data[-1] = np.append(np.zeros((padding, data[-1].shape[1]), dtype=np.int),
                                              data[-1], axis=0)
        yield np.array(data), np.array(labels).reshape((batch_size, num_class))

if __name__ == '__main__':
    if not os.path.exists("models"):
        os.mkdir("models")

    batch_size = 32
    window_width = 64
    dg = data_gen("train.list", batch_size, window_width)
    num_classes = next(dg)

    dgv = data_gen("val.list", batch_size, window_width)
    assert num_classes == next(dgv)

    # Make a model
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model = make_model(window_width, 3, num_classes)
    model.compile(loss='categorical_crossentropy',
                  optimizer=opt,
                  metrics=['accuracy'])

    # Train the model
    hist = model.fit_generator(dg, steps_per_epoch=100, 
                               validation_data=dgv, validation_steps=100,
                               epochs=TRAIN_EPOCH).history

    # Write training history into a file
    with open("history.log", mode="a") as f:
        d = datetime.datetime.today()
        f.write("#"+d.strftime("%Y-%m-%d %H:%M:%S")+"\n")
        f.write("#Optimizer :"+str(type(opt))+str(opt.get_config())+"\n")
        f.write("#Data in the last line are calculated by the last model"
                " and not calculated in model.fit()\n")
        f.write("#epoch train-loss train-acc val-loss val-acc\n")
        for v in range(0, TRAIN_EPOCH):
            f.write("{0} {1:10.6f} {2:10.6f} {3:10.6f} {4:10.6f}\n"
                    .format(v, hist["loss"][v], hist["acc"][v],
                            hist["val_loss"][v], hist["val_acc"][v]))
        f.write("\n\n")

生成モデルと識別モデル


機械学習の話です。
分類問題における生成モデルと識別モデルの違いについて、説明のしかたが色々あるようですので、まとめてみます。

まとめ


分類問題の解き方は、2つではなく、3つに分類されます。シンプルな順に、
  • 識別関数
  • 識別モデル
  • 生成モデル
という3つのアプローチが考えられます。

識別関数は、確率は関係なく、入力空間に識別面を直接引く方法です。2クラスなら、いい感じに分割できそうな面を線形・非線形を問わず引き、面で区切られたの一方の領域にあるサンプルをクラスA、もう一方をクラスBとする、という方法です。

識別モデルは、事後確率\(P(C_i|x)\)を直接モデル化する方法です。入力サンプル\(x\)に対応する事後確率をモデルを使って計算してから、最も事後確率が高いクラスを識別したクラスとします(\(C_i\)は\(i\)番目のクラスです)。

生成モデルは、事前確率\(P(C_i)\)と確率密度関数\(P(x|C_i)\)をモデル化し、それをベイズの定理を使って事後確率\(P(C_i|x)\)を計算し、最も事後確率が高いクラスを識別したクラスとする方法です。確率密度関数をつかって、クラス\(C_i\)に属するサンプルをいくらでも生成できることから、生成モデルと呼ばれます。

これでは簡略化しすぎていて分からないという場合は、全文をどうぞ。

2017/09/21

日常の活動を記録した加速度データセット

手首に付けた加速度センサーで日常の活動を記録したデータセットが公開されています。

https://archive.ics.uci.edu/ml/datasets/Dataset+for+ADL+Recognition+with+Wrist-worn+Accelerometer

大量のデータというわけではありませんが、合計839ファイルあり、ファイル単位でラベルが付いています。 クラス数は14です。具体的には、

  • 歯磨き
  • 階段を上る
  • 髪を梳く
  • 階段を下りる
  • 飲み物を飲む
  • 肉を食べる
  • スープを飲む
  • ベッドから起きる
  • ベッドで横になる
  • 水を注ぐ
  • 椅子に座る
  • 椅子から立ち上がる
  • 電話する
  • 歩く
です。バランス良く収録されているわけではなく、クラスによってファイル数に偏りがあります。

このデータとKerasで色々と試してみたいので、学習用データと検証用データに分けてみました。 きちんと色々と調べるには交差検証で行ったほうが良いのですが、少し試してみたいだけですので、4/5を学習用データに、1/5を検証用データとします。

分割には次のスクリプトを用いました。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#!/bin/bash
DATADIR="ADL_Dataset/HMP_Dataset"

# Collect all files
find $DATADIR -iname "Accelerometer*.txt" | grep -v MODEL | sed 's#^.*/c/#/#' > all.list

# Make a category list
sed 's/.*Accelerometer-[0-9]\{4\}-[0-9]\{2\}-[0-9]\{2\}-[0-9]\{2\}-[0-9]\{2\}-[0-9]\{2\}-//' all.list | \
  sed 's/-.*//' | sort | uniq | awk 'BEGIN{c=0}{print c" "$1; c+=1}' > id2cat.list

# Make a category-ID mapper
declare -A CAT2ID
exec < id2cat.list
while read id cat
do
    CAT2ID[$cat]=$id
done
echo ${CAT2ID[getup_bed]}

# Add category ID into all.list and randomize it
{
while read line
do
    category=${line%-*}
    category=${category##*-}
    echo ${CAT2ID[$category]}" "$line 
done < all.list
} | sort -R > all.random.list

# Make training set and validation set
N=$(wc -l all.list | awk '{print $1}')
TRAIN=$((N*4/5))
VAL=$((N-TRAIN))
head -n $TRAIN all.random.list > train.list
tail -n $VAL all.random.list > val.list

実際に実行してみると、少ないクラスのファイルが検証用リスト(val.list)に含まれませんでしたので、何度が実行して、train.list、val.listともに全てのクラスが含まれるようにしました。

生成したファイルリストの中身については train.listval.list を参照ください。

train.listとval.listの各クラスのファイル数は、

Class ID # of train
files
# of val
files
0 10 2
1 74 28
2 27 4
3 35 7
4 83 17
5 2 3
6 1 2
7 80 21
8 22 6
9 83 17
10 82 18
11 85 17
12 11 2
13 76 24

となりました。クラス5と6はファイル数が学習データのほうが少ないので、学習が難しそうです。

2017/09/17

Pythonのyield構文

Pythonのyield構文について、どういう動作をしているのか簡単に調べてみました。Python 3.5で調べています。

forループで使用


3回yieldで返す関数genを作成し、forループで値を取得してみました。
1
2
3
4
5
6
7
8
9
def gen():
    yield 1
    yield 3
    yield 2
    
g = gen()
for v in g:
    print(v)
print("End")
実行すると、
1
3
2
End
が得られました。yieldで返した値が順に表示されます。

次に、無限に出力できるようにしてみました。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def gen():
    x = 0
    while True:
        x += 1
        yield x

g = gen()
for v in g:
    print(v)
print("End")
実行すると、1, 2, 3,...と順に1ずつ値を増やしながら停止するまで無限に整数を出力し続けます。

nextを使用


forループを使わず、1回ずつ呼び出してみます。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
def gen():
    yield 1
    yield 3
    yield 2
    
g = gen()
try:
    c = 1
    print(next(g))
    c = 2
    print(next(g))
    c = 3
    print(next(g))
    c = 4
    print(next(g)) # raise a StopIteration exception
    c = 5
    print(next(g))
    c = 6
except StopIteration:
    print("End")
print("c =", c)
実行結果は、
1
3
2
End
c = 4
となります。yieldが返すことが出来るのは3回だけなので、4回目のnextで例外が発生します。

次に、グローバル変数aの値を変更・参照しながら、gen()を呼び出してみます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
a = 0
def gen():
    c = 1
    yield c + a
    c += 2
    yield c + a
    c -= 1
    yield c + a
    
g = gen()
try:
    c = 1
    print(next(g))
    a = 10; c = 2
    print(next(g))
    a = 100; c = 3
    print(next(g))
    c = 4
    print(next(g)) # raise a StopIteration exception
    c = 5
    print(next(g))
    c = 6
except StopIteration:
    print("End")
print("c =", c)
実行結果は、
1
13
102
End
c = 4
となります。nextを呼び出す前にaを変更していますが、反映されていますね。cは固定値ではないですが、正しく前の値を覚えてくれています。

動作の仕組みは?


内部の仕組みを調べたわけではありませんので、以下は、事実と異なる可能性が高いのでご注意を。

上のコードのように、nextを呼び出すとgenが実行され、yieldに到達すると、yieldで指定されている値を返す直前に、yieldの次に実行すべき行番号や関数内でのみ使われている変数が保存されているのでしょう。次にnext経由で当該関数が呼び出されると、保存していた行番号や変数の値を復元し、続きを実行しているのだと思われます。

2017/09/15

ダイソーの乾電池

ずいぶん昔に買った、ダイソーの乾電池を使おうと思ってみてみたら
色々残念なことになっていました。

型番は
スーパーアルカリ乾電池
単4型AAA LR03
KW-13-80
T097
アルカリ電池No.14
で、使用推奨期限が2018年8月。

6本入りのうち、4本が元のパッケージに入った状態で新品として残っていました。
しかし、保管状況が悪かったのか、1本は液漏れ、2本は放電してしまっていて
電圧が低すぎて使えず、なぜか1本だけ正常でした。

品質のばらつきが大きすぎですね。

安い電池に文句を言っても仕方がありませんが、
ストック用に買うのはいまいちなようです。

ニューラルネットワークで正弦波にフィット(その2)

前回の正弦波フィットでは、フィット具合がいまいちでした。そこで、どういうパラメータにすると、どの程度まで正弦波にフィットできるのかを調べてみました。

前回と同様、\(y=\sin(2x)+n\) (\(n\)がノイズ)からサンプルしたデータにフィッティングをして、\(y=\sin(2x)\)を学習できるか試してみました。\(x\)の値の生成範囲は[-10, 10)です。モデルのパラメータは下表の値の組み合わせで、合計72パターン試しました。

Parameter nameValues
Activation typetanh, relu
# of units for each layer1, 2, 4, 8, 16, 32
# of hidden layers1, 2, 3, 4, 5, 6

活性化関数が tanh の場合


まずは、Activationをtanhにした場合です。u1は各レイヤーのユニット数が1で、f2は元の関数が\(\sin(2x)\)であることを表します。図中の各系列名のlayerの数値は隠れ層-1を表します。背景の薄い青色の正弦波が学習データの元のになった正弦波です。これに近いほど、うまく学習できていることになります。
ユニット数が1なので、layerが1だと、 \[y=w_2 \tanh(w_1 x + b_1) + b_2\] になるため、単に\(\tanh\)の平行移動と拡大縮小ができるだけです。

ユニット数が2つになると、次の図のように、原点の前後の正弦波の山に近づくようにフィットできるようになります。

隠れ層が1(図中ではlayer=0)の場合、\(\tanh\)を2つ、拡大縮小移動して重ね合わせて作ることができます。 隠れ層が2なら、\(\tanh\)を4つ重ね合わせていることになります。隠れ層が1増えるたびに使える\(\tanh\)は2倍になるので、フィットしやすくなります。しかし、この数では原点の前後の正弦波の山にフィットできるのみでした。

ユニット数が4つになると、原点から2つ目までの山と谷にフィットできそうになります(下図)。

ユニット数が8つになると、隠れ層が6層なら、原点の前後3個の山と谷にフィットでき、4つめはもう少しでフィットできそうです(下図)。

ユニット数が16個になると、隠れ層が6層なら、原点の前後5個の山と谷にフィットできます。

ユニット数が32個になると、隠れ層が4層以上なら原点の前後7個の山と谷にフィットできます。これは、学習データの存在範囲の正弦波にフィットできていることを表します。

活性化関数が relu の場合

次に、reluの場合です。ユニット数が1の場合は、\(\tanh\)の場合と同じで、ReLUそのものの形を拡大縮小と平行移動させたグラフになります。

ユニット数が2に増えると、隠れ層の数によっては折れ曲がる箇所が増えます。

ユニット数が4に増えると、山が大きめになります。

ユニット数が8に増えると、隠れ層の数によっては、原点の前後の1つ目の山または谷に近い形になります。

ユニット数が16に増えると、隠れ層の数が6の場合、学習データの\(x\)の範囲である[-10, 10)の範囲では、不完全ながらも正弦波に追従しています。

ユニット数が32の場合、隠れ層の数が5または6なら、正弦波をかなり近似できているように見えます。

ユニット数と隠れ層の数がそこそこ大きくなれば、活性化関数が\(\tanh\)でもReLUでも今回の条件の正弦波からサンプリングしたデータから元の正弦波に近いモデルを学習できることが分かりました。

ソースコード


学習時のソースコードは以下の通りです。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# This script trains a neural network model using MNIST.
import keras, os, math
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation

TRAIN_EPOCH = 100

# Model for fitting
#  nunit : # of units for each dense layer
#  nlayer: # of hidden layers - 1
#  act   : activation type
def approx_model(nunit, nlayer, act):
    model = Sequential()
    model.add(Dense(nunit, input_shape=(1,)))
    model.add(Activation(act))
    for v in range(nlayer):
        model.add(Dense(nunit))
        model.add(Activation(act))
    model.add(Dense(1))
    print(model.summary())
    return model

# Make samples on sine curve with noise
def sine_dist(freq):
    xd = np.random.rand(10000,)*20-10
    yd = [math.sin(x*freq)+np.random.normal(0, 1) for x in xd]
    return xd, yd

if __name__ == '__main__':
    if not os.path.exists("models"):
        os.mkdir("models")

    for act in ['tanh', 'relu']:
        for nunit in [1, 2, 4, 8, 16, 32]:
            for nlayer in [0, 1, 2, 3, 4, 5]:
                # Make a model
                model = approx_model(nunit, nlayer, act)
                opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                              epsilon=1e-08, decay=1e-4)
                model.compile(optimizer=opt, loss='mean_squared_error')
                model.save_weights("initial-weights.hdf5")

                for freq in [2]:
                    # Make data for training and validation
                    x_train, y_train = sine_dist(freq)
                    x_val, y_val = sine_dist(freq)
                    
                    # Train and save a model
                    model.load_weights("initial-weights.hdf5")
                    model.fit(x_train, y_train, validation_data=(x_val, y_val),
                              epochs=TRAIN_EPOCH, batch_size=32)
                    model.save("models/sine-{0}-u{1}-l{2}-f{3}.hdf5".format(act, nunit, nlayer, freq))

また、Jupyter notebookでグラフを表示させるコードは次の通りです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
%matplotlib inline
import keras, math
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(-20, 20, 0.001)
for act in ['tanh', 'relu']:
    for u in [1, 2, 4, 8, 16, 32]:
        plt.figure(figsize=(15, 10))
        plt.title('{0}-u{1}-f2'.format(act, u))
        plt.grid()
        yc = np.array([math.sin(xi*2) for xi in x])
        plt.plot(x, yc, color='lightblue')
        for l in range(6):
            model = keras.models.load_model(
                'models/sine-{0}-u{1}-l{2}-f2.hdf5'.format(act, u, l))
            y = model.predict(x)
            plt.plot(x, y, label='layer={0}'.format(l))
        plt.legend(loc='upper left', frameon=False)
        plt.show()

2017/09/13

ニューラルネットワークで正弦波にフィット

ニューラルネットワークを使って、正弦波のデータ\(y=\sin(fx)\)にノイズをのせた、\(y=\sin(fx)+n\) (\(n\)がノイズ)からサンプルしたデータにフィッティングをして、\(y=\sin(fx)\)を学習できるか試してみました。\(f\)は周波数の大小を表します。

モデルは、直線フィッティングのときと同じです。

学習データを[-10,10)の範囲で一様乱数で10000サンプルとって、100エポック学習しました。 学習できたモデルの出力は\(f\)の値ごとに下図のようになりました。横軸がx、縦軸がyです。 薄い水色の線は、ノイズを加える前の学習データの生成に使った正弦波です。紺色の線が、学習したニューラルネットワークが生成した曲線です。

\(f=0.5\)

\(f=1\)

\(f=2\)

\(f=4\)

\(f=8\)

\(f=16\)

\(f=0.5\)や\(f=1\)はともかく、それ以降はフィットできていません。特に\(f=4\)以降は完全に崩れています。\(f=16\)では、そもそも周波数が高すぎて、図では元の正弦波がつぶれてしまっています。きちんとフィットさせるためにはニューラルネットワークをもう少し工夫する必要があることが分かりました。

\(f=4\)のときの学習データに使ったデータの生成と同じ条件で、ただし乱数違いで生成したデータをプロットしたものが下図です。横軸の範囲が異なることに注意してください。

\(f=4\) sampled with noise
まだ正弦波を見て取れるのですが、今回使ったニューラルネットワークでは難しいようです。

同様に、\(f=16\)のときの学習データに使ったデータの生成と同じ条件で、ただし乱数違いで生成したデータをプロットしたものが下図です。一応、正弦波の名残は残っていますが、元が正弦波とは言いがたいです。

\(f=16\) sampled with noise

今回の実験のモデルの学習時に使ったコードは以下です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# This script trains a neural network model using MNIST.
import keras, os, math
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation

TRAIN_EPOCH = 100

def approx_model():
    model = Sequential()
    model.add(Dense(32, input_shape=(1,)))
    model.add(Activation('tanh'))
    model.add(Dense(32))
    model.add(Activation('tanh'))
    model.add(Dense(1))
    print(model.summary())
    return model

def sine_dist(freq):
    xd = np.random.rand(10000,)*20-10
    yd = [math.sin(x*freq)+np.random.normal(0, 1) for x in xd]
    return xd, yd

if __name__ == '__main__':
    if not os.path.exists("models"):
        os.mkdir("models")

    # Make a model
    model = approx_model()
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model.compile(optimizer=opt, loss='mean_squared_error')
    model.save_weights("initial-weights.hdf5")

    for freq in [0.5, 1, 2, 4, 8, 16]:
        # Make data for training and validation
        x_train, y_train = sine_dist(freq)
        x_val, y_val = sine_dist(freq)
        
        # Train and save a model
        model.load_weights("initial-weights.hdf5")
        model.fit(x_train, y_train, validation_data=(x_val, y_val),
                  epochs=TRAIN_EPOCH, batch_size=32)
        model.save("models/sine-{0}.hdf5".format(freq))

グラフの描画に使ったコードは以下です(\(f=1\)の場合)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
%matplotlib inline
import keras, math
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(-30, 30, 0.001)
model = keras.models.load_model('models/sine-1.hdf5')
y = model.predict(x)
yc = np.array([math.sin(xi*1) for xi in x])
plt.grid()
plt.plot(x, yc, color='lightblue')
plt.plot(x, y, color='midnightblue')
plt.show()

2017/09/11

ニューラルネットワークで直線にフィット

ニューラルネットワークを使って、直線のデータ\(y=x\)にノイズをのせた、\(y=x+n\) (\(n\)がノイズ)からサンプルしたデータにフィッティングをして、\(y=x\)を学習できるか試してみました。

Kerasを使って学習します。モデルは、

model = Sequential()
model.add(Dense(32, input_shape=(1,)))
model.add(Activation('tanh'))
model.add(Dense(32))
model.add(Activation('tanh'))
model.add(Dense(1))
としました。ただの直線フィッティングにしてはゴージャスな構成です。

学習データを[-10,10)の範囲で一様乱数で10000サンプルとって、100エポック学習したところ、学習できたモデルの出力は下図のようになりました。横軸が\(x\)、縦軸が\(y\)です。

学習データのある\(x\)が-10から10の範囲では目で見る限りにおいてノイズに惑わされず直線になっていますが、その範囲外は明らかに\(y=x\)ではありません。

学習用のコードは以下です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# This script trains a neural network model by using y=x with noise.
import keras, os
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation

TRAIN_EPOCH = 100

def approx_model():
    model = Sequential()
    model.add(Dense(32, input_shape=(1,)))
    model.add(Activation('tanh'))
    model.add(Dense(32))
    model.add(Activation('tanh'))
    model.add(Dense(1))
    print(model.summary())
    return model

if __name__ == '__main__':
    if not os.path.exists("models"):
        os.mkdir("models")

    # Make data for training and validation
    x_train = (np.random.rand(10000,)*20-10).reshape(10000, 1)
    y_train = np.array([x+np.random.normal(0, 1) for x in x_train]).reshape(10000, 1)
    x_val = (np.random.rand(10000,)*20-10).reshape(10000, 1)
    y_val = np.array([x+np.random.normal(0, 1) for x in x_val]).reshape(10000, 1)

    # Make a model
    model = approx_model()

    # Train the model
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model.compile(optimizer=opt, loss='mean_squared_error', metrics=['mse'])
    hist = model.fit(x_train, y_train, validation_data=(x_val, y_val),
                     epochs=TRAIN_EPOCH, batch_size=32).history

    # Evaluate the last model by train and validation data
    score_train = model.evaluate(x_train, y_train, batch_size=32)
    score_val = model.evaluate(x_val, y_val, batch_size=32)
    model.save("models/linear.hdf5")

図の表示用のコードは、Jupyter notebook用のコードで、以下の通りです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
%matplotlib inline
import keras
import matplotlib.pyplot as plt
import numpy as np
x = np.array([v/10 for v in range(-300,300)])
model = keras.models.load_model('models/linear.hdf5')
y = model.predict(x)
plt.grid()
plt.scatter(x, y, s=1)
plt.show()