TIM Labs
前回、28x28の手書き数字のAutoEncoderプログラムの説明をした。
そして走らせて、やっと(笑)画像が集まったので、ご覧いただこう。

10000エポックでは、この程度の画像にしか再現されない。
mnistae10000SGD.png1万エポック程度では、何が何だか分からない。
それで、一気に10万エポックまで飛ぼう。

mnistae100000SGD.png10万エポックになると、半数くらいは判別できそうだが、まだとっても苦しいところだ。

30万エポックになると、かなり数字がはっきり見えてくるようになる。
mnistae300000SGD.png30万エポックになると、かなり正確に元のイメージを再現するようになってくる。

そして、100万エポックになると、ここまで再現できるようになる。

mnistae1000000.png100万エポックになると、文句なく十分な再生が行われたと判断できるだろう。

しかし、まずい。
エポック数が100万とは、100万回も学習を繰り返さないと使えないということだ。つまり、超バカだ。

もっとエポック数が少なくて、ちゃんと再現される方法はないものだろうか。


プログラムの初期化部分に optimizerの選択があるようだ。
# Initialize model        
model = MyAE()
optimizer = optimizers.SGD()
optimizer.setup(model)
SGDというオプティマイザーが選択されているようだ。
SGDとは stochastic gradient descent(確率的勾配降下法) のことだが、面倒なので説明は省略。まあ、最初に説明されることが多い最適化方法である。

もしかして、ここを書き換えると、もっと上手くいくのだろうか?
MNISTの28x28の手書き数字は、現実の処理でも使うくらいのドット数がある。
これでAutoEndoderを作って、どのくらい再現性があるか調べてみよう。

まず、クラスMyAeのノード数の部分を書き換えよう。
中間層は、40ノードとする。つまり、748 --> 40 --> 748 と変換していこう。

class MyAE(Chain):
    def __init__(self):
        super(MyAE, self).__init__(
            l1=L.Linear(784,40),
            l2=L.Linear(40,784),
        )
変更はたったこれだけで、28x28の画像に対応できるようになる。
後は、書くだけだ。

最初に、データを読み込む。

# http://yann.lecun.com/exdb/mnist/
train, test = chainer.datasets.get_mnist()
xtrain = train._datasets[0]
ytrain = train._datasets[1]
xtest = test._datasets[0]
ytest = test._datasets[1]
今回は、学習の途中で、一定エポック毎にAutoEncoderの結果の最初の48枚の画像を1つの画像ファイルまとめて出力した。

# Learn
losslist = []
for j in range(1000000):
    x = Variable(xtrain[:10000])
    model.cleargrads()             # model.zerograds() 非推奨
    loss = model(x)
    if j%10000 == 9999:
        print( "%6d   %10.6f" % (j+1, loss.data) )
        xx = Variable(xtrain[:48], volatile='on')
        yy = model.fwd(xx)
        plotresults( yy, "mnistaeout/mnistae%d.png" % (j+1) )

    losslist.append(loss.data)     # 誤差をリストに追加
    loss.backward()
    optimizer.update()
10000エポック毎にスナップショット画像を吐き出して、全体で100万エポックまでやったのだが、このくらいやるとしっかり時間がかかり、走らせて結果は翌日確認することになった。(まだGPUは使っていない)

といことで、今回はプログラムの紹介だけで、結果は次回に示す。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第5章、GLMの尤度比検定と検定の非対称性についてのまとめの三回目です。

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 なので、パラメトリックブートストラップ法でどのようにデータがシミュレートされているか見るためのコードを用意しました。 コードはRで書きました。

linkFunction <- function(b1, b2) {
  function(x) {
    exp(b1 + b2 * x)
  }
}

createData <- function(n, b1, b2) {
  x <- runif(n, 1, 20)
  y <- rpois(n, lambda = linkFunction(b1, b2)(x))
  data.frame(x, y)
}

plotGraph <- function(data, fit1, fit2, title) {
  plot(data$x, data$y, type = "p", main = title, xlab = "x", ylab = "y")
  lines(0:20, linkFunction(fit1$coefficient[1], 0)(0:20), lty = "solid")
  lines(0:20, linkFunction(fit2$coefficient[1], fit2$coefficient[2])(0:20), lty = "dashed")
  nl <- paste0("y = exp(", round(fit1$coefficient[1], 3), ")")
  xl <- paste0("y = exp(", round(fit2$coefficient[1], 3), " + ", round(fit2$coefficient[2], 3), " * x)")
  legend("topleft", legend = c(nl, xl), lty = c("solid", "dashed"))
}

plotBootstrap <- function(b1, b2) {
  simulatedData <- createData(100, b1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = simulatedData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = simulatedData)
  title <- paste0("Bootstrap (b1 = ", b1, ", b2 = ", b2, ")")
  plotGraph(simulatedData, fit1, fit2, title)
}

expData <- createData(20, 0.1, 0.07)
fit1 <- glm(formula = y ~ 1, family = poisson, data = expData)
fit2 <- glm(formula = y ~ x, family = poisson, data = expData)
plotGraph(expData, fit1, fit2, "Experimental data")

plotBootstrap(fit1$coefficient[1], 0)
plotBootstrap(fit2$coefficient[1], fit2$coefficient[2])

コードを実行するとグラフが3枚プロットされます。 ぜひ実行してみてください。 実行結果の解説は次回にまわします。

CUDAをインストールしたら、GPUが動作しているかどうか、とりあえず確認しよう。
nvidia-smi というコマンドが入ったはずで、それを実行してみよう。

$ nvidia-smi
Sat Mar 18 18:46:16 2017       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 375.26                 Driver Version: 375.26                    |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|===============================+======================+======================|
|   0  GeForce GTX 105...  Off  | 0000:01:00.0     Off |                  N/A |
| 35%   37C    P0    35W /  75W |    337MiB /  4038MiB |      1%      Default |
+-------------------------------+----------------------+----------------------+
                                                                               
+-----------------------------------------------------------------------------+
| Processes:                                                       GPU Memory |
|  GPU       PID  Type  Process name                               Usage      |
|=============================================================================|
|    0      1129    G   /usr/lib/xorg/Xorg                             164MiB |
|    0      1919    G   compiz                                         127MiB |
|    0     11883    C   gimp-2.8                                        43MiB |
|    0     12492    G   unity-control-center                             1MiB |
+-----------------------------------------------------------------------------+

上図のように、2段に別れて表示される。
現在はGPUが1個だけなので、GPU0として表示されている。たくさん同時に使うと、0,1,2,...となるはず。
FANや、温度、GPUのメモリ容量など色々表示されるので、GPUは機能しているようだ。

下段の表を見ると、すでにいくつかのソフトが勝手にGPUを使っていることが分かる。
明示的にGPUを使わなくても、様々なアプリがGPUを使うことがあるようだ。

とりあえず、動いているらしいことが分かったら、Cのプログラムをコンパイルして動作確認すべきだ。
といっても、Cのプログラムをいきなり書くのは難しい。

こういうときの場合に、CUDAのサンプルプログラムが用意されている。

普通にインストールすると、/usr/local/cuda/samples というディレクトリの中に分野別に分けてプログラムが入っているみたいだ。

fuji@fujigpu:/usr/local/cuda/samples$ ls
0_Simple     2_Graphics  4_Finance      6_Advanced       EULA.txt  bin     err
1_Utilities  3_Imaging   5_Simulations  7_CUDALibraries  Makefile  common

その中に、Makefileというのがあるので、さっそくmakeしてみよう。 すると、binの下に実行可能ファイルがいっぱいできる。
どんなプログラムかは、CUDA code sampleで概要が分かるだろう。
とにかく、適当にクリックして実行してみよう。

これは、3Dのグラフィックスで、実際には刻々と形が変化し続ける。
 
CudaMarchingCubes.png
他にもいろいろあるのだが、説明するのは面倒だし、コンパイルし実行できることが一応分かったので、CUDAのインストールの確認はこれで終わることにする。


このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第5章、GLMの尤度比検定と検定の非対称性についてのまとめの二回目です。

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 前回はパラメトリックブートストラップ法がどのように動くのかを見るためのコードを紹介しました。 コードを実行すると結果は次のようになります。 何度か実行してみて、だいたいの傾向をつかんでみてください。

data size 5 , lambda = exp(0.01 + 0.03x) 
P = 0.956 
Failling to reject the null hypothesis.

data size 5 , lambda = exp(0.01 + 0.1x) 
P = 0.098 
Failling to reject the null hypothesis.

data size 20 , lambda = exp(0.01 + 0.03x) 
P = 0.128 
Failling to reject the null hypothesis.

data size 20 , lambda = exp(0.01 + 0.1x) 
P = 0 
The null hypothesis is rejected.

だいたいの傾向としてサンプル数が多くて x の相関が強いほど帰無仮説は棄却されやすくなります。 データを生成する真の分布は x と相関があるので、帰無仮説は棄却されて欲しいですが、そうならない場合があるわけです。 その様子をグラフにもしてあるので、紹介します。 まずは帰無仮説が棄却されて、無事正しいモデルが支持された場合です。

bootstrap_sample2.png

検定に使われたデータが丸、実線が帰無仮説、点線が対立仮設です。 対立仮設の x の係数は 0.098 と推定されています。 おかげで対立仮設は帰無仮説と比べて x と相関のあるデータによく当てはまり、二つのモデルの逸脱度の差が大きくなってP値が低くなり尤度比検定で正しいモデルが選択されたわけです。

一方、帰無仮説が棄却できなかった場合はこのような図になります。

bootstrap_sample1.png

一見して帰無仮説も対立仮設も違いが見えません。 対立仮設の x の係数は 0.008 と推定されていて 0 と違いはほとんどありません。 とすると、帰無仮説でも対立仮設でもデータへの当てはまりの良さに違いはほとんどなく、十分なP値が得られなかったということです。

だいたいの感覚として、帰無仮説が棄却されたりされなかったりする時はこういうことが起きています。

MNISTの手書き数字の読み込みは、きわめて簡単である。
最初におまじないを並べたあと、次の1行だけで、トレーニングデータとテストデータが読み込まれ、2つのオブジェクトに入る。
train, test = chainer.datasets.get_mnist()

トレーニングデータがどのように入っているか、確認しよう。
>>> type(train)
<class 'chainer.datasets.tuple_dataset.tupledataset'>
>>> len(train)
60000
>>> train[0]
(array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,

.........中略..........
0. , 0.11764707, 0.14117648, 0.36862746, 0.60392159, 0.66666669, 0.99215692, 0.99215692, 0.99215692, 0.99215692, 0.99215692, 0.88235301, 0.67450982, 0.99215692, 0.94901967, 0.76470596, 0.25098041, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.19215688, 0.9333334 , 0.99215692, 0.99215692, 0.99215692, 0.99215692, 0.99215692, .........中略.......... 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], dtype=float32), 5) >>> train[0][0].shape (784,)
これから、trainは要素数60000個のリストで、各リストは、画像データと、数値(0から9)のタプル。
画像データは、要素数784個の1次元配列。

今回は、最初の48個の画像だけを表示するので、最初の48個だけを取り出す。
xtrain = train._datasets[0][:48]
ytrain = train._datasets[1][:48]
あとは、subplotsを使って、画面を分割して表示するだけである。
fig,ax = plt.subplots(nrows=6,ncols=8,sharex=True,sharey=True)
ax = ax.flatten()
for i in range(48):
    img = xtrain[i].reshape(28,28)
    ax[i].imshow(img,cmap='Greys',interpolation='none')
ax.flatten()により、forループ中で、表示枠を指定するのに、単にax[i]で済ますことができる。
画像データは、サイズが784の1次元配列になっているので、28x28の2次元配列に直している。
画像データは賢く表示してくれると困るので、ありのまま表示するようにinterpolation='none'を指定している。
後は、ファイルにセーブし表示している。
ax[0].set_xticks([])
ax[0].set_yticks([])
plt.tight_layout()
plt.savefig("mnistdisp48.png")
print(ytrain.reshape(6,8))
plt.show()
plt.subplots()で、分割された枠のx軸y軸の情報を共有するために、sharex、shareyをTrueにしている。そして、最後のところで、set_xticks([])、set_yticks([])により、目盛りなど目障りなものを表示しないようにしている。

ということで、最後にプログラム全体を示す。
NVIDIAのGPU、GEFORCEの装着については先日書いた。
次にすべきことは、ちゃんとGPUをプログラムでいじれるか確認することである。

GPUのソフト開発ツールを入れてあれこれ試せばよいのだが、その前にOSを入れないといけない。
AI、 Deep Learning などに使おうとしている、そしてPython, 国産Deep LearningフレームワークのChainer を使おうとしているので、こういう組み合わせの場合よく選択されるLinux/Ubuntuを入れることにした。

Download Ubuntu Desktopに行って、Ubuntu 16.04 LTSを他のマシンで落としてきてDVDに焼いて、それを使ってインストールした。
このあたりのインストール手順についてはネット上に溢れているので、省略する。

しかし、今回、Ubuntuのインストールに何故か失敗した。
もう何回インストールしたか記憶にないのだが、そのため、インストールしながら、他の雑用なども並行しておこなっていた。
そして、インストールが終わって、リブートしたら、動きがおかしい。画面がもさっとした動きしかできない。ディスプレイはGPUボードではなく、マザーボードからやっていて、マザーのディスプレイが超ダメだからこんなに遅いのかとも思ったが、インストール中の反応は悪くなかった。
UbuntuSystemSettings.png
そして、どうしようもないのが、システム設定のアイコンを押しても全く反応がないのである。
これでは、なんの設定変更も困難である。

どうも、ながらインストールしたのがいけないらしい。手抜きをし過ぎて、Ubuntuが嫌がらせをしたようであるので、心を入れ替えて、再インストールすることにした。


すると、今度は何の問題もなく、ちゃんとインストールできた。
ここまでは単なる準備である。

次に、NVIDIAのGPU、GeForceのソフトウエア開発環境をインストールをしないといけない。
AIのソフトはPythonで書くつもりなのだが、PythonからGPUを使う場合には、Cで書かれたプログラムが呼ばれる。まあ、PythonなどでGPUを動かしていたら、せっかくのGPUの速度が生かせないので、GPUを使うにはCを使うのだ。

ということで、NVIDIAが用意しているNVIDIA(自社)のGPUの開発環境であるCUDAをインストールし、GPUがちゃんと動いているか確認することにした。

まず、CUDAをダウンロードするなり、オンライン・インストールしないといけない。
なにはともあれ、CUDAのダウンロードページに行くしかないのだ。

今回は、実行環境は次のように選択した。
CudaSelectPlatform.png
全部選択し終えると、次のように3つのインストール選択肢が出てくる。 そして、サイズは1.9GBだと。巨大だ!
Installation Instructions:

 1. `sudo dpkg -i cuda-repo-ubuntu1604-8-0-local-ga2_8.0.61-1_amd64.deb`
 2. `sudo apt-get update`
 3. `sudo apt-get install cuda`
どれを選んでも良いと思うが、今回は1.で、一気に1.9GBダウンロードしてからインストールしてみた。
このインストールについても、ネット上に説明が転がっているし、失敗するような面倒なことはない。ということで、インストールの詳細は省略する。

さて、インストールは済んだはずなので、動作確認をしよう。(次回につづく)

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第5章、GLMの尤度比検定と検定の非対称性についてのまとめの一回目です。

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 なので、パラメトリックブートストラップ法がどのように動くのかを見るためのコードを用意しました。 コードはRで書きました。

linkFunction <- function(b1, b2) {
  function(x) {
    exp(b1 + b2 * x)
  }
}

createData <- function(n, b1, b2) {
  x <- runif(n, 1, 20)
  y <- rpois(n, lambda = linkFunction(b1, b2)(x))
  data.frame(x, y)
}

plotGraph <- function(data, fit1, fit2) {
  plot(data$x, data$y, type = "p", main = "Bootstrap", xlab = "x", ylab = "y")
  lines(0:20, linkFunction(fit1$coefficient[1], 0)(0:20), lty = "solid")
  lines(0:20, linkFunction(fit2$coefficient[1], fit2$coefficient[2])(0:20), lty = "dashed")
  nl <- paste0("y = exp(", round(fit1$coefficient[1], 3), ")")
  xl <- paste0("y = exp(", round(fit2$coefficient[1], 3), " + ", round(fit2$coefficient[2], 3), " * x)")
  legend("topleft", legend = c(nl, xl), lty = c("solid", "dashed"))
}

bootstrap <- function(b1) {
  simulatedData <- createData(100, b1, 0)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = simulatedData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = simulatedData)
  fit1$deviance - fit2$deviance
}

showBootstrap <- function(dataSize, b1, b2) {
  d <- createData(dataSize, b1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = d)
  fit2 <- glm(formula = y ~ x, family = poisson, data = d)
  dev <- fit1$deviance - fit2$deviance

  devs <- replicate(1000, bootstrap(fit1$coefficients[1]))
  cat("data size", dataSize, paste0(", lambda = exp(", b1, " + ", b2, "x)"), "\n")
  summary(d$y)
  P <- sum(devs >= dev) / 1000
  cat("P =", P, "\n")
  if (P < 0.05) {
    cat("The null hypothesis is rejected.\n\n")
  } else {
    cat("Failling to reject the null hypothesis.\n\n")
  }
  plotGraph(d, fit1, fit2)
}

showBootstrap(5, 0.01, 0.03)
showBootstrap(5, 0.01, 0.1)
showBootstrap(20, 0.01, 0.03)
showBootstrap(20, 0.01, 0.1)

コードを実行すると、ポアソン分布からサンプルデータをサンプリングして、二つのモデルでポアソン回帰をし、尤度比検定を行って帰無仮説が棄却できるかを見ます。 サンプルデータを用意して、ポアソン回帰をしているのは以下のコードです。 createData は平均 exp(b1 + b2 * x) を持つポアソン分布からデータをサンプリングします。 帰無仮説は formulay ~ 1 を使っている方で、対立仮説は y ~ x を使っている方です。 最後に逸脱度の差を計算しています。

  d <- createData(dataSize, b1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = d)
  fit2 <- glm(formula = y ~ x, family = poisson, data = d)
  dev <- fit1$deviance - fit2$deviance

パラメトリックブートストラップ法では、推定した帰無仮説のパラメータを使ってデータをシミュレートしますが、それをやっているのは以下のコードです。 b1 には推定したパラメータが渡されてきます。 そして平均 exp(b1) を持つポアソン分布からデータをサンプリングして、そのデータについてポアソン回帰をやり直し、逸脱度の差を計算します。

bootstrap <- function(b1) {
  simulatedData <- createData(100, b1, 0)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = simulatedData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = simulatedData)
  fit1$deviance - fit2$deviance
}

あとはこのシミュレーションを繰り返し、最初の逸脱度の差がどれだけ珍しいか確率を計算して帰無仮説を棄却できるかどうか結論します。 コードを何回か実行してみると、だいたいの傾向としてサンプル数が多くてx の相関が強いほど帰無仮説は棄却されやすくなるのが分かると思います。 実行結果の解説は次回やります。

sklearnの8x8ドットの手書き文字を使うのをChainer:iris以外のデータでDeep Learning で紹介したが、ここで最初の48文字の画像を示す。
digitsdisp48.png これは、以下の数字を示しているのだが、相当無理があろう。
[[0 1 2 3 4 5 6 7]
 [8 9 0 1 2 3 4 5]
 [6 7 8 9 0 1 2 3]
 [4 5 6 7 8 9 0 9]
 [5 5 6 5 0 9 8 9]
 [8 4 1 7 7 3 5 1]]

ということで、もっと良い数字画像データセットを使うことにする。

よく使われるものに、THE MNIST DATABASE of handwritten digits がある。
The MNIST database (Mixed National Institute of Standards and Technology database) は、28x28の手書き数字のデータ・セットで、トレーニング用が6万文字、テスト用が1万文字ある大規模なものである。
全部使うとシステムが重くなったりするので、一部だけを利用することも多い。

とりあえず、トレーニングセットの最初の48文字を上と同じ形式で示す。
mnistdisp48.png
[[5 0 4 1 9 2 1 3]
 [1 4 3 5 3 6 1 7]
 [2 8 6 9 4 0 9 1]
 [1 2 4 3 2 7 3 8]
 [6 9 0 5 6 0 7 6]
 [1 8 7 9 3 9 8 5]]
全然画像の細かさが違うのが分かるだろう。

さて、これをどうやって読み込み、並べて表示したプログラムについては、次回に説明しよう。
AIをいじっていると、どうしても計算時間で苦しめられるようになる。ちょっと勉強するにしても、何とか計算時間を短縮できないものかと思うようになる。
AIでは、大量の計算をCPUではなく、GPU(グラフィックプロセッサ)で行うことが普通だ。
また、大量のGPUが使えるクラウドサービスも色々ある。
東京工業大学のTSUBAMEのように、スパコンを細かく分けて使えるようにしているサービスもある。詳しくは、TSUBAME計算サービスを参照のこと。運用状況なども公開されている。

しかし、もうちょっと身近に、かつ勝手に使ってみないと、細かいことを実感できないなと思っていると、こんなものが先日机の上に置かれていた。

GPUといえば、今は圧倒的にNVIDIAが有名だ。多くのスパコンでも使われていて、Deep Learning は、NVIDIAのGPUの上で動いていると言ってもだいたい当たっている。
....などと考えていると、

GTX0150Ti-1.jpg
マザーボードと、NVIDIAのGPUである GEFORCCE GTX 1050Ti と、真っ赤なメモリが机の上に置かれていた。

CPUは、そのへんのをつければ動く。GPUでも超省エネの場合は普通のミニタワー型パソコンでも大丈夫だが、これは2スロット分の幅があり、電気もそれなりに食うので、電源がしっかりしているボックスでないとダメだ。
とりあえず、750ワットのケースの中身を入れ替えて、そこいらのi5を挿してみた。

GTX1050Ti-3.jpgGTX1050Ti-2.jpg
GTX1050Ti-4.jpgGTX1050Ti-5.jpg

電源を入れると、暗いケースの中で、GPUボードのLEDが怪しく光るようになった。ちゃんと動いているということだろうか。

ちゃんと動作しているかは、起動時にDELキーを押してBIOS設定画面に入ると分かりやすい。
どんどんグラフィカルな制御盤みたいになっていくな。

GTX1050-Ti-6.jpgということで、インストールについては、次回(GPU編)、書こう。

最近のコメント