FaxOCR手書き数字データの認識 その2

%matplotlib inline
import pylab as plt
import pandas as pd
import numpy as np

概要

前回、FaxOCRという手書き数字のデータの認識をやった。 認識自体はぼちぼちできたが、MNISTデータで学習させたCNNで認識を行うといまいちだったのが気になった。 バグやミスの可能性を潰してもう一度やってみたけど、同様にうまくいかなかった。 データを見ていると文字のサイズが異なっていることに気づいた。 サイズを統一してやってみると、MNISTデータで学習させたCNNである程度正しく予測することができた。

はじめに

前回、FaxOCRという手書き数字のデータの認識をやった。 学習データを回転させるなどしてデータを増やして、CNNを使って学習させると、96%ぐらいの精度で予測することができた。 一方で、MNISTのデータを使って学習させたCNNでは70%弱でしか当てることができなかった。 自分のプログラムや計算にミスがある可能性も考えながら、色々やる。 今回はMNISTデータとの違いを見ることが目的なので、FaxOCRのデータは全て前処理されていない元画像データ(numbers-sample, mustread)を使った。 FaxOCRについてはこちら

ミスの可能性をなくす

全部画像データに変換する

前回はMNISTのデータをバイナリで読み込んでCNNやt-SNEに突っ込んでいた。やらかしているならここだなと思って、とりあえず全部pngにした。

import pylab as plt
from sklearn.datasets import fetch_mldata
import numpy as np


def save(image, name):
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    imgplot = ax.imshow(image, cmap=plt.cm.Greys)
    imgplot.set_interpolation('nearest')
    ax.xaxis.set_ticks_position('top')
    ax.yaxis.set_ticks_position('left')
    plt.imsave(name)
    # plt.savefig(name) # savefigだとグラフの軸も描画される

mnist = fetch_mldata('MNIST original', data_home=".")
y = mnist.target
X = - mnist.data.reshape(len(y), 28, 28) + 255

counter = np.zeros(10)
from itertools import izip
for image, label in izip(X, y):
    label = int(label)
    plt.imsave("%d_%d.png" % (label, counter[label]), image, cmap=plt.cm.gray)
    counter[label] += 1

プレビューが死ぬほど重いけど、特に問題なさそう。

FaxOCR

f:id:ksknw:20160430171915p:plain

MNIST

f:id:ksknw:20160430172107p:plain

FaxOCRのデータは以下のようにImageMagicを使って、28x28に変えた。

for i in *
do              
convert -resize 28x28! $i ../28/$i                                       
done

CNNのコードを書き直す

ネットにあったコードを行き当りばったりな感じで編集してコードを書いていた。 コメントアウトで条件変えたり、ミスしてそうなところがあったので、そこそこちゃんと書き直す。(書きなおしたあと色々あってまた行き当りばったり的なコードになっているけど気にしない)

learn.py

# -*- coding: utf-8 -*-
import numpy as np
import glob
import cv2 as cv
from itertools import izip
import random

from cnn import cnn


def read_imgs(dirname, labelpos=1):
    imgs = []
    labels = []
    for img_file in glob.glob(dirname + "/*.png"):
        imgs.append((255 - cv.imread(img_file, flags=0)) / 255.0)
        labels.append(int(img_file[len(dirname) + labelpos]))
    return np.array(imgs), np.array(labels)


def learn(train="faxocr", imsize="28"):
    X_test, y_test = read_imgs("./data/faxocr/test/%s" % imsize)

    # if train == "mnist":
    #     assert(int(imsize) == 28)
    X_train, y_train = read_imgs("./data/%s/train/%s" % (train, imsize))

    size = tuple(np.array([X_train[0].shape[1], X_train[0].shape[0]]))

    if train == "faxocr":
        new_imgs = []
        new_labels = []
        for img, label in izip(X_train, y_train):
            for i in range(20):
                rad = (random.random() - 0.5) * 0.5
                pos1 = (random.random() - 0.5) * 5
                pos2 = (random.random() - 0.5) * 5
                mat = np.float32([[np.cos(rad), -1 * np.sin(rad), pos1],
                                  [np.sin(rad), np.cos(rad), pos2]])
                dst = cv.warpAffine(img, mat, size, flags=cv.INTER_LINEAR)
                new_imgs.append(dst)
                new_labels.append(label)
        X_train = np.r_[X_train, new_imgs]
        y_train = np.r_[y_train, new_labels]
    cnn(X_train, y_train,
        X_test,  y_test,
        #        "./results/" + train + "_%s_" % imsize, size=imsize)
        "./results/" + train + "_%s_" % imsize, size=28)

if __name__ == '__main__':
    learn(train="mnist", imsize="trim_28")

cnn.py

# coding: utf-8
import numpy as np
import chainer
from chainer import cuda
import chainer.functions as F
from chainer import optimizers
import time


def cnn(train_data, train_label,
        test_data,  test_label,
        resultname_header,
        n_epoch=50, batchsize=100,
        size=28):

    cuda.check_cuda_available()
    xp = cuda.cupy

    N = train_label.size
    N_test = test_label.size

    train_data = train_data.reshape(len(train_label), -1)
    train_data = train_data.astype(xp.float32)
    train_label = train_label.astype(xp.int32)
    test_data = test_data.reshape(len(test_label), -1)
    test_data = test_data.astype(xp.float32)
    test_label = test_label.astype(xp.int32)

    train_data = train_data.reshape((len(train_data), 1, size, size))
    test_data = test_data.reshape((len(test_data), 1, size, size))

    print test_data.shape
    print train_data.shape

    print test_data.mean()
    print train_data.mean()

    if size == 28:
        model = chainer.FunctionSet(conv1=F.Convolution2D(1, 20, 5),
                                    conv2=F.Convolution2D(20, 50, 5),
                                    l1=F.Linear(800, 500),
                                    l2=F.Linear(500, 10))

    else:
        model = chainer.FunctionSet(conv1=F.Convolution2D(1, 20, 3),
                                    conv2=F.Convolution2D(20, 50, 3),
                                    l1=F.Linear(6050, 800),
                                    l2=F.Linear(800, 10))

    cuda.get_device(0).use()
    model.to_gpu()

    def forward(x_data, y_data, train=True):
        x, t = chainer.Variable(x_data), chainer.Variable(y_data)
        h = F.max_pooling_2d(F.relu(model.conv1(x)), 2)
        h = F.max_pooling_2d(F.relu(model.conv2(h)), 2)
        h = F.dropout(F.relu(model.l1(h)), train=train)
        y = model.l2(h)
        if train:
            return F.softmax_cross_entropy(y, t)
        else:
            return F.accuracy(y, t)

    optimizer = optimizers.Adam()
    # optimizer = optimizers.RMSprop()
    optimizer.setup(model)

    fp1 = open(resultname_header + "accuracy_row.txt", "w")
    fp2 = open(resultname_header + "loss_row.txt", "w")

    fp1.write("epoch\ttest_accuracy\n")
    fp2.write("epoch\ttrain_loss\n")

    
    start_time = time.clock()
    for epoch in range(1, n_epoch + 1):
        print "epoch: %d" % epoch

        perm = np.random.permutation(N)
        sum_loss = 0
        for i in range(0, N, batchsize):
            x_batch = xp.asarray(train_data[perm[i:i + batchsize]])
            y_batch = xp.asarray(train_label[perm[i:i + batchsize]])

            optimizer.zero_grads()
            loss = forward(x_batch, y_batch)
            loss.backward()
            optimizer.update()
            sum_loss += float(loss.data) * len(y_batch)

        print "train mean loss: %f" % (sum_loss / N)
        fp2.write("%d\t%f\n" % (epoch, sum_loss / N))
        fp2.flush()

        sum_accuracy = 0
        for i in range(0, N_test, batchsize):
            x_batch = xp.asarray(test_data[i:i + batchsize])
            y_batch = xp.asarray(test_label[i:i + batchsize])

            acc = forward(x_batch, y_batch, train=False)
            sum_accuracy += float(acc.data) * len(y_batch)

        print "test accuracy: %f" % (sum_accuracy / N_test)
        fp1.write("%d\t%f\n" % (epoch, sum_accuracy / N_test))
        fp1.flush()

    end_time = time.clock()
    print end_time - start_time

    fp1.close()
    fp2.close()

    import cPickle
    model.to_cpu()
    cPickle.dump(model, open(resultname_header + "model_cnn_row.pkl", "wb"), -1)

学習結果

ミスしそうなところはだいたい直したので、もう一度CNNを使って学習させてみた。以下は結果。

FaxOCR -> FaxOCR (28x28)

accuracy = pd.read_csv("./results/faxocr_28_accuracy_row.txt", sep="\t")
loss = pd.read_csv("./results/faxocr_28_loss_row.txt", sep="\t")

fig = plt.figure(figsize=[10,10])
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot()
plt.title("FaxOCR->FaxOCR")
plt.show()

f:id:ksknw:20160430172219p:plain

これは前回と同様な感じ。だいたいOK。

MNIST -> FaxOCR (28x28)

accuracy = pd.read_csv("./results/mnist_28_accuracy_row.txt", sep="\t")
loss = pd.read_csv("./results/mnist_28_loss_row.txt", sep="\t")

fig = plt.figure(figsize=[10,10])
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot()
plt.title("MNIST->FaxOCR")
plt.show()

f:id:ksknw:20160430172235p:plain

残念ながらこれも前回と同じ。どうもミスとかではなく、普通にダメそう。 loss自体は下がっているので、MNISTとFaxOCRのデータが何かしら違うことが原因っぽい。

データをみる

個々のデータを目で見ていても、特に不自然なところはないように感じたので、色々絵を描いてみて考えることにした。

from learn import read_imgs

mnist_data,  mnist_label  = read_imgs("./data/mnist/train/28")
faxocr_data, faxocr_label = read_imgs("./data/faxocr/train/28")

print mnist_data.mean()
print faxocr_data.mean()
0.131017957897
0.079225616316

画素の平均値が違っているのが少し気になる。

t-SNE

まずは僕の大好きなt-SNEで絵を描く。 前回はMNISTを1000個とFaxOCRのテストデータを使って可視化したけど、今回はMNISTデータ7000個とFaxOCRの学習データ6709個を使った。

num_mnist = 7000

import random
indecies = random.sample(range(len(mnist_data)), num_mnist)

data = np.r_[mnist_data[indecies].reshape(num_mnist, -1), faxocr_data.reshape(len(faxocr_data),-1)]

from sklearn.manifold import TSNE
model = TSNE(n_components=2)
tsned = model.fit_transform(data)
label = np.r_[["b" for i in range(num_mnist)], ["r" for i in range(len(faxocr_data))]]
plt.figure(figsize=(30,30))
plt.scatter(tsned[:,0], tsned[:,1], c=label, linewidths=0)
plt.show()

f:id:ksknw:20160430172255p:plain

青がMNIST、赤がFaxOCR。分離してるなぁーって感じの図。わかりにくいので、数字ごとに図を描いてみる。

fig = plt.figure(figsize=(30,40))
for num in range(10):
    fig.add_subplot(4,3,num+1)
    label = np.r_[["b" if i==num else "w" for i in mnist_label[indecies]], ["r" if i==num else "w" for i in faxocr_label]]
    plt.scatter(tsned[:,0], tsned[:,1], c=label, linewidths=0, alpha=0.6, marker=".")
    plt.title(str(num))
plt.show()

f:id:ksknw:20160430172315p:plain

まあだめでしょうねって感じの図になった。

ちょっと気になったので、FaxOCRデータとMNISTデータそれぞれでt-SNEして図を描いてみる。

model_faxocr = TSNE(n_components=2)
tsned_faxocr = model_faxocr.fit_transform(faxocr_data.reshape(len(faxocr_data),-1))
model_mnist = TSNE(n_components=2)
tsned_mnist = model_faxocr.fit_transform(mnist_data[indecies].reshape(num_mnist,-1))

plt.figure(figsize=(20,10))
plt.subplot(121)
plt.title("FaxOCR")
plt.scatter(tsned_faxocr[:,0], tsned_faxocr[:,1], c=faxocr_label, linewidths=0, marker=".")
plt.subplot(122)
plt.title("MNIST")
plt.scatter(tsned_mnist[:,0], tsned_mnist[:,1], c=mnist_label[indecies], linewidths=0, marker=".")
plt.show()

f:id:ksknw:20160430172358p:plain

なんだこれは。MNISTの方はすごいきれいに分かれているのに。

ちなみに回転させたデータを入れたFaxOCRのデータを可視化すると以下。

new_imgs = []
new_labels = []
size = tuple(np.array([faxocr_data[0].shape[1], faxocr_data[0].shape[0]]))

from itertools import izip
import cv2 as cv
for img, label in izip(faxocr_data, faxocr_label):
    for i in range(20):
        rad = (random.random() - 0.5) * 0.5
        pos1 = (random.random() - 0.5) * 5
        pos2 = (random.random() - 0.5) * 5
        mat = np.float32([[np.cos(rad), -1 * np.sin(rad), pos1],
                          [np.sin(rad), np.cos(rad), pos2]])
        dst = cv.warpAffine(img, mat, size, flags=cv.INTER_LINEAR)
        new_imgs.append(dst)
        new_labels.append(label)
many_data = np.r_[faxocr_data, new_imgs]
many_label = np.r_[faxocr_label, new_labels]

fax_indecies = random.sample(range(len(many_data)), num_mnist) 

model_many = TSNE(n_components=2)
tsned_many = model_many.fit_transform(many_data.reshape(len(many_data),-1)[fax_indecies])

plt.figure(figsize=(10,10))
plt.scatter(tsned_many[:,0], tsned_many[:,1], c=many_label[fax_indecies], linewidths=0, marker=".")
plt.show()

f:id:ksknw:20160430172413p:plain

このデータ分離できるというのはCNNがすごいのかt-SNEがいまいちなのかなんなんだ。 というかMNISTのデータはなんであんなに綺麗に描けるんだ。

どの点がどの画像なのかを見る。

数字ごとに分けて書いた図を見ているとどうもMNISTとFaxOCRで被っている点もある。その点がどの点なのかを見ることで、何かわかるんじゃないかと思って、以下のように可視化した。

MNISTとFaxOCRの点が重なっているところで緑とかになっているのは、直すのがめんどうなだけなので気にしないでほしい。(これ系の図のもっと楽な書き方を知っている人がいたら教えて欲しいです…)

img_size = 28 * 100
label = np.r_[mnist_label[indecies].reshape(num_mnist, -1), faxocr_label.reshape(len(faxocr_data),-1)]
positions = (tsned - tsned.min()) *img_size/(tsned.max() - tsned.min())
plt.figure(figsize=(50*2,50*5))
for num in range(10):
    plt.subplot(5,2,num+1)
    img = np.ones((img_size, img_size, 3))
    for i, pos in enumerate(positions):
        if label[i] != num:
            continue
        temp = data[i].reshape(28,28)
        if i < num_mnist:
            temp = np.c_[ np.zeros([784]), data[i], data[i]]
        else:
            temp = np.c_[data[i], data[i],  np.zeros([784])]
        temp = temp.reshape(28,28,3)
      
        
        if pos[0]-14<0 or pos[0]+14>img_size or pos[1]-14<0 or pos[1]+14>img_size:
            continue
        img[pos[0]-14:pos[0]+14, pos[1]-14:pos[1]+14, :] -= temp
        
    plt.imshow(img)
    plt.title(num)
    #plt.savefig("./results/tsne%d.png"%num)
    #plt.savefig("./results/tsne%d.eps"%num)
plt.show()

f:id:ksknw:20160430173224j:plain

図をぼんやり眺めていると、「これ字のサイズが違うだけじゃね」って思い始めた。 (図が縮小されてわからないと思うので、こちらに元サイズの画像をおいた。ちなみに5492x13993ある。)

画像中の文字の大きさを統一する

FaxOCRのデータは元データをそのまま入力しているので、大きさが統一されていない。MNISTのデータも色々な大きさの数字が混ざっているんだろうと思い込んでいたんだけど、どうもそうでもないみたい。 ちゃんと公式サイトみると、

"The digits have been size-normalized and centered in a fixed-size image."

って書いてあった。

というわけで同様の処理をFaxOCRのデータにも行う。 トリミングして数字を画像の中心にもってきてってやるの、ちゃんとプログラム書くとそこそこめんどうだなぁと思っていたけど、 ImageMagickで探してみたら意外とあったので、以下のようにコマンドを叩いてぱぱっとやる。 こちらこちらを参考にした。

for i in *
do              
convert -fuzz %60 -trim $i ../trim/$i
done
cd ../trim
for i in *
do              
convert $i -background white -gravity center -thumbnail 28x28 -extent 28x28 ../trim_28/$i
done

余白の設定がめんどうだったので、FaxOCRだけでなくMNISTのデータにも適応した。 できた画像は以下のような感じ f:id:ksknw:20160430174546p:plain

再びt-SNE

サイズを調整した画像を再びt-SNEに突っ込んで可視化する。

from learn import read_imgs

mnist_data,  mnist_label  = read_imgs("./data/mnist/train/trim_28")
faxocr_data, faxocr_label = read_imgs("./data/faxocr/train/trim_28")

print mnist_data.mean()
print faxocr_data.mean()
num_mnist = 7000

import random
indecies = random.sample(range(len(mnist_data)), num_mnist)

data = np.r_[mnist_data[indecies].reshape(num_mnist, -1), faxocr_data.reshape(len(faxocr_data),-1)]
0.298007055179
0.206816194953
from sklearn.manifold import TSNE
model = TSNE(n_components=2)
tsned = model.fit_transform(data)
label = np.r_[["b" for i in range(num_mnist)], ["r" for i in faxocr_data]]
plt.figure(figsize=(30,30))
plt.scatter(tsned[:,0], tsned[:,1], c=label, linewidths=0)
plt.show()

f:id:ksknw:20160430174622p:plain

img_size = 28 * 100
label = np.r_[mnist_label[indecies].reshape(num_mnist, -1), faxocr_label.reshape(len(faxocr_data),-1)]
positions = (tsned - tsned.min()) *img_size/(tsned.max() - tsned.min())
plt.figure(figsize=(50*2,50*5))
for num in range(10):
    plt.subplot(5,2,num+1)
    img = np.ones((img_size, img_size, 3))
    for i, pos in enumerate(positions):
        if label[i] != num:
            continue
        temp = data[i].reshape(28,28)
        if i < num_mnist:
            temp = np.c_[ np.zeros([784]), data[i], data[i]]
        else:
            temp = np.c_[data[i], data[i],  np.zeros([784])]
        temp = temp.reshape(28,28,3)
      
        
        if pos[0]-14<0 or pos[0]+14>img_size or pos[1]-14<0 or pos[1]+14>img_size:
            continue
        img[pos[0]-14:pos[0]+14, pos[1]-14:pos[1]+14, :] -= temp
        
    plt.imshow(img)
    plt.title(num)
    #plt.savefig("./results/tsne%d.png"%num)
    #plt.savefig("./results/tsne%d.eps"%num)
plt.show()

f:id:ksknw:20160430174653p:plain

元サイズ画像

model_faxocr = TSNE(n_components=2)
tsned_faxocr = model_faxocr.fit_transform(faxocr_data.reshape(len(faxocr_data),-1))

plt.figure(figsize=(30,30))
plt.scatter(tsned_faxocr[:,0], tsned_faxocr[:,1], c=faxocr_label, linewidths=0)
plt.show()

f:id:ksknw:20160430174724p:plain

あ、これいけるわ。

再びCNN

いけそうなのでCNNに突っ込んだ。結果は以下。

fig = plt.figure(figsize=[20,10])
plt.subplot(121)
accuracy = pd.read_csv("./results/mnist_trim_28_accuracy_row.txt", sep="\t")
loss = pd.read_csv("./results/mnist_trim_28_loss_row.txt", sep="\t")
print accuracy
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot()
plt.title("MNIST trim -> FaxOCR trim")
plt.subplot(122)
accuracy = pd.read_csv("./results/mnist_28_accuracy_row.txt", sep="\t")
loss = pd.read_csv("./results/mnist_28_loss_row.txt", sep="\t")
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot()
plt.title("MNIST->FaxOCR")
plt.show()
        epoch  test_accuracy
    0       1       0.787149
    1       2       0.867470
    2       3       0.895582
    3       4       0.947791
    4       5       0.931727
    5       6       0.923695
    6       7       0.931727
    7       8       0.955823
    8       9       0.935743
    9      10       0.939759
    10     11       0.927711
    11     12       0.959839
    12     13       0.955823
    13     14       0.955823
    14     15       0.931727
    15     16       0.931727
    16     17       0.923695
    17     18       0.951807
    18     19       0.935743
    19     20       0.951807
    20     21       0.951807
    21     22       0.959839
    22     23       0.951807
    23     24       0.955823
    24     25       0.939759
    25     26       0.963855
    26     27       0.951807
    27     28       0.959839
    28     29       0.963855
    29     30       0.955823
    30     31       0.963855
    31     32       0.951807
    32     33       0.967871
    33     34       0.967871
    34     35       0.967871
    35     36       0.963855
    36     37       0.979920
    37     38       0.959839
    38     39       0.967871
    39     40       0.971888
    40     41       0.951807
    41     42       0.967871
    42     43       0.959839
    43     44       0.967871
    44     45       0.971888
    45     46       0.971888
    46     47       0.947791
    47     48       0.927711
    48     49       0.959839
    49     50       0.927711

MNISTだけで学習したCNNを使って、無事90%を超えるぐらいの精度は出すことができた。 右は最初にやった正規化していないデータ。 f:id:ksknw:20160430180707p:plain

おわりに

normalize大事という意識は今までもあったつもりだったけど、正直ここまでとは思ってなかった。 今回は特にnormalizeされたデータであるMNISTのデータを使って、normalizeされてないFaxOCRのデータを認識しようとしていたのが良くなかった。実際にFaxOCRのデータでFaxOCRのデータを予測すると、それなりに上手くいっていた。CNNがきちんと学習してくれていたんだと思う。

一方でt-SNEを、正規化していないFaxOCRのデータに対して適用すると、かなりまずいことになっていた。今までなんとなくt-SNEに突っ込んで学習できそうかどうか見るというのをよくやっていたけど、もう少し気をつけたほうが良さそう。まずはt-SNEの論文をちゃんと読もうと思った。

おまけ

正規化したFaxOCRのデータを使って正規化したFaxOCRのデータを当てにいった。結果は以下。

28x28

fig = plt.figure(figsize=[10,10])
accuracy = pd.read_csv("./results/faxocr_trim_28_accuracy_row.txt", sep="\t")
loss = pd.read_csv("./results/faxocr_trim_28_loss_row.txt", sep="\t")
print accuracy
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot()
plt.title("FaxOCR trim -> FaxOCR trim 28x28")
plt.show()
        epoch  test_accuracy
    0       1       0.975904
    1       2       0.979920
    2       3       0.983936
    3       4       0.975904
    4       5       0.983936
    5       6       0.975904
    6       7       0.991968
    7       8       0.971888
    8       9       0.971888
    9      10       0.975904
    10     11       0.979920
    11     12       0.971888
    12     13       0.979920
    13     14       0.987952
    14     15       0.983936
    15     16       0.987952
    16     17       0.979920
    17     18       0.987952
    18     19       0.979920
    19     20       0.979920
    20     21       0.983936
    21     22       0.983936
    22     23       0.971888
    23     24       0.979920
    24     25       0.983936
    25     26       0.987952
    26     27       0.983936
    27     28       0.987952
    28     29       0.975904
    29     30       0.987952
    30     31       0.983936
    31     32       0.979920
    32     33       0.979920
    33     34       0.979920
    34     35       0.967871
    35     36       0.975904
    36     37       0.971888
    37     38       0.983936
    38     39       0.983936
    39     40       0.963855
    40     41       0.979920
    41     42       0.967872
    42     43       0.967871
    43     44       0.983936
    44     45       0.983936
    45     46       0.987952
    46     47       0.983936
    47     48       0.987952
    48     49       0.975904
    49     50       0.983936

何気に記録更新だった。たまたま感あるけど。 f:id:ksknw:20160430180726p:plain

import cPickle as pickle
import chainer
from chainer import cuda
import chainer.functions as F

def forward(x_data, y_data):
    x, t = chainer.Variable(x_data), chainer.Variable(y_data)
    h = F.max_pooling_2d(F.relu(model.conv1(x)), 2)
    h = F.max_pooling_2d(F.relu(model.conv2(h)), 2)
    h = F.dropout(F.relu(model.l1(h)), train=False)
    y = model.l2(h)

    return y, t,F.accuracy(y,t)
    
with open("./results/faxocr_trim_28_model_cnn_row.pkl", 'rb') as i:
    model = pickle.load(i)
from learn import read_imgs
test_data, test_label = read_imgs("./data/faxocr/test/trim_28")

test_data = test_data.reshape((len(test_data), 1, 28, 28))
test_data = test_data.astype(np.float32)
test_label = test_label.astype(np.int32)

y,t,acc = forward(test_data, test_label)
plt_num = 1
plt.figure(figsize=(10,10))
for i,(temp_y,temp_t,temp_test_data) in enumerate(izip(y.data,t.data, test_data)):
    if np.argmax(temp_y)!=temp_t:
        print "No.%d 正解:%d 出力:%d (%s)"%(i,temp_t, np.argmax(temp_y),temp_y)
        plt.subplot(2,2,plt_num)
        plt_num+=1
        plt.imshow(temp_test_data.reshape(28,28), cmap=plt.cm.gray_r)
plt.show()
    No.133 正解:9 出力:3 ([ -93.78523254  -76.56691742  -77.28510284   15.78890038  -87.52314758
      -58.99074554 -176.7293396  -101.68743896  -67.08155823   14.66318226])
    No.205 正解:9 出力:3 ([-44.00131989 -30.13937569 -40.71391678  17.24973297 -49.49590683
     -39.89061356 -82.53543091 -53.00856781   3.72428441   3.67775631])
    No.223 正解:9 出力:5 ([-24.54581642 -29.07642365 -31.50553131 -19.87841034 -23.92304039
      22.30206299 -18.75444031  -1.30475843 -20.22147751 -21.81653404])
    No.239 正解:5 出力:6 ([-12.27904129 -23.67368317 -23.48480606 -15.10187721 -11.83357048
       2.34730673   4.19125843 -13.38466549  -3.61155844 -24.05160713])

f:id:ksknw:20160430180739p:plain

48x48

fig = plt.figure(figsize=[10,10])
accuracy = pd.read_csv("./results/faxocr_trim_48_accuracy_row.txt", sep="\t")
loss = pd.read_csv("./results/faxocr_trim_48_loss_row.txt", sep="\t")
print accuracy
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot()
plt.title("FaxOCR trim -> FaxOCR trim 48x48")
plt.show()
        epoch  test_accuracy
    0       1       0.967871
    1       2       0.967871
    2       3       0.971888
    3       4       0.971888
    4       5       0.959839
    5       6       0.975904
    6       7       0.967871
    7       8       0.967871
    8       9       0.975904
    9      10       0.971888
    10     11       0.975904
    11     12       0.979920
    12     13       0.963855
    13     14       0.975904
    14     15       0.975904
    15     16       0.971888
    16     17       0.967871
    17     18       0.963855
    18     19       0.967871
    19     20       0.975904
    20     21       0.967871
    21     22       0.963855
    22     23       0.975904
    23     24       0.971888
    24     25       0.967871
    25     26       0.967871
    26     27       0.971888
    27     28       0.971888
    28     29       0.963855
    29     30       0.963855
    30     31       0.975904
    31     32       0.967871
    32     33       0.963855
    33     34       0.979920
    34     35       0.971888
    35     36       0.967871
    36     37       0.979920
    37     38       0.979920
    38     39       0.967871
    39     40       0.975904
    40     41       0.975904
    41     42       0.975904
    42     43       0.975904
    43     44       0.979920
    44     45       0.975904
    45     46       0.979920
    46     47       0.979920
    47     48       0.975904
    48     49       0.975904
    49     50       0.975904

f:id:ksknw:20160430180751p:plain

import cPickle as pickle
import chainer
from chainer import cuda
import chainer.functions as F

def forward(x_data, y_data):
    x, t = chainer.Variable(x_data), chainer.Variable(y_data)
    h = F.max_pooling_2d(F.relu(model.conv1(x)), 2)
    h = F.max_pooling_2d(F.relu(model.conv2(h)), 2)
    h = F.dropout(F.relu(model.l1(h)), train=False)
    y = model.l2(h)

    return y, t,F.accuracy(y,t)
    
with open("./results/faxocr_trim_48_model_cnn_row.pkl", 'rb') as i:
    model = pickle.load(i)
from learn import read_imgs
test_data, test_label = read_imgs("./data/faxocr/test/trim_48")

test_data = test_data.reshape((len(test_data), 1, 48, 48))
test_data = test_data.astype(np.float32)
test_label = test_label.astype(np.int32)

y,t,acc = forward(test_data, test_label)

plt_num = 1
plt.figure(figsize=(10,10))
for i,(temp_y,temp_t,temp_test_data) in enumerate(izip(y.data,t.data, test_data)):
    if np.argmax(temp_y)!=temp_t:
        print "No.%d 正解:%d 出力:%d (%s)"%(i,temp_t, np.argmax(temp_y),temp_y)
        plt.subplot(3,2,plt_num)
        plt_num+=1
        plt.imshow(temp_test_data.reshape(48,48), cmap=plt.cm.gray_r)
plt.show()
No.31 正解:7 出力:3 ([-37.6387825  -39.16486359 -42.47499466  16.941082   -35.50929642
 -28.37680244 -50.37080765 -18.17589951 -43.64836502 -26.19575882])
No.185 正解:3 出力:7 ([-21.7053051  -26.57409286 -22.67599106   5.41281748 -29.81308937
 -15.2016449  -49.91247177  11.7005167  -18.8066597  -14.08572674])
No.223 正解:9 出力:5 ([-26.99477196 -24.89096451 -43.02868652 -15.52783775 -25.02332115
  22.98989105 -17.69327545 -12.22776604 -12.18619633 -25.91065407])
No.228 正解:8 出力:9 ([ -4.70918608 -52.99452972 -28.40055275 -29.08706474   6.42821693
 -24.2097435  -52.58614731 -46.19430542 -14.41009426  12.67389202])
No.229 正解:2 出力:8 ([-22.11865997 -47.65965271  -1.94437969 -54.52325821 -24.92860031
 -19.47817802 -17.64678574 -40.01361084  26.35980225 -27.82418633])
No.239 正解:5 出力:8 ([ -3.87564874 -33.8475914  -27.65610313 -16.90502548 -28.28848457
 -23.67333221  -2.17319727 -35.55740356  22.8358345  -31.50007057])

f:id:ksknw:20160430180804p:plain

間違えちゃいけないデータを間違えているような気もするけど、前処理が適当で画像ぼけてるのがあれかも。 あとMNISTのデータとFaxOCRのデータをくっつけると汎化性能とか上がっていい感じかも。

あとは、データのaugmentationをもうちょっとちゃんとやるとか、複数のCNNでアンサンブル的なやつとかと思うけど、そのへんはよくわかってない

参考

FaxOCR手書き数字データの認識 その1

概要

FaxOCRという手書き数字認識の問題に挑戦した。 mnistで学習させたCNNでテストデータを判別すると正答率70%弱と低かった。 FaxOCRのデータだけで学習させたCNNでは96%程度の正答率だった。 mnistのデータとFaxOCRのデータはどうも違うようだけど、何が違うのかよくわからない。

以下はipython notebookの出力をちょこちょこいじったので、変なところがいくつかある。

はじめに

ツイッターを眺めているとこんなツイートを見つけた。

どうもFaxで送られてきた手書き数字を認識したいらしい。 はじめから電子データでいいんちゃうかとか、お役所も色々大変なんだろうなぁとか思いつつ。 手書き数字認識とかCNNに突っ込んだら終わりっしょぐらいの感じで始めた。

FaxOCR

サイトはこちら

バイナリ形式で画像とラベルのセットが用意されている。mnistと同じ形式らしい。 やり始めた当時は学習データが1711画像だった(たぶん)。 mnistが70000枚とかなのを考えても、CNNに入れるにはだいぶ少ないなぁという印象だったので、mnistで学習させたCNNを使ってFaxOCRのテストデータを分類することにした。

mnistで学習させたCNN

mnistでCNNを学習させるというのはTensorFlowのチュートリアルにあるぐらい鉄板なので、ググるとたくさんヒットする。個人的にchainerの使い方なら、なんとなく習得しているので、こちら のchainerの実装を使わせてもらうことにした。 CNNやらのコードは完全にコピペなのでここには書かないが、エポック数だけ20から50に変更した。

学習結果はtest_accuracy=0.694779とむっちゃ低かった。 以下はエポック毎のaccuracy(青)と誤差関数(緑)。学習は進んでいるのに、accuracyが上がっていない。

%matplotlib inline
import pylab as plt
import pandas as pd

accuracy = pd.read_csv("./accuracy_mnist2fax.txt", sep="\t")
loss = pd.read_csv("./loss_mnist2fax.txt", sep="\t")
fig = plt.figure(figsize=[10,10])

accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot(secondary_y=True)
plt.title("mnist->faxor")
plt.show()

f:id:ksknw:20160424232539p:plain

mnistのデータをテストデータにすると以下のようになる。 lossは同じぐらいなのに、accuracyは1エポックの時点で0.98を超えている。

accuracy = pd.read_csv("./accuracy_mnist2mnist.txt", sep="\t")
loss = pd.read_csv("./loss_mnist2mnist.txt", sep="\t")
fig = plt.figure(figsize=[10,10])
accuracy["test_accuracy"].plot()
plt.ylim([0,1])#plt.ylim([0,1])
loss["train_loss"].plot(secondary_y=True)
plt.title("mnist->mnist")
plt.show()

f:id:ksknw:20160424232639p:plain

一応、用意された学習データで学習すると以下のようになって、やっぱりデータ足りてない感がある。

accuracy = pd.read_csv("./accuracy_fax2fax.txt", sep="\t")
loss = pd.read_csv("./loss_fax2fax.txt", sep="\t")
fig = plt.figure(figsize=[10,10])
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot(secondary_y=True)
plt.title("faxor->faxor")
plt.show()

f:id:ksknw:20160424232651p:plain

データを眺める

手書き数字といえばmnistでおっけーというイメージだったので、ちょっとショックだった。 精度が出ていない原因として、データが質的に異なっているという可能性があるので、色々と可視化してみることにした。

画像データ

何はともあれ学習データの画像を見る。以下は可視化のコード。全部表示すると多すぎるので、適当に10件ずつだけ。

%matplotlib inline
from read_data import read,show
from sklearn.datasets import fetch_mldata
import numpy as np

mnist = fetch_mldata('MNIST original', data_home=".")
X = mnist.data
y = mnist.target
X = X.astype(np.float32)
y = y.astype(np.int32)

X /= X.max()
X_train = X
y_train = y

data = read()

test_data = read(dataset="testing")
X_test = test_data[0]
y_test = test_data[1]
X_test = X_test.reshape(len(y_test), -1) 
X_test = X_test / float(X_test.max())

import random
indecies = random.sample(range(len(X_train)), 1000)

for i in range(10):
    show(X_test[i].reshape(28,28))
    
print "###################################################"
print "############## mnist ここから######################"
print "###################################################"
for i in range(10):
    show(X_train[indecies[i]].reshape(28,28))

f:id:ksknw:20160424232712p:plain

f:id:ksknw:20160424232715p:plain

f:id:ksknw:20160424232728p:plain

f:id:ksknw:20160424232743p:plain

f:id:ksknw:20160424232753p:plain

f:id:ksknw:20160424232803p:plain

f:id:ksknw:20160424232808p:plain

f:id:ksknw:20160424232813p:plain

f:id:ksknw:20160424232819p:plain

f:id:ksknw:20160424232827p:plain


mnist ここから


f:id:ksknw:20160424232847p:plain

f:id:ksknw:20160424232852p:plain

f:id:ksknw:20160424232903p:plain

f:id:ksknw:20160424232932p:plain

f:id:ksknw:20160424232927p:plain

f:id:ksknw:20160424232908p:plain

f:id:ksknw:20160424232912p:plain

f:id:ksknw:20160424232918p:plain

f:id:ksknw:20160424232921p:plain

f:id:ksknw:20160424232924p:plain

ぱっと見た感じ、mnistのほうが太い線で書かれたものが多い気がする。とはいえ、細い線のものもあるし、認識してくれてもいいんじゃないかという印象。

t-SNEによる可視化

なんかよくわからんときは、とりあえずt-SNEにぶちこむというのが最近のマイブーム。 このあたりが詳しい。 これとかをみると、PCAよりええんちゃうって思う。 scikit-learnに関数があるので、使うのはとても簡単。mnist全データを突っ込むとメモリが足りないと怒られたので、ランダムに1000点選んで描画した。

data = np.r_[X_train[indecies], X_test]

from sklearn.manifold import TSNE
model = TSNE(n_components=2)

tsned = model.fit_transform(data)
import pylab as plt
label = np.r_[["b" for i in X_train[:1000]], ["r" for i in X_test]]
plt.figure(figsize=(30,30))
plt.scatter(tsned[:,0], tsned[:,1], c=label, linewidths=0)
plt.show()

f:id:ksknw:20160424232933p:plain

青(mnist)と赤(FaxOCR)のデータが明らかに分離している。こりゃあかんわって感じ。

データを増やしてCNN

そうこうしているうちに、こちらに精度を抜かれてしまっていた。

"適当に拡大縮小や回転をして画像データの枚数を11倍に(1枚から10枚生成)しました。"

とあって、すごく妥当だと思うし、なんで自分はこんなわけわからんことやってんだろうと思う。 とはいえ、精度で負けてるのはなんか悔しいので48x48のデータを51倍にしてCNNに突っ込んだ。

一応画像を適当に増やすところのコードは以下。

size = tuple(np.array([X_train[0].shape[1], X_train[0].shape[0]]))

new_imgs = []
new_labels = []
import random
import cv2
from itertools import izip
for img, label in izip(X_train, y_train):
    for i in range(50):
        rad = (random.random() - 0.5) * 0.5
        pos1 = (random.random() - 0.5) * 5
        pos2 = (random.random() - 0.5) * 5
        mat = np.float32([[np.cos(rad), -1 * np.sin(rad), pos1],
                          [np.sin(rad), np.cos(rad), pos2]])
        dst = cv2.warpAffine(img, mat, size, flags=cv2.INTER_LINEAR)
        new_imgs.append(dst)
        new_labels.append(label)
X_train = np.r_[X_train, new_imgs]
y_train = np.r_[y_train, new_labels]

accuracy(青)と誤差関数(緑)は以下。とりあえず96%とかいっているのでまあまあというところ。

accuracy = pd.read_csv("./accuracy_fax2fax_copied_48.txt", sep="\t")
loss = pd.read_csv("./loss_fax2fax_copied_48.txt", sep="\t")
fig = plt.figure(figsize=[10,10])
accuracy["test_accuracy"].plot()
plt.ylim([0,1])
loss["train_loss"].plot(secondary_y=True)
plt.title("faxor48->faxor48")
plt.show()

f:id:ksknw:20160424234032p:plain

ちなみに間違っていた画像は以下。これはしょうがないんじゃないかと思うものが多い。(というか学習データのほうは大丈夫なんだろうか…)

%matplotlib inline

import cPickle as pickle
import numpy as np
import chainer
from chainer import cuda
import chainer.functions as F

xp=np


def forward(x_data, y_data):
    x, t = chainer.Variable(x_data), chainer.Variable(y_data)
    h = F.max_pooling_2d(F.relu(model.conv1(x)), 2)
    h = F.max_pooling_2d(F.relu(model.conv2(h)), 2)
    h = F.dropout(F.relu(model.l1(h)), train=False)
    y = model.l2(h)

    return y, t,F.accuracy(y,t)
    
with open("model_cnn_48.pkl", 'rb') as i:
    model = pickle.load(i)
    
from read_data import read, show
test_data = read(dataset="testing", size=48)
X_test = test_data[0].astype(xp.float32)
y_test = test_data[1].astype(xp.int32)
X_test = X_test.reshape(len(y_test), -1)
X_test = X_test / float(X_test.max())
X_test = X_test.reshape((len(X_test), 1, 48, 48))

y,t,acc = forward(X_test, y_test)
print "#################################"
print "accuracy: " + str(acc.data)
print "#################################"

from itertools import izip

for i,(temp_y,temp_t,temp_X_test) in enumerate(izip(y.data,t.data, X_test)):
    if np.argmax(temp_y)!=temp_t:
        print "No.%d 正解:%d 出力:%d (%s)"%(i,temp_t, np.argmax(temp_y),temp_y)
        show(temp_X_test.reshape(48,48))
    accuracy: 0.967871487141
    No.14 正解:9 出力:8 ([  3.53644633 -66.54345703 -44.21648026 -30.54708862 -47.66264725 -61.4640274  -36.32198715 -54.34354401  28.37440872 14.29025173])

f:id:ksknw:20160424234123p:plain

No.116 正解:9 出力:5 ([-35.68978119 -23.23931885 -76.64511108 -26.64510727 -42.89046097
      50.29698944 -16.02383232 -52.14741135 -22.0790844  -29.76072311])

f:id:ksknw:20160424234135p:plain

    No.127 正解:9 出力:8 ([ -0.36938047 -72.61532593 -38.76506424 -25.25551796 -40.20273972
     -45.48267365 -41.19197464 -39.32862854  26.44895935  11.87366581])

f:id:ksknw:20160424234144p:plain

    No.170 正解:7 出力:1 ([-13.41542912   8.12284565 -28.43426895  -6.35744619 -42.47330093
     -31.51161766 -36.15800858  -0.18305674 -23.2899704  -10.80175686])

f:id:ksknw:20160424234153p:plain

    No.172 正解:9 出力:1 ([-18.31829453   8.41508961 -24.00697708 -18.74674988 -15.38268757
     -19.12284851 -21.66376686 -29.5259037   -9.97081757   5.56778383])

f:id:ksknw:20160424234202p:plain

    No.196 正解:9 出力:4 ([-25.27557564 -20.27404976 -33.58036041 -38.26721573  18.712677
     -17.12530899 -26.8935318  -34.70022964 -20.10196686  11.45114803])

f:id:ksknw:20160424234210p:plain

   No.231 正解:5 出力:8 ([ -3.76606822 -27.21845436 -40.48172379 -41.74909592 -20.90390015
     -27.25065613   8.70675373 -39.89279938  28.23112106 -19.23669815])

f:id:ksknw:20160424234219p:plain

    No.247 正解:2 出力:4 ([-15.3788166   -8.5814476  -14.29022121 -16.61594963   5.68172646
      -8.83567333   0.14922404 -22.43881035 -17.46813393 -14.75987816])

f:id:ksknw:20160424234228p:plain

mnistデータとFaxOCRデータの違い

データを公開された方の本来の目的からすると、パラメータチューニングとかして性能を上げたほうがいいのかもしれないけど、正直そっちにはあまり興味がない。Kaggleガチ勢の方がこういうの とか出してくれているので、参考にするといいのかもしれない。

個人的に気になったのは今回のデータとmnistの違い。 見た目は同じような手書き文字なのに、t-sneで可視化すると明らかに分離している。 Faxのデータはどうも線の細さを統一したり、回転させたりと前処理を結構しているらしいので、それが影響しているのかなと思った。 なので、生データを可視化してみる。FaxOCRの画像サイズはまちまちだったので、ImageMagickで28x28にリサイズした。アスペクト比は保存していないのでややまずい。

import cv2 as cv
import glob
mnist = fetch_mldata('MNIST original', data_home=".")
X = mnist.data
y = mnist.target
X = X.astype(np.float32)
y = y.astype(np.int32)

X /= X.max()
X_train = X
y_train = y

X_test = []
y_test = []
for img_file in glob.glob("./data/mustread/28/*.png"):
    y_test.append(int(img_file[16]))
    X_test.append(255 - cv.imread(img_file, flags=0))
    
data_row = np.r_[X_train[indecies], np.array(X_test).reshape(len(X_test),-1)]
from sklearn.manifold import TSNE
model = TSNE(n_components=2)

tsned_row = model.fit_transform(data_row)
label = np.r_[["b" for i in range(1000)], ["r" for i in X_test]]
plt.figure(figsize=(30,30))
plt.scatter(tsned_row[:,0], tsned_row[:,1], c=label, linewidths=0)
plt.show()

f:id:ksknw:20160424234255p:plain

なんでや…

ミス訂正(2016/4/26)

FaxOCRのほうが0~1に補正されていなかった。正しいt-SNEの結果は以下。 f:id:ksknw:20160426233426p:plain これなら、線が細いやつが固まっていると思えば、まだありえる(?)

一旦終わり

正直なにがダメなのかよくわかってない。というかバグなんじゃないのか、僕のコードがなにかやらかしてるんじゃないのか。 見た目おんなじに見えるんだけど、何か本当に違うのか… 誰かバグとか根本的な間違いとか見つけたら教えてください…

ミス訂正(2016/4/26)

コードにミスがあったので、元画像データはmnistの一部がある空間にありそうだとわかった。時間があるときにもうちょっとちゃんとやって追記します。

続き

ksknw.hatenablog.com

参考

omake+latexdiff-gitで作る最強のTeX環境

この記事のやつはもうやってない

最近はlatexmkで自動コンパイルして,査読をお願いするときだけlatexdiffを使って差分ファイルを作っている

デモ

これがぼくのかんがえたさいきょうのてふかんきょうだ。

概要

omakeでディレクトリの変更を監視して、TeXファイルを自動コンパイルする。 コンパイル時にlatexdiff-gitで前回のコミットとの差分を出力し、それもコンパイルする。

はじめに

最近周りでTeXを書いている人があまりいない。Word派の主張を聞いているとどうもコンパイルがイケてないらしい。 実際、TeXを書いているときにいちいちmakeとか叩くのは正直だるい。 (コマンドとか無理という主張はどうしようもないので無視した)

探しているとこんな最高の記事を見つけた。 導入するついでに以前書いたものを合体させて、 自動で本体と差分のpdfを出力するようにする。

  • OMake 0.9.8.5
  • latexdiff 1.0.2

omake

omakeについては、こちらを参考にした。 omakeはなんかよくわからんけど、makeの代わりに使えるもの。 まだ全く調べてないけど、今回使ったディレクトリ監視だけでなく他にも色々できそう。

aptでインストールする。 famというのはディレクトリの監視に必要らしい。

 $ sudo apt-get install omake fam

適当なTeXファイルのあるディレクトリで、

 $ omake --install

とすると、OMakefileとOmakerootというファイルが生成される。

latexdiff

こちらを参考にして日本語に対応する必要がある。

latexdiffについては以前こんなものを書いた。 以前の記事では以下のようなコマンドを叩いて差分ファイルを生成していたが、

git ldiff HEAD~1 > diff.tex

色々いじっているうちに、 latexdiff-gitというコマンドが用意されていることに気づいた。 このコマンドを使うと以下のように簡単に前回のコミットとの差分をTeX形式で出力できる。

latexdiff-git -r -e utf8 test.tex

合体させる

こちらこちらの記事にあるOMakefileにlatexdiffの部分を追記して、以下のようなOmakefileを作った。

TARGET = test
DIFF = test-diff
IMAGE_DIR = figs

LATEX = platex
BIBTEX = pbibtex
DVIPDFM = dvipdfmx
LATEXDIFF = latexdiff-git

DIFFFLAGS = -r -e utf8 --force

# Bounding Box生成コマンド
EBB = extractbb

# グロブ展開に失敗したときに空の文字列を返すようにする
GLOB_OPTIONS = i

# Bounding Boxの生成ルール
.SUBDIRS: $(IMAGE_DIR)
    %.xbb: %.png
        $(EBB) $<
    %.xbb: %.jpg
        $(EBB) $<
    %.xbb: %.pdf
        $(EBB) $<

TEX_FILES = $(glob *.tex */*.tex)
BIB_FILES = $(glob *.bib)
EPS_IMAGE_FILES = $(glob $(IMAGE_DIR)/*.eps)
OTHER_IMAGE_FILES = $(glob $(IMAGE_DIR)/*.png $(IMAGE_DIR)/*.jpg $(IMAGE_DIR)/*.pdf)
IMAGE_FILES = $(EPS_IMAGE_FILES) $(OTHER_IMAGE_FILES)
XBB_FILES = $(addsuffix .xbb, $(removesuffix $(OTHER_IMAGE_FILES)))

# コンパイルに必要なファイル
TEXDEPS[] = $(TEX_FILES) $(BIB_FILES) $(IMAGE_FILES) $(XBB_FILES)

LaTeXDocument($(TARGET), $(TARGET))
LaTeXDocument($(DIFF), $(DIFF))

.DEFAULT: $(TARGET).pdf $(DIFF).pdf

$(DIFF).tex:$(TARGET).tex
  $(LATEXDIFF) $(DIFFFLAGS) $<

.PHONY: clean
clean:
    rm $(glob *.toc *.log *.pdf *.dvi *.fls *.aux *.maf *.mtc *.bbl *.blg) $(XBB_FILES)

Emacsユーザーなら

after-save-hook にmakeを設定するだけでもいいじゃないかという気もする。

参考

raspberry pi2 + カメラ付きサーボで特定の色を追っかけるカメラを作る

概要

raspberry pi2 を使ってカメラ付きサーボモータを制御した。 python+OpenCVで特定の色の中心を追っかけるようにした。 かくかくしてるけどだいたい良い。

はじめに

学生の頃、ロボットに触る機会はよくあったが、だいたい「これが見えたらあっちに歩け」ぐらいの抽象度で、モータの角度を何度にするとか、電気回路がどうとかはさっぱりわからない。 仮にも機械工学科目を修めておいてこれではいけないなと思って、とりあえず部屋に転がっていたサーボモータを動かしてみることにした。

使用したもの

  • 秋月のwebカメラ付きサーボ
  • raspberry pi2
  • ブレッドボードなど

画像中の特定の色の中心を求める

OpenCVを使って何も考えずに書くと以下のようになる。 今回は青色を追っかけるようにした。

import cv2
import numpy as np

capture = cv2.VideoCapture(0)
capture.set(3, 320)
capture.set(4, 240)

dst = np.zeros((320, 240), dtype=np.uint8)

while True:
    _, img = capture.read()
    dst = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    conditions = (dst[:,:,0]<150)*1 * ((dst[:,:,0]>90)*1) * (dst[:,:,1]>100)

    sum_x = 0
    sum_y = 0
    sum_item = 0
    for x,row in enumerate(conditions):
        for y,item in enumerate(row):
            sum_x += x * item
            sum_y += y * item
            sum_item += item
    mean_x = sum_x / sum_item
    mean_y = sum_y / sum_item

    cv2.circle(img, (mean_y, mean_x),10, (0,0,255), -1)
    cv2.imshow("camera", img)
    cv2.waitKey(1)

capture.release()
cv2.destroyAllWindows()

pythonを使っている人にはわかると思うが、このプログラムはとんでもなく遅い。 実際にラズパイで実行してみると体感2fpsぐらいだった。 というわけで、以下のように書き換えた。

import cv2
import numpy as np

capture = cv2.VideoCapture(0)
capture.set(3, 320)
capture.set(4, 240)

while True:
    _, img = capture.read()
    dst = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    conditions = (dst[:,:,0]<150)*1 * ((dst[:,:,0]>90)*1) * (dst[:,:,1]>100)

    sum_item = np.sum(conditions)
    t_conditions = np.transpose(conditions)

    temp = range(len(conditions))
    mean_x = np.sum([temp * t_condition for t_condition in t_conditions])/sum_item
    temp = range(len(t_conditions))
    mean_y = np.sum([temp * condition for condition in conditions])/sum_item

    cv2.circle(img, (mean_y, mean_x),10, (0,0,255), -1)
    cv2.imshow("camera", img)
    cv2.waitKey(1)

capture.release()
cv2.destroyAllWindows()

これでまあ満足できる程度の早さになる。 実行するとこのようになる。

f:id:ksknw:20160228102336p:plain

サーボモータを制御する

回路を適当に組む。 サーボの動作電圧が4.8Vだったので、電池を3つつなげて4.5Vで使用した。 サーボの信号線をラズパイのPWMができるピンに接続する。

f:id:ksknw:20160228003450j:plain

ラズパイのピンのレイアウトはこちらを参考にした。 PWMで制御するために今回はGPIO13と18を使う。

こちら を参考にしてモータを動かしてみる。

import wiringpi2

PWM_PIN = 13
#PWM_PIN = 18
DUTY_MAX = 123 # 90°
DUTY_MIN = 26  # -90°
DUTY_HOME = 74 # 0°
duty = 0

wiringpi2.wiringPiSetupGpio()
wiringpi2.pinMode(PWM_PIN, wiringpi2.GPIO.PWM_OUTPUT)
wiringpi2.pwmSetMode(wiringpi2.GPIO.PWM_MODE_MS)
wiringpi2.pwmSetClock(375)

wiringpi2.pwmWrite(PWM_PIN, DUTY_HOME)
wiringpi2.delay(100)

def move(degree):
    duty = int((DUTY_MAX-DUTY_MIN)/180.0 * degree + DUTY_HOME)
    wiringpi2.pwmWrite(PWM_PIN, duty)

for degree in [0,45,90,0,-45,-90,0,30,-30,0]:
    move(degree)
    wiringpi2.delay(500)

実行するとこんな感じで動く。

f:id:ksknw:20160228103343g:plain

青色を追っかける

2つのプログラムを組み合わせる。

import wiringpi2
import cv2
import numpy as np

PIN_X = 13
PIN_Y = 18
DUTY_MAX = 123 # 90°
DUTY_MIN = 26  # -90°
DUTY_HOME = 74 # 0°
duty = 0

wiringpi2.wiringPiSetupGpio()
wiringpi2.pinMode(PIN_X, wiringpi2.GPIO.PWM_OUTPUT)
wiringpi2.pinMode(PIN_Y, wiringpi2.GPIO.PWM_OUTPUT)
wiringpi2.pwmSetMode(wiringpi2.GPIO.PWM_MODE_MS)
wiringpi2.pwmSetClock(375)

wiringpi2.pwmWrite(PIN_X, DUTY_HOME)
wiringpi2.pwmWrite(PIN_Y, DUTY_HOME)
wiringpi2.delay(100)

def move(degree_x, degree_y):
    duty_x = int((DUTY_MAX-DUTY_MIN)/180.0 * degree_x + DUTY_HOME)
    duty_y = int((DUTY_MAX-DUTY_MIN)/180.0 * degree_y + DUTY_HOME)
    wiringpi2.pwmWrite(PIN_X, duty_x)
    wiringpi2.pwmWrite(PIN_Y, duty_y)
    wiringpi2.delay(100)

def detect_blue_center(img):
    dst = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    conditions = (dst[:,:,0]<150)*1 * ((dst[:,:,0]>90)*1) * (dst[:,:,1]>100)
    sum_item = np.sum(conditions)
    t_conditions = np.transpose(conditions)

    temp = range(len(conditions))
    mean_x = np.sum([temp * t_condition for t_condition in t_conditions])/sum_item
    temp = range(len(t_conditions))
    mean_y = np.sum([temp * condition for condition in conditions])/sum_item
    return mean_x, mean_y

capture = cv2.VideoCapture(0)
capture.set(3, 320)
capture.set(4, 240)

now_angle_x = 0
now_angle_y = 0

while True:
    _, img = capture.read()
    mean_x, mean_y = detect_blue_center(img)
    cv2.circle(img, (mean_y, mean_x),10, (0,0,255), -1)

    angle_x = (mean_x - 160)/160.0*7
    angle_y = -(mean_y - 120)/120.0*7
    move(now_angle_y + angle_y, now_angle_x + angle_x)
    now_angle_y = now_angle_y + angle_y
    now_angle_x = now_angle_x + angle_x

    cv2.imshow("camera", img)
    cv2.waitKey(1)

capture.release()
cv2.destroyAllWindows()

実行するとこんな感じ。 かくかくしてるし、ラグがあって微妙だけど動くからいいかな。 f:id:ksknw:20160228142524g:plain

おわりに

かくかくしてるのとかラグとか直すべきところはいっぱいあるけど、特に目標があって始めたわけでもないからやる気がでない。

追記(3月4日)

なんとなくdelayを入れていたけど、とってみるとカクカクしなくなった。金麦飲んでたせいで気づかなかったわー金麦のせいだわー。 ついでに画像サイズ小さくしたらかなりよくなった。

参考

Ubuntu14.04での自分的環境構築をちゃんと書く

概要

先日なんとなくgnome3を使ってみようと思い、色々やっているうちにデスクトップが死んだ。 復旧するより入れなおしたほうが早い感じだったので、昔書いたものを見ながらインストールしようと思ったが、あまりに適当なことしか書いておらず絶望した。 今後のために、クリーンインストールされたubuntu14.04から、ちゃんと自分の設定を復元できるように書いておく。

インストールするものは以下とその他細かい諸々。

  • cuda
  • chainer
  • TensorFlow
  • Emacs

細かいこと色々

$ LANG=C xdg-user-dirs-gtk-update
$ sudo apt-get update
$ sudo apt-get upgrade
$ sudo apt-get install freeglut3-dev build-essential libx11-dev libxmu-dev  libxi-dev libgl1-mesa-glx libglu1-mesa libglu1-mesa-dev linux-headers-`uname -r`
$ sudo apt-get install zsh
$ chsh -s /usr/bin/zsh
$ sudo apt-get install git
$ sudo apt-get install gimp ipython ipython-notebook

cudaのインストール

nvidiaのサイトから自分にあったグラフィックボードのドライバとcudaをダウンロードする。 TensorFlowは7.0をインストールしろと言っているので、最新版の7.5ではなく、7.0をダウンロードする。

まだよくわかってないけど、こちら にやれと書いてあるので、ブラックリストに登録をする。たぶんドライバが干渉するのを防ぐんだと思う。

$ sudo nano /etc/modprobe.d/blacklist-nouveau.conf

ファイルを開いて以下を記述。

blacklist nouveau
blacklist lbm-nouveau
options nouveau modeset=0
alias nouveau off
alias lbm-nouveau off

その後色々

$ echo options nouveau modeset=0 | sudo tee -a /etc/modprobe.d/nouveau-kms.conf
$ update-initramfs -u
$ sudo update-initramfs -u
$ sudo reboot

再起動したらCUIに切り替えて以下を実行する。

$ sudo service lightdm stop
$ chmod +x NVIDIA-Linux-x86_64-352.63.run
$ sudo ./NVIDIA-Linux-x86_64-352.63.run
$ chmod +x cuda_7.0.28_linux.run
$ sudo ./cuda_7.0.28_linux.run

cudaの方でグラフィックドライバをインストールするかと聞かれるが、別途入れているのでnoを選択する。 ちなみにyesを選択すると、黒画面に白カーソルが表示されるだけで操作できず、しかもCUIにも入れないという悲しい事態になった。

$ sudo service lightdm start

ここで正しくGUIが表示されなかったら諦める。

.zshrcに以下を追記して、CUDAにパスを通す。

export CUDA_HOME=/usr/local/cuda-7.0
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${CUDA_HOME}/lib64
export PATH=$PATH:${CUDA_HOME}/bin

chainerのインストール

とりあえずpipを入れる。最近はget-pip.pyを実行して入れるらしい。 こちらを参考にした。

chainerを入れるときにhdf5がなんとかかんとかと言われたが、こちら を見てごにょごにょするとうまくいった。 実行したのは以下。

$ curl -kL https://raw.github.com/pypa/pip/master/contrib/get-pip.py | sudo python
$ sudo apt-get install python-dev
$ sudo pip install -U cython
$ sudo apt-get install libblas-dev liblapack-dev gfortran
$ sudo apt-get install g++
$ sudo pip install  numpy scipy
$ sudo apt-get install libhdf5-dev
$ sudo pip install chainer

ついでに使いそうなものを色々入れる

$ sudo pip install sklearn
$ sudo pip install pandas
$ sudo pip install xgboost

Emacsのインストール

とりあえずgit hubから自分の設定ファイルを落としてくる。

$ ssh-keygen

できた id_rsa.pub をgithubの自分のレポジトリに登録して、以下を実行。

$ git clone git@github.com:kskkwn/.emacs.d.git

依存しているものを色々入れる。

$ sudo apt-get install cmigemo pyflakes ruby
$ sudo pip install jedi epc autopep8

Emacsをインストールする。24.4用に設定ファイルを作っているので、面倒なことを避けるためにバージョンは固定。

$ sudo apt-get install build-essential
$ sudo apt-get build-dep emacs24
$ wget http://ftp.jaist.ac.jp/pub/GNU/emacs/emacs-24.4.tar.gz
$ tar -xf emacs-24.4.tar.gz
$ ./configure
$ make
$ sudo make install

ちゃんとできていればここまでで、いつものEmacsが起動するはず。 前に書いたやつからjedi入れたりtabbar入れたりhelm入れたりして、だいぶ変わっているし、またなんか書こう。

TensorFlow

基本的に公式サイト の言うとおりに実行すればいい。

GPU使いたいので、まずはcuDNNをインストールする。 cuDNNはNVIDIAのDeveloperサイトに登録しないとダウンロードできない。 ユーザー登録の際に色々聞かれるけど、「Tensor Flow使いたいねん。しゃす。」ぐらいのことを英語で書いて、数時間待ったらあっさり登録できた。これも最新版でないものをダウンロードする。

$ tar xvzf cudnn-6.5-linux-x64-v2.tgz
$ sudo cp cudnn-6.5-linux-x64-v2/cudnn.h /usr/local/cuda-7.0/include
$ sudo cp cudnn-6.5-linux-x64-v2/libcudnn* /usr/local/cuda-7.0/lib64

続いてTensorFlowのインストール。

$ sudo pip install --upgrade https://www.tensorflow.org/versions/master/get_started/os_setup.htmlhttps://storage.googleapis.com/tensorflow/linux/gpu/tensorflow-0.6.0-cp27-none-linux_x86_64.whl

やってみるとこんな感じのメッセージが表示された。

>>> import tensorflow
I tensorflow/stream_executor/dso_loader.cc:101] successfully opened CUDA library libcublas.so.7.0 locally
I tensorflow/stream_executor/dso_loader.cc:101] successfully opened CUDA library libcudnn.so.6.5 locally
I tensorflow/stream_executor/dso_loader.cc:101] successfully opened CUDA library libcufft.so.7.0 locally
I tensorflow/stream_executor/dso_loader.cc:101] successfully opened CUDA library libcuda.so locally
I tensorflow/stream_executor/dso_loader.cc:101] successfully opened CUDA library libcurand.so.7.0 locally

successfullyって言ってるし、大丈夫やろう。

今後の課題

  • pythonの仮想環境の構築 正直みんな何のためにやっているのかよくわからないので、今はやってないけど、そのうちやったほうがいいのかもしれない。
  • TeX環境の構築 家でTeX書くことないしなぁと思っていれていない。
  • C++環境の構築 ここ半年ぐらいC++書いてない。できれば今後も書きたくない。

参考

python_tips_備忘録

はじめに

コードを書いていると、「あ、これ前もググって実装したやつだ」ということがよくある。 しかもどうやったか全く覚えていない。悲しい。 あげく検索ワードすらわからないときが頻繁にある。悲しい。

これ以上、そのような悲しいことを起こさないために、たまに使う細かい処理について書いていく。 間違っている部分や黒魔術的な部分もあるかもしれないので注意。 随時追加するつもり。 番号は参考サイトを示している。サンプルコードについて非常に参考にした。

for文が正しく終了したときに、処理を実行する。(for...else) [1]

正常終了するとelseというのが、なんとなく慣れない。いちいちフラグとかたてなくてもいい。

for i in range(10):
    print i
else:
    print "hoge"
0
1
2
3
4
5
6
7
8
9
hoge
for i in range(10):
    print i
    break
else:
    print "hoge"
0

コマンドライン引数を処理する。(argparser, click) [2]、[3]

コマンドライン引数を処理するモジュールとして、argparserとclickがある。まずはargpaeser。

import argparse

parser = argparse.ArgumentParser()
parser.add_argument('-n', '--n_iter', default=1000, type=int, help="num_iteration")
parser.add_argument('-g', '--gpu', action='store_true')
args = parser.parse_args()

print args
print args.n_iter
print args.gpu
$ python temp.py -n 10 --gpu
Namespace(gpu=True, n_iter=10)
10
True

続いてclick。最近知ったので、使い方がまだ良くわかっていないが、こちらのほうがコンパクトに書けて便利な気がする。

import click

@click.command()
@click.argument("name")
@click.option("-n", "--num_iter", type=int, help="num iteration")
def main(name, num_iter):
    for i in range(num_iter):
        msg = "hello %s %d" % (name, i)
        print msg

if __name__ == '__main__':
    main()
$ python temp2.py world -n 10
hello world 0
hello world 1
hello world 2
hello world 3
hello world 4
hello world 5
hello world 6
hello world 7
hello world 8
hello world 9

ファイル名に時間を埋め込む 。(datetime)[4]

結果の出力など、うっかり上書きしないようにするときによく使う。

from datetime import datetime
filename = "results/" + datetime.now().strftime("%Y%m%d_%H%M%S") + ".csv"
print filename
results/20151130_225340.csv

ファイルを開く。(with open... as f)

普通にopenだけで開くと後でcloseしないといけない。忘れるので、withを使う。 自作のクラスもこういうものに対応させておくと後で便利かも。

with open("temp.csv") as f:
    for i in f:
        print i

プログラムの進行状況などをprintするときに改行せず上書きする。 (sys.stdout.write)

ループの中で、print i などと書くと、いちいち改行されて鬱陶しい。 以下のコードでいい感じに進行状況を表示できる。

import sys
import time
for i in range(100):
    sys.stdout.write("\riter=%d" % i)
    sys.stdout.flush()
    time.sleep(0.1)
iter=99

ファイル名で検索してリストを取得する。(glob)

paramsディレクトリ内に存在する全てのパラメータファイルを読んでいって、順番に実行するとか、resultsディレクトリ以下の全てのファイルの平均値をとるとか、やりたいときに使う。

import glob

print glob.glob("./*")
print glob.glob("./*.py")
['./tips.org', './temp.csv', './Untitled.ipynb', './temp1.py', './temp2.py', './temp10.py']
['./temp1.py', './temp2.py', './temp10.py']

文字中に含まれる数字でソート。[5]

ファイル名などに数字が付いている時、普通にsortメソッドを使うと、2の前に10が来たりしてしまう。はじめから%05dとかで書いておけば不要。

import re
files = glob.glob("./*.py")
temp = [(re.search("[0-9]+", x).group(), x) for x in files]
temp.sort(cmp = lambda x, y: cmp(int(x[0]), int(y[0])))
print  [x[1] for x in temp]
['./temp1.py', './temp2.py', './temp10.py']

別のディレクトリのファイルをインポートする。(sys.path.append)

昔作ったコードを再利用したいとき、コピーしてきてもいいけど、これでもOK。

import sys
sys.path.append("../kaggle/mnist/")
import nn
net = nn.Perceptron(10,2)

ある条件を満たす数列を生成する。(リスト内包表記、ジェネレータ)

リスト内包表記をごちゃごちゃすると色々できる。 あまり使いすぎるとわかりにくいので、ほどほどに使う。

x = [i**2 for i in range(10) if i%2==0]
print x
[0, 4, 16, 36, 64]

ジェネレータでもいける。yieldは遅延評価なので、1つ値を生成するコストが大きいかつ、いくつ必要かわからないときとかいいかもしれない。

def generator():
    counter = 0
    while(True):
        yield counter**2
        counter += 2
for i in generator():
    print i
    if i>=64:
        break
0
4
16
36
64

学習データからランダムなミニバッチを作る。( np.random.permutation(N) ) [6]、[7]

SGDとかやるときに使う。

    def train(self, inputs, labels, tests):
        N = len(inputs)
        inputs = np.array(inputs).astype(np.float32)
        labels = np.array(labels).astype(np.int32)
        tests = np.array(tests).astype(np.float32)

        for epoch in range(self.num_epoch):
            print "epoch: %d" % epoch
            perm = np.random.permutation(N)
            for i in range(0, N, self.batchsize):
                x_batch = cuda.to_gpu(inputs[perm[i:i + self.batchsize]])
                y_batch = cuda.to_gpu(labels[perm[i:i + self.batchsize]])

KFold Cross Validation(sklearn.cross_validation.KFold) [8]

一瞬自分でつくろうかと思ったけど、探してみたらやっぱりあったscikit_learn。 自分で書くコードは最小限にする方針からこれを使う。

kf = cross_validation.KFold(len(all_vec), n_folds=10, shuffle=True)
for train_index, test_index in kf:
    test_vec = all_vec[test_index]
    test_label = labels[test_index]
    train_vec = all_vec[train_index]
    train_label = labels[train_index]

特徴ベクトル(の一部)を正規化する。(sklearn.preprocessing.MinMaxScaler) [9]

sklearnに正規化の関数があるので、それを使う。1ofKの変数とかは正規化しなくてもいいと思うので、一部だけやる例。

min_max_scaler = sklearn.preprocessing.MinMaxScaler()
no_normalize_index = 8
min_max_scaler.fit(train_vec[:, no_normalize_index:])
test_vec[:, no_normalize_index:] = min_max_scaler.transform(test_vec[:, no_normalize_index:])
train_vec[:, no_normalize_index:] = min_max_scaler.transform(train_vec[:, no_normalize_index:])

bool型の掛け算順

以前、一度バグを産んだ原因コード。演算子の優先度は、(符号としての) - > * > (引き算としての)- なので、こういうことをすると、一見引き算の順序を変えただけなのに、結果がかわる。これからはおとなしくastypeでキャストしようと心に誓った。

import numpy as np
bool_list = np.array([True, False])

print [0,0] - bool_list * 0.5
print -bool_list * (0.5) + [0,0]

print - bool_list.astype(int) * 0.5 + [0,0]
[-0.5  0. ]
[ 0.   0.5]
[-0.5  0. ]

クラスタリングなどの結果をscatterでプロットする(plt.scatter)

だいたい特徴量的にx, yを管理しているので、特徴ベクトル行列のrow, colを入れ替えてプロットする。 colorの変数に適当に突っ込むとよしなにやってくれる。 ついでにPCAもする。

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pcaed_data = pca.fit_transform(data[:,1:]/255.0)
pcaed_data = zip(*pcaed_data)
plt.scatter(pcaed_data[0], pcaed_data[1], marker=".",c=data[:,0], linewidths=0, alpha=0.7)

csvファイルからデータを読み込む。(pandas)

普通にopenでファイルを開いてもいいけど、pandasでやると欠損値とかの扱いがかなり楽。

import pandas as pd
data = pd.read_csv("train.csv")

data = data[~data.Cabin.isnull()] #"Cabin"がNanでないものだけを抽出する。
print data
     PassengerId  Survived  Pclass  \
1              2         1       1
3              4         1       1
6              7         0       1
10            11         1       3
..           ...       ...     ...
879          880         1       1
887          888         1       1
889          890         1       1

                                                  Name     Sex   Age  SibSp  \
1    Cumings, Mrs. John Bradley (Florence Briggs Th...  female  38.0      1
3         Futrelle, Mrs. Jacques Heath (Lily May Peel)  female  35.0      1
6                              McCarthy, Mr. Timothy J    male  54.0      0
10                     Sandstrom, Miss. Marguerite Rut  female   4.0      1
..                                                 ...     ...   ...    ...
879      Potter, Mrs. Thomas Jr (Lily Alexenia Wilson)  female  56.0      0
887                       Graham, Miss. Margaret Edith  female  19.0      0
889                              Behr, Mr. Karl Howell    male  26.0      0

     Parch       Ticket      Fare        Cabin Embarked
1        0     PC 17599   71.2833          C85        C
3        0       113803   53.1000         C123        S
6        0        17463   51.8625          E46        S
10       1      PP 9549   16.7000           G6        S
..     ...          ...       ...          ...      ...
879      1        11767   83.1583          C50        C
887      0       112053   30.0000          B42        S
889      0       111369   30.0000         C148        C

[204 rows x 12 columns]

リストで重複した要素を削除する。(pandas)

例えばIDが重複していて、どちらかの要素がいらないとき、これを見つけて削除する。 オプションで上にあるものを残すか、下にあるものを残すかを選べる。

print data[~data.Cabin.duplicated()]
     PassengerId  Survived  Pclass  \
1              2         1       1
3              4         1       1
6              7         0       1
10            11         1       3
..           ...       ...     ...
879          880         1       1
887          888         1       1
889          890         1       1

                                                  Name     Sex   Age  SibSp  \
1    Cumings, Mrs. John Bradley (Florence Briggs Th...  female  38.0      1
3         Futrelle, Mrs. Jacques Heath (Lily May Peel)  female  35.0      1
6                              McCarthy, Mr. Timothy J    male  54.0      0
10                     Sandstrom, Miss. Marguerite Rut  female   4.0      1
..                                                 ...     ...   ...    ...
879      Potter, Mrs. Thomas Jr (Lily Alexenia Wilson)  female  56.0      0
887                       Graham, Miss. Margaret Edith  female  19.0      0
889                              Behr, Mr. Karl Howell    male  26.0      0

     Parch       Ticket      Fare        Cabin Embarked
1        0     PC 17599   71.2833          C85        C
3        0       113803   53.1000         C123        S
6        0        17463   51.8625          E46        S
10       1      PP 9549   16.7000           G6        S
..     ...          ...       ...          ...      ...
879      1        11767   83.1583          C50        C
887      0       112053   30.0000          B42        S
889      0       111369   30.0000         C148        C

[147 rows x 12 columns]

特定の文字を含むものを取り出す(pandas)

数字は簡単だけど、文字の一部に含まれているというのがよくわからなかった。これもpandasを使えば簡単。

print data[data.Cabin.str.contains("A")]
     PassengerId  Survived  Pclass  \
23            24         1       1
96            97         0       1
174          175         0       1
185          186         0       1
209          210         1       1
284          285         0       1
445          446         1       1
475          476         0       1
556          557         1       1
583          584         0       1
599          600         1       1
630          631         1       1
647          648         1       1
806          807         0       1
867          868         0       1

                                                  Name     Sex  Age  SibSp  \
23                        Sloper, Mr. William Thompson    male   28      0
96                           Goldschmidt, Mr. George B    male   71      0
174                            Smith, Mr. James Clinch    male   56      0
185                              Rood, Mr. Hugh Roscoe    male  NaN      0
209                                   Blank, Mr. Henry    male   40      0
284                         Smith, Mr. Richard William    male  NaN      0
445                          Dodge, Master. Washington    male    4      0
475                        Clifford, Mr. George Quincy    male  NaN      0
556  Duff Gordon, Lady. (Lucille Christiana Sutherl...  female   48      1
583                                Ross, Mr. John Hugo    male   36      0
599       Duff Gordon, Sir. Cosmo Edmund ("Mr Morgan")    male   49      1
630               Barkworth, Mr. Algernon Henry Wilson    male   80      0
647                Simonius-Blumer, Col. Oberst Alfons    male   56      0
806                             Andrews, Mr. Thomas Jr    male   39      0
867               Roebling, Mr. Washington Augustus II    male   31      0

     Parch    Ticket     Fare Cabin Embarked
23       0    113788  35.5000    A6        S
96       0  PC 17754  34.6542    A5        C
174      0     17764  30.6958    A7        C
185      0    113767  50.0000   A32        S
209      0    112277  31.0000   A31        C
284      0    113056  26.0000   A19        S
445      2     33638  81.8583   A34        S
475      0    110465  52.0000   A14        S
556      0     11755  39.6000   A16        C
583      0     13049  40.1250   A10        C
599      0  PC 17485  56.9292   A20        C
630      0     27042  30.0000   A23        S
647      0     13213  35.5000   A26        C
806      0    112050   0.0000   A36        S
867      0  PC 17590  50.4958   A24        S

おわりに

書いたら覚えるだろう、そう思っていた。

参考

  1. for文の終了時の処理(for...else) - 繰り返し - Python入門
  2. argparseを使ってみた - そこはかとなく書くよ。
  3. Python: コマンドラインパーサの Click が便利すぎた - CUBE SUGAR CONTAINER
  4. 現在時刻を取得・ファイル名用に整形 - 週末京都
  5. スーパーIT戦士になりたい!: Python: ファイル名の番号でソートする
  6. Chainerによる多層パーセプトロンの実装 - 人工知能に関する断創録
  7. numpy.random.permutation — NumPy v1.10 Manual
  8. sklearn.cross_validation.KFold — scikit-learn 0.17 documentation
  9. API Reference — scikit-learn 0.17 documentation

kaggle2回目 タイタニック号の生存者予測

概要

kaggle2回目。タイタニック号の乗員乗客の生存者予測。 データを眺めて思いついたことをやってみると、ベースラインは超えられたけれど、今ひとつうまくいかなかった。 今後、続きをやるか、別の課題をやるか、全然違うことをやるかは考え中。

はじめに

kaggle初挑戦ということで、前回手書き文字認識をやった。 課題自体やり尽くされている感じであまりおもしろくないので、今回はタイタニック乗客の生存予測に挑戦することにする。 手法は諸々の事情からニューラルネットワークにした。

どういう課題か

タイタニック号の沈没事故によって、多くの人が亡くなっている。 今回の課題は、この沈没事故によって死亡した人と生き残った人とを分類する問題である。 kaggleから与えられる情報は以下。

  • PassengerId: kaggleが振ったただの連番。役に立たない気がする。
  • PassengerClass: 乗客の等級。1から3まで。金持ちが生き残りそう。
  • Name: 名前。"Mr."とか"Don."とか色々情報がある。うまくやれば家族の情報も抽出できるかも。
  • Sex: 性別。映画では女性、子どもが先に救命ボートに乗ったらしい。関係ありそう。
  • Age: 年齢。同上。
  • Sibsp: 海外に住んでいる兄弟、配偶者の数。なんやこれ
  • Parch: 海外に住んでいる両親、子どもの数。なんやこれ
  • Ticket: チケットの番号。"113803"とか"A/5 21171"とか規則がよくわからない。
  • Fare: 運賃。金持ちが生き残りそう。
  • Cabin: 部屋番号。"C85"みたいな感じ。ぱっと見た感じほとんど欠損している。
  • Embarked: 乗船した港。Cherbourg、Queenstown、Southamptonの3種類。

ぱっとみると、関係ありそうなデータからそうでもないデータまで色々ある。 課題用のものというより、事後処理のときの名簿か何かなのかもしれない。

kaggleから提示されるベースラインは4つ。

  1. Assume All Perished [0.62679] 全員死んだとしたときの精度。
  2. Gender Based Model [0.76555] 女性が全員助かって、男性が全員死んだとしたときの精度。
  3. My First Random Forest [0.77512] Name、Ticket、Cabinのデータ以外を使ってランダムフォレスト。コードも公開。
  4. Gender, Price and Class Based Model [0.77990] 性別、チケットのクラス、チケットの値段の割合に基づいて予測(?)

とりあえず、これら全てのベースラインを超えることを目標とする。

どういうデータか

とりあえず、生データを眺めて以下のことに気づいた。

  • 欠損値が多い。: Ageは結構抜けていて、Cabinにいたっては欠損している方が多い。
  • Fareが0の人がいる。: 欠損しているのか、来賓的なやつなのか。どう扱っていいものか迷う。
  • 同じ名字、部屋番号の人がそれなりにいる。: 家族なら生死を共にしている可能性が高い気がする。

欠損値の補完の方法によって精度が結構変わりそう。 kaggleのコードでは文字列的な情報は使っていないので、欠損値の補完は、単に中央値を用いていた。 補完せずにデータを捨てると、いくつぐらいになるのかを計算してみるとデータ数が891→183になった。 もともと多くはないデータ数をかなり削ることになってしまうので、この案はなし。

とりあえず、色々ヒストグラムを書きながら考えることにする。 まず、データを生存か死亡かで分割する。

import pandas as pd
import pylab as plt

data = pd.read_csv("train.csv").replace("male",0).replace("female",1)
data = data.replace("C",0).replace("Q",1).replace("S",2)
split_data = []
for did_survive in [0,1]:
    split_data.append(data[data.Survived==did_survive])

以下のヒストグラムは全て、青色が死亡した人、緑が生き残った人。 まずはAge。

temp = [i["Age"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=16)

思ったよりも5歳以下の子どもが助かる確率が高い。優先的に逃げれたのだろうか。 一方で、5歳から10歳になると突然低くなる。

60歳超えるとだんだん厳しくなると思いきや、75歳以上の人が1人助かっている。 少ない訓練データにこういうのがあると、すごい高齢の人は助かるみたいな謎の学習しそうで怖い。 なんにせよ年齢はかなり使えそう。

f:id:ksknw:20151129151059p:plain

次にPassenger Class。

temp = [i["Pclass"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=3)

やはり一等客の生存率が高い。お金の力か。 f:id:ksknw:20151129151101p:plain

続いて性別。

temp = [i["Sex"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=3)

男(左)はだいたい死んでいて、女性(右)はかなり助かっている。 救命ボートの数が足りなかったらしいタイタニックにおいて、レディファーストとかやってたんだろうか。 この値もかなり使えそう。 f:id:ksknw:20151129151103p:plain

次に運賃。

temp = [i["Fare"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=10)

これもある程度出していたほうが、生存率高くなりそう。 f:id:ksknw:20151129151105p:plain

ここから、役に立つかよくわからないデータ。 一つ目、国外に住む兄弟または配偶者の数

temp = [i["SibSp"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=8)

グラフを見る限り、0人よりも1人のほうが生存率高そうなんだけど、この値が何と関係しているのかわからない。 国外に家族がいたほうが、収入が高くなって、高い運賃を払いやすくなるんだろうか。 f:id:ksknw:20151129151107p:plain

次、国外に住む親または子どもの数 一方で、こちらはパッとしない。

temp = [i["Parch"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=8)

f:id:ksknw:20151129151109p:plain

乗船した港。 これも意外と役に立ちそうな値。 地域による所得差とか、売られていた部屋の分布が偏っているとかそういうのだろうか。

temp = [i["Embarked"].dropna() for i in split_data]
plt.hist(temp, histtype="barstacked", bins=3)

f:id:ksknw:20151129151112p:plain

My First Neural Network

自分的ベースラインとして、とりあえず、ニューラルネットに突っ込んでみる。 使った特徴量はPassengerClass,Sex,Age,Sibsp,Parch,Fare,Embarked。 欠損値はなんとなく平均値で補完した。 PassengerClassとEmbarkedは1ofKにした。 結果60%ぐらいで、これは先の男女だけで判別したモデルよりも低い。雑魚すぎる。 (実は正規化のところでやらかしていたので普通にミスだった)

欠損値の補完とはなにか

欠損していた場合、モデルにとってニュートラルな信号を入力すべきと思うんだけど、 例えば年齢を欠損していた時、これを平均値もしくは中央値で補完すると、死亡と判定する率が上がりそうな気がする。

調べてみると、欠損値を補完する手法が色々あるらしいが、面倒なので欠損の場所によってモデルを分けることにした。 上記の特徴量に関して、Testデータで欠損している部分を調べると、Ageが欠損しているパターンとFareが欠損しているパターンがあった。 よって3つのモデルを学習させることにする。

  1. Age、Fareが欠損していないデータを学習データにして、同様のデータを予測。
  2. Fareが欠損していないデータを学習データにして、Ageが欠損しているデータを予測。
  3. Ageが欠損していないデータを学習データにして、Fareが欠損しているデータを予測。

3つのネットワークを学習させて結果をsubmitすると、 0.77990でベースラインに並んだ。

f:id:ksknw:20151129150618p:plain

ここまで適当にしていたパラメータ調整をちゃんとやればベースラインを超えられる気がしたので、 パラメータ色々いじってsubmitしまくると、 0.78469になり、無事ベースラインを超えることができた。

f:id:ksknw:20151129150625p:plain

試行回数をふやす

何度かパラメータ調整をしてsubmitを繰り返していると、スコアが結構ブレることに気づく。 また、先ほどのモデルの1.に使う学習データは結構少ないので、精度が本当にでるのかやや不安でもある。 そこで、学習を複数回行い、また、1.のモデルで判別していたデータを2.のモデルでも判別した。 それらの結果の平均を最終的な判別結果とする。

submitすると、0.77990で、うーんという感じ。 結果の分散が減ったぶん、たまたまうまくいくのもなくなっているという印象。

別の特徴量を考える。

このままモデルをいじっていてもジリ貧な気がしたので、別な特徴量を考えることにした。

  • PassengerId: なんとなく削っていたけれども、もしかすると何かあるかもしれないので足す。
  • 欠損値の有無: そもそも欠損しているのはなぜか。死んだからではないのか。
  • Nameの敬称: Mr.とかMs.とかだけだと思っていたら、こちらによるともっと色々あるらしいので、突っ込んでみることにした。

以上を加えてモデルを学習させると 0.76077でいまいちだった。ぐぬぬ

ちなみに現在のコードのだいたい。

#-*- coding: utf-8 -*-
import pandas as pd
import numpy as np
from nn import Perceptron
from chainer import cuda
import sklearn.preprocessing


def to_feature(data):

    emberked = []
    for i in data["Embarked"]:
        temp = np.zeros(3)
        if not np.isnan(i):
            temp.put(i, 1)
        emberked.append(temp)
    feature = emberked

    pclass = []
    for i in data["Pclass"]:
        temp = np.zeros(3)
        if not np.isnan(i):
            temp.put(i - 1, 1)
        pclass.append(temp)
    feature = np.c_[feature, pclass]

    feature = np.c_[feature, data["Age_na"]]
    feature = np.c_[feature, data["Cabin_na"]]

    for keyword in ["Master.", "Col.", "Mrs.",
                    "Ms.", "Miss.", "Rev.",
                    "Mr.",  "Dr.", "Major."]:

        feature = np.c_[feature, data.Name.str.contains(keyword).astype(int)]

    try:
        feature = np.c_[feature, data["Age"]]
    except:
        pass
    try:
        feature = np.c_[feature, data["Fare"]]
    except:
        pass
    try:
        feature = np.c_[feature, data["SibSp"]]
    except:
        pass
    try:
        feature = np.c_[feature, data["Parch"]]
    except:
        pass
    feature = np.c_[feature, data["Sex"]]
    feature = np.c_[feature, data["PassengerId"]]

    return feature


def read_data(data_filename):
    data = pd.read_csv(data_filename).replace("male", 0).replace("female", 1)
    data = data.replace("C", 0).replace("Q", 1).replace("S", 2)
    return data


def add_is_na_columns(data):
    for label in ["Age",  "Cabin"]:
        temp = pd.DataFrame(pd.isnull(data[label]).astype(int))
        temp.columns = [label + "_na"]
        data = pd.concat([data, temp], axis=1)
    return data

if __name__ == '__main__':
    test_data = add_is_na_columns(read_data("test.csv"))
    train_data = add_is_na_columns(read_data("train.csv"))

    label_categories = [['PassengerId', 'Name', 'Pclass', 'Sex', 'Age', 'SibSp', 'Parch', 'Fare',  'Embarked', 'Cabin_na', 'Age_na'],
                        ['PassengerId', 'Name', 'Pclass', 'Sex', 'SibSp', 'Parch', 'Fare',  'Embarked', 'Cabin_na', 'Age_na']]

    results = []

    for i in range(1):
        print i
        for label in label_categories:
            temp_train_data = train_data.ix[:, train_data.columns.isin(label + ["Survived"])].dropna()
            temp_test_data = test_data.ix[:, test_data.columns.isin(label)].dropna()

            train_vec = to_feature(temp_train_data)
            test_vec = to_feature(temp_test_data)

            min_max_scaler = sklearn.preprocessing.MinMaxScaler()

            no_normalize_index = 8 + 9
            min_max_scaler.fit(train_vec[:, no_normalize_index:])

            test_vec[:, no_normalize_index:] = min_max_scaler.transform(test_vec[:, no_normalize_index:])
            train_vec[:, no_normalize_index:] = min_max_scaler.transform(train_vec[:, no_normalize_index:])

            nn = Perceptron(len(train_vec[0]), 2, 700, 45, 45)
            results.append([temp_test_data.PassengerId,
                            nn.train(train_vec,
                                     temp_train_data["Survived"],
                                     test_vec)])
            # from sklearn.ensemble import RandomForestClassifier
            # random_forest = RandomForestClassifier()
            # random_forest.fit(train_vec, temp_train_data["Survived"])
            # results.append([temp_test_data.PassengerId,
            #                 random_forest.predict(test_vec)])

    results.append([[1044], [1, 0]])  # Fareが欠損している一人は死んだことにする
    # results.append([[1044],  0])  # Fareが欠損している一人は死んだことにする

    print results
#    f = open("result_random_forest.csv", "w")
    f = open("result_nn_10.csv", "w")
    f.write("PassengerId,Survived\n")

    for result in results:
        for i in zip(*result):
            print i
            f.write("%d,%d\n" % (i[0], np.argmax(i[1])))
#            f.write("%d,%d\n" % (i[0], i[1]))
    f.close()

終わりに

なんか適当にパラメータいじっていたら0.78947になった。 たまたま感ある。 f:id:ksknw:20151129150629p:plain

色々やってみたけれども、なかなか難しい。 特徴量を増やしてみても、よくなるわけでもなく。 というかそもそも、どれぐらいいけるものなのかもわからない。 ランキングの一番上は1.000なので、萎える。

もうちょっとアイデイアがあるので、もしかすると続くかもしれないし、 飽きてきている感があるので、別の課題をやるかもしれない。

参考

kaggleで予測モデルを構築してみた (1) - kaggleって何? - About connecting the dots.