TIM Labs

2017年3月アーカイブ

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第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"))
}

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
}

for(i in 1:10) {
  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)
  dev <- fit1$deviance - fit2$deviance

  P <- mean(replicate(1000, bootstrap(fit1$coefficient[1])) > dev)
  title <- paste("P value", P)
  plotGraph(expData, fit1, fit2, title)
}

コードを実行するとグラフが10枚プロットされます。 それぞれにポアソン回帰で使われるデータと、推定されたモデルをプロットしてあります。 あと、パラメトリックブートストラップ法でP値を計算してタイトルに入れてあります。 典型的には下図のようなグラフがプロットされます。

P_value_0.png

ですが、時々は下図のようなグラフがプロットされます。

P_value_0191.png

有意水準として5%を設定するなら、P値が0.05以上だと帰無仮説を棄却できません。 棄却できる時は対立仮説は帰無仮説とはだいぶ形の違う曲線になりますが、棄却できない時は直線に近い形になるのが分かります。 データは x に相関のある分布から生成されていますが、たまたまデータが偏っていると最尤推定で x と相関のないようなパラメータを推定してしまう場合があるわけですね。

28x28のオートエンコーダだが、ちゃんとエンコードできるまでには100万エポック程度教えこまないといけない。
つまり、とても学習がのろく、馬鹿である。
もうちょっと賢く、さっさと学習できるようにしたいものである。
# Initialize model        
model = MyAE()
optimizer = optimizers.SGD()
optimizer.setup(model)
ここで、Chainerのマニュアルの該当部分を読もう。
Docs ≫ chainer Reference Manual ≫ Optimizers

どうやら何種類も用意されているようで、何を選ぶべきか困る。 いままでは、optimiers.SGD()により、SGD(Vanilla Stochastic Gradient Descent)というタイプのOptimizerを生成していたようだ。 でも、これは頭が悪そうなので、別のにしよう。 Deep Learningの本で時々出てくるAdamというのがあるので、とりあえずこれを使ってみよう。
変更は、以下の1行のみである。

optimizer = optimizers.Adam()

Adamの説明が、"Adam optimization algorithm." となっているが、これではどんな最適化が行われるのかさっぱり分からない。 Adamの説明についてもネット上に情報があるので、ここでは説明を省略し、さっそく実行してみよう。

とりあえず、前回と同じ10000エポック時の出力を示す。
mnistae10000Adam.png
SGDでは10000エポックでは何も分からない状態だったが、今回のAdamでは、もう十分な再現性があり、数字が判別できる。SGDの100万エポックに近い状態である。 ならば、1000エポックでどこまで判別できるか調べてみた。
mnistae1000Adam.pngたった1000エポックでも、そんなに悪くない。 ということで、Optimizerの選択は非常に重要なようだ。

Optimizerは、AutoEncoderだけでなく、もっと他の場合でも有効であるに違いない。
とうことで、普通のDeep Learning の学習に戻ろう。
トークイベント

岩波データサイエンス Vol.5 刊行記念
特集 「スパースモデリングと多変量データ解析」

ちょっぴり協力、寄稿しているので、こちらで紹介します。
詳しくは、最終ページをご覧ください。

4月13日(木)19:00から、東京駅八重洲口の南にあるグラントウキョウサウスタワーで第5巻執筆者などによるトークショーです。

参加費無料です。

事前に「岩波データサイエンス Vol.5」を今すぐ入手し、事前に読んでおくと、トークショーの内容の理解がさらに深まるかと思います。

イベント詳細:https://connpass.com/event/53926/

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

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

コードを実行するとグラフが3枚プロットされます。 1枚目は以下のような図です。 createData を使ってサンプリングしたデータと、帰無仮説と対立仮説に対してポアソン回帰して得られた曲線がプロットされています。 このデータは何かの実験データのような、検定にかけたいデータだとします。 帰無仮説は formulay ~ 1 を使っている方で、対立仮説は y ~ x を使っている方です。 データは x と正の相関のあるデータなので、対立仮説は右肩上がりの曲線になっています。

ExperimentalData.png

2枚目は以下のような図です。 パラメトリックブートストラップ法で使う、ポアソン回帰で得られた帰無仮説のパラメータを使ってシミュレートしたデータと、そのデータに対して新たにポアソン回帰して得られた帰無仮説と対立仮説の曲線がプロットされています。 帰無仮説は x と相関のないモデルなので、これを真の分布だと仮定してシミュレートすると、得られるデータは x に相関のないデータになります。 すると、帰無仮説も対立仮説も、ほとんど同じような直線になってしまいます。

Bootstrap_null.png

3枚目は以下のような図です。 パラメトリックブートストラップ法ではこういうシミュレーションは行いませんが、参考までに帰無仮説ではなく、対立仮説の方のパラメータを使ってデータをシミュレートして、2枚目と同じようなプロットを作ってみました。 対立仮説は x と相関のあるモデルなので、そのパラメータを使うと x に相関のあるデータがシミュレートされます。 すると、ポアソン回帰で得られる曲線は、実験データの時と同じような曲線になります。

Bootstrap_x.png

このように、検定にかけたい元の実験データと、パラメトリックブートストラップで行うシミュレーションで生成されたデータは性質が違うデータです。 この違いを利用して尤度比検定は帰無仮説を棄却できるのですね。

前回、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編)、書こう。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 前回はAICを使ってモデル選択の様子を見るためのコードの内容を解説しました。 今回は実行結果の内容を解説します。

コードを実行すると結果は次のようになります。

indipendent data
                    null        x           c           xc
log likelihood       -208.4236   -207.8885   -207.5525   -207.5126 
mean log likelihood  -209.1079   -209.7084   -209.7     -210.2147 
bias                 0.6843544   1.819908    2.14749     2.702091 
aic                  418.8471    419.7771    419.1049    421.0252 

correlated X data
                    null        x           c           xc
log likelihood       -244.4405   -191.2695   -243.681    -190.9884 
mean log likelihood  -246.9158   -193.4544   -248.0273   -193.9776 
bias                 2.475276    2.184949    4.346326    2.989236 
aic                 490.881      386.5389    491.3619    387.9768 

correlated C data
                    null        x           c           xc
log likelihood       -245.6007   -244.4792   -191.0656   -190.7746 
mean log likelihood  -246.9547   -247.9338   -193.4138   -193.9289 
bias                 1.354043    3.454623    2.348202    3.154229 
aic                 493.2014     492.9584    386.1313    387.5493 

correlated XC data
                    null        x           c           xc
log likelihood       -584.6693   -410.6131   -410.5213   -244.1359 
mean log likelihood  -591.7481   -421.7135   -421.5719   -246.7352 
bias                 7.078751    11.10039    11.05061    2.599335 
aic                  1171.339    825.2262    825.0427    494.2717 

結果の見方ですが、上から順に xc の両方に関係のないデータ、x と関係したデータ、c と関係したデータ、xc の両方と関係したデータについて計算した場合です。 それぞれの場合において、4つのモデルに対する最大対数尤度、平均対数尤度、バイアス、AICが計算されています。 モデルは左から順にリンク関数に exp(b1)exp(b1 + b2 * x)exp(b1 + b3 * c)exp(b1 + b2 * x + b3 * c) を使った場合です。

どの場合においても、対数尤度が最大になるのは、xc の両方を説明変数に含むモデルです。 ここには複雑なモデルほどデータによく当てはまるという一般的な性質が表れています。 データが xc の両方に関係する場合には望ましい結果ですが、実際にはデータが xc に関係のない場合にも、関係あるとしてモデルを作った方が特定のデータへの当てはまりは良くなるわけです。 しかし、平均対数尤度は大きくなっていません。 これは、無用に複雑なモデルは推定のために与えられたデータだけに過度に当てはまって、未知のデータへの当てはまりが悪くなっていることを意味します。 尤度だけを使ってモデルの良さを測ることはできないわけです。

そこで、AICの出番です。 AICはモデルのパラメータの個数を当てはまりの良さから割り引いて計算するので、無用なパラメータを増やしても良い数値になりません。 実行結果でも、AICが最小になるモデルを選べば正しくモデルを選択できることが分かります。

AutoEncoder1回分の処理を行う関数は前回説明した。

今回は、そのexeconce()を呼び出して脳内情報を得て、matplotlibを利用して図示しよう。

9回呼び出して、3x3の形に並べてみよう。
まず、図のサイズを決める。デフォルトでは小さいので、サイズ指定した。

plt.figure(figsize=(16,12))
figure()により、図に関するさまざまな事を指定できる。
figsizeは、全体のサイズをインチ単位で指定する。
dpiもあり、デフォルトはdpi=100になっている。
なので、上記は、100dpiで、16インチx12インチになる。
つまり、ドットで言えば、1600ドットx1200ドットとなる。

準備ができたら、9回ループして、実行し、結果を毎回画面の指定区分の中に表示するだけ。

for idx in range(9):
    print("exec ",idx)
    ans = execonce()

    ansx1 = ans[0:50,0]
    ansy1 = ans[0:50,1]
    ansx2 = ans[50:100,0]
    ansy2 = ans[50:100,1]
    ansx3 = ans[100:150,0]
    ansy3 = ans[100:150,1]

    plt.subplot(3,3,idx+1)
    plt.scatter(ansx1,ansy1,c="r",alpha=0.5,linewidth=0)
    plt.scatter(ansx2,ansy2,c="g",alpha=0.5,linewidth=0)
    plt.scatter(ansx3,ansy3,c="b",alpha=0.5,linewidth=0)
ansx1,ansy1,...の部分は、3種類のIrisに合わせて、3つのx座標、y座標に分けているだけ。

plt.subplot()で、画面分割している。
行分割数、列分割数の順に指定し、最初から何番目枠に表示するかを指定する。枠番号は、1から始まる。

plt.scatter()で散布図の形で示した。
最初の2つが、x座標値、y座標値んびなる。
c=色、alpha=透過率、linewith=0 で表示している。
マークを使ったり、色を細かく指定したり、いろいろできるが説明省略。

これだけで、Deep Learning の脳内がちょっと表示できたと言えるかな。

Irisデータで実験したAutoEncoderの、中間層(ノード数2)の中身は次のようになっていた。
各グラフが1回の実験を示し、9回実験したものを並べてみた。
もちろん同じデータで同じ条件で実験したのだが、こんなに違う。
目盛りも図毎にかなり異なる。
つまり、

同じデータを与え、ほとんど似た結果が出ても、中間層の状態は同じにはならない。

中間層は2ノードの層が1つあるだけ。この2つのノードのそれぞれをx軸、y軸とし、Irisの種類にによって、RGBの3色で中間層のノードの値をプロットしてみた。
Irisの種類によって、偏りがあるのがわかる。
とくに、R(赤)は、G、Bとははっきり分かれている。
irisautoencode.pngでは、この図をどうやって作ったか、プログラムの説明をしよう。

今までは、学習し、評価し、それで終わりのプログラムであったが、今回は全データを学習に用い、中間層の中身を取り出して、その中間層の中身を返す関数にオートエンコーダの1回分の処理をまとめ、execonce() とした。

def execonce():
    # Initialize model        
    model = MyAE()
    optimizer = optimizers.SGD()
    optimizer.setup(model)

    # Learn
    n = 150
    for j in range(3000):
        x = Variable(xtrain)
        model.cleargrads()    # model.zerograds() は古い
        loss = model(x)
        loss.backward()
        optimizer.update()
                                                                                                                        
    # get middle layer data
    x = Variable(xtrain, volatile='on')
    yt = F.sigmoid(model.l1(x))
    ans = yt.data

    return ans
学習が終わったところで、今度はmodelではなく、その中の中間層の値が決まる所まで実行し、中間層のノードが入っているdataだけを取り出して返す。 これで、execonce()を呼ぶと、各入力データに対する中間層の2ノードの値がまとめて帰ってくる。

では、このexeconce() を直接実行してみよう。
>>> ans = execonce()
>>> ans.shape
(150, 2)
>>> ans[:5]
array([[ 0.99554908,  0.00865224],
       [ 0.99200833,  0.01623198],
       [ 0.99302506,  0.01248077],
       [ 0.992621  ,  0.02188915],
       [ 0.99591887,  0.00837815]], dtype=float32)
>>> 
ちゃんとデータ数150個に対し、2つの値が求まっている。
これを、Irisの種類に応じて色分けして表示すると、何か分かるかも知れない。

ということで、最初に示した図が得られるのだが、そのあたりは次回説明しよう。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 前回はAICを使ってモデル選択の様子を見るためのコードを紹介しました。 今回はそのコードの内容を解説します。

このコードを実行すると、4種類ある分布からサンプリングしたデータに対して、4種類のモデルを使ってポアソン回帰をし、AICやバイアスを計算して表示します。 計算されたAICを使って複数あるモデルから最適な物を選択する様子を見ることができます。

回帰に使うデータを作っているのは以下のコードです。 xc0 から 10 までの一様分布をします。 yxccreateY を通して何らかの関係を持つ数です。 xc が説明変数で y が応答変数です。

createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 10)
  c <- runif(dataSize, 0, 10)
  y <- createY(x, c)
  data.frame(x, c, y)
}

関係としては、以下の4つがあります。 xc と関係のある無しで4種類です。

function(x, c) rpois(dataSize, lambda = 4)
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x))
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * c))
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x + 0.2 * c))

ポアソン回帰をやっているのは以下のコードです。 f でリンク関数を指定して、 genData によって生成されたデータ d に対してポアソン回帰をやっています。 推定されたパラメータを使いやすくするために、n という名前を付けています。

calculateAic <- function(d, genData, f, n) {
  fit <- glm(formula = f, family = poisson, data = d)
  co <- append(fit$coefficients, replicate(2, 0))
  names(co) <- n

このようにして得たパラメータを使って、最大対数尤度、平均対数尤度、バイアス、AICを以下のように計算しています。

  meanLogLik <- meanLogLikelihood(genData, co["b1"], co["b2"], co["b3"])
  bias <- logLik(fit) - meanLogLik
  aic <- -2 * (logLik(fit) - length(fit$coefficients))
  result <- c(logLik(fit), meanLogLik, bias, aic)
  names(result) <- c("logLik", "meanLogLik", "bias", "aic")
  result
}

これで1つのモデルに対して、1回分の実験が終わります。 以下のコードで実験を何度も繰り返して、結果を平均しています。

calculateMeanAic <- function(genData, f, n) {
  genData <- match.fun(genData)
  m <- replicate(expSize, calculateAic(genData(), genData, f, n))
  apply(m, 1, mean)
}

以下のコードで4種類のモデルに対して上記の計算をしています。 リンク関数は exp(b1)exp(b1 + b2 * x)exp(b1 + b3 * c)exp(b1 + b2 * x + b3 * c) の4種類です。

showAic <- function(genData) {
  r1 <- calculateMeanAic(genData, y ~ 1, c("b1","b2","b3"))
  r2 <- calculateMeanAic(genData, y ~ x, c("b1","b2","b3"))
  r3 <- calculateMeanAic(genData, y ~ c, c("b1","b3","b2"))
  r4 <- calculateMeanAic(genData, y ~ x + c, c("b1","b2","b3"))

後はサンプリングのための関数をこの関数に渡して、結果を表示しています。 得られた結果の解説は次回に回します。

前回MyAEというクラスを説明した。
このクラスからオブジェクトを作ると、用意ができる。

# Set data
from sklearn import datasets
iris = datasets.load_iris()
xtrain = iris.data.astype(np.float32)

# Initialize model        
model = MyAE()
optimizer = optimizers.SGD()
optimizer.setup(model)
ここまで用意ができたら、ループを回して延々と学習するだけであり、その処理は今までとほとんど同じである。

# Learn
losslist = []
for j in range(5000):
    x = Variable(xtrain)
    model.cleargrads()             # model.zerograds() 非推奨
    loss = model(x)
    if j%1000 == 999:
        print( "%6d   %10.6f" % (j, loss.data) )
    losslist.append(loss.data)     # 誤差をリストに追加
    loss.backward()
    optimizer.update()
今回は、誤差が徐々に縮まっていくところをグラフ化してみよう。
そのため、modelに学習データを与えて実行(学習)させると、誤差(loss)が帰ってきて、それをリストlossに加えていく。
正確には、loss.dataが誤差である。

これで、誤差がlosslistに溜まっているので、表示し、画像ファイルとしてセーブした。

# 誤差グラフの表示
plt.plot(losslist)
plt.savefig("irisautoencodeloss.png")
plt.show()

irisautoencodeloss.png誤差がどんどん小さくなっているいることが分かる。
つまり、学習が進むと、入力と出力がどんどん一致していっていることが分かる。

さて、次回は、もっと突っ込んで、ネットワークの内部のデータを調べてみよう。


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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 そこでAICを使ったモデル選択の様子を見るためのコードを用意しました。 コードはRで書きました。

dataSize <- 100
expSize <- 1000

createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 10)
  c <- runif(dataSize, 0, 10)
  y <- createY(x, c)
  data.frame(x, c, y)
}

createIndipendentData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = 4))
}

createCorrelatedXData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x)))
}

createCorrelatedCData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * c)))
}

createCorrelatedXCData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x + 0.2 * c)))
}

logLikelihood <- function(genData, b1, b2, b3) {
  genData <- match.fun(genData)
  data <- genData()
  prob <- vector(length = dataSize)
  for (j in 1:dataSize) {
    prob[j] <- log(dpois(data$y[j], lambda = exp(b1 + b2 * data$x[j] + b3 * data$c[j])))
  }
  sum(prob)
}

meanLogLikelihood <- function(genData, b1, b2, b3) {
  mean(replicate(dataSize, logLikelihood(genData, b1, b2, b3)))
}

calculateAic <- function(d, genData, f, n) {
  fit <- glm(formula = f, family = poisson, data = d)
  co <- append(fit$coefficients, replicate(2, 0))
  names(co) <- n

  meanLogLik <- meanLogLikelihood(genData, co["b1"], co["b2"], co["b3"])
  bias <- logLik(fit) - meanLogLik
  aic <- -2 * (logLik(fit) - length(fit$coefficients))
  result <- c(logLik(fit), meanLogLik, bias, aic)
  names(result) <- c("logLik", "meanLogLik", "bias", "aic")
  result
}

calculateMeanAic <- function(genData, f, n) {
  genData <- match.fun(genData)
  m <- replicate(expSize, calculateAic(genData(), genData, f, n))
  apply(m, 1, mean)
}

showAic <- function(genData) {
  r1 <- calculateMeanAic(genData, y ~ 1, c("b1","b2","b3"))
  r2 <- calculateMeanAic(genData, y ~ x, c("b1","b2","b3"))
  r3 <- calculateMeanAic(genData, y ~ c, c("b1","b3","b2"))
  r4 <- calculateMeanAic(genData, y ~ x + c, c("b1","b2","b3"))

  cat("\t\t\tnull\t\tx\t\tc\t\txc\n")
  cat("log likelihood\t\t", r1["logLik"], "\t", r2["logLik"], "\t",
    r3["logLik"], "\t", r4["logLik"], "\n")
  cat("mean log likelihood\t", r1["meanLogLik"], "\t",
     r2["meanLogLik"], "\t", r3["meanLogLik"], "\t", r4["meanLogLik"], "\n")
  cat("bias\t\t\t", r1["bias"], "\t", r2["bias"], "\t", r3["bias"], "\t", r4["bias"], "\n")
  cat("aic\t\t\t", r1["aic"], "\t", r2["aic"], "\t", r3["aic"], "\t", r4["aic"], "\n\n")
}

cat("indipendent data\n")
showAic(createIndipendentData)

cat("correlated X data\n")
showAic(createCorrelatedXData)

cat("correlated C data\n")
showAic(createCorrelatedCData)

cat("correlated XC data\n")
showAic(createCorrelatedXCData)

少し長くなったので、解説は次回に回します。 ぜひ実行してみてください。

Deep Learningでクラス分けをやってきたが、今回はちょっと違った処理をしてみる。
クラス分けというのは、よく分からないデータを分類するもので、Irisの場合は4種のデータからIrisの種類を推測するものであり、手書き数字の場合は数字の画像データから数字を当てるものであった。

今回は、そういうのとは全然違う AutoEncoder というものを作ってみる。
今回は、簡単に済ますために、またIrisのデータを使うことにする。
すると、入力は4つの値の組である。
クラス分けでは、出力はIrisの種類であったのだが、今回は出力は入力と同じものになるように指示してみよう。つまり、入力と出力が同じである。

ネットワーク図いうと、入力と出力は4ノードである。
しかし、中間層は、2ノードの層を1つだけにする。
つまり、次図のような構成だ。

AutoEncoder.pngつまり、4つのデータを、途中で2つのデータに変換(圧縮)し、その2つのデータから元のデータを復元し、4つのデータにすることを考える。
これを実現するのに、いままで説明してきた仕組みをそのまま使う。
違うところは、誤差を最小化するところで、今回は、入力がこのネットワークで変換され出力された値と、元の値との2乗誤差を求め、その誤差が最小化するようにする。

とりあえず、今回は、誤差がどんどん小さくなっていることだけを確認してみよう。
class MyAE(Chain):
    def __init__(self):
        super(MyAE, self).__init__(
            l1=L.Linear(4,2),               # 中間層は2ノード    
            l2=L.Linear(2,4),
        )
        
    def __call__(self,x):
        bv = self.fwd(x)
        return F.mean_squared_error(bv, x)  # 結果と元データとの2乗誤差
        
    def fwd(self,x):
        fv = F.sigmoid(self.l1(x))
        bv = self.l2(fv)
        return bv

これを利用して実行すると、次第に2乗誤差が小さくなる。
Chainer$ python -i ae1.py 
   999     0.602435
  1999     0.340246
  2999     0.208172
  3999     0.156505
  4999     0.136456

今回は、数値による結果だけを示し、次回にプログラムも含めて説明する。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 それを理解するために、ポアソン回帰のバイアスを見るためのコードを用意しました。 コードはRで書きました。

createData <- function() {
  y <- rpois(10, lambda = 4)
  data.frame(y)
}

logLikelihood <- function(data, b1) {
  sum(log(dpois(data, lambda = exp(b1))))
}

printLogLikelihood <- function(data, i) {
  cat(paste0("Input data is data", i, "\n"))
  print(data[i]$y)

  fit <- glm(formula = y ~ 1, family = poisson, data = data[i])
  b1 <- fit$coefficients[1]

  cat("Estimated b1 =", b1, "\n")
  cat("Log likelihood for\n")
  cat("data1\t\tdata2\t\tdata3\t\n")
  cat(logLikelihood(data[1]$y, b1), "\t",
    logLikelihood(data[2]$y, b1), "\t",
    logLikelihood(data[3]$y, b1), "\n\n")
}

printData <- function(data, i) {
  cat(paste0("data", i))
  print(data[i]$y)
}

data <- c(createData(), createData(), createData())

for(i in 1:3) {
  printData(data, i)
}
cat("\n")
for(i in 1:3) {
  printLogLikelihood(data, i)
}

コードを実行すると、まずポアソン分布からサンプリングした 10 個のデータからなるデータセットを3個作ります。 それぞれのデータに対してポアソン回帰をして、推定されたパラメータを使って他の2つのデータセットに対して対数尤度を計算します。 対数尤度を計算しているコードは以下のコードです。

logLikelihood <- function(data, b1) {
  sum(log(dpois(data, lambda = exp(b1))))
}

b1 が最尤推定されたパラメータです。 data が対数尤度を計算する対象のデータです。 平均 exp(b1) を持つポアソン分布からデータが得られる確率の対数を計算しています。 実行結果は以下のようになります。

data1 [1] 5 1 2 4 9 3 3 2 3 4
data2 [1] 7 4 3 3 3 8 4 6 3 1
data3 [1] 9 1 5 6 1 3 3 8 3 2

Input data is data1
 [1] 5 1 2 4 9 3 3 2 3 4
Estimated b1 = 1.280934 
Log likelihood for
data1       data2       data3   
-20.59338    -21.43294   -24.32331 

Input data is data2
 [1] 7 4 3 3 3 8 4 6 3 1
Estimated b1 = 1.435085 
Log likelihood for
data1       data2       data3   
-21.04396    -20.95861   -24.00313 

Input data is data3
 [1] 9 1 5 6 1 3 3 8 3 2
Estimated b1 = 1.410987 
Log likelihood for
data1       data2       data3   
-20.91147    -20.97071   -23.99113 

傾向として、ポアソン回帰に使ったデータセットへの尤度が一番大きくて、他のデータセットへの尤度は小さくなります。 時々は、他のデータセットへの尤度が大きくなります。 最尤推定は与えられたデータに最もよく当てはまるパラメータを返すので、まだ見ぬデータにもよく当てはまるかどうかは分からないということです。

LearningPython5th.gif
Learning Python, 5th Edition
Powerful Object-Oriented ProgrammingDSCF5158 (400x300).jpg
By Mark Lutz
Publisher: O'Reilly Media
Final Release Date: June 2013
Pages: 1648

ISBN-13: 978-1449355739

厚さ: 57mm



Pythonの本で、1冊非常に分厚いのがある。厚さ57mm、重さ不明、ページ数は本文だけで1417ページもある。
それでいて、本書はタイトルにLearning とついているように、入門書なのだ。
持ち歩きたい場合には、紙ではなく、電子ブック版を入手した方が良いであろう。すくなくとも、持ち歩ける。

本書の和訳も出ているのだが、本書と比べると、とても薄い。
理由は、和訳は3th Editionなのだ。これでは、ちょっと躊躇してしまう。Python 3 対応ではないので、オススメできない。

本書、これだけ分厚くても、NumPyなどの人工知能などでよく使う拡張モジュールについては、存在することが紹介されているだけで、その説明はないのだ。
つまり、本書は、Pythonの基本部分だけを、延々と丁寧に説明している。

本体は、8 パート, 41章から成り立っている。1章は30-40ページ程度と、かなりコンパクトになっていて、こういう単位で読みきるには丁度良いかもしれない。

普通のPythonの本にはちょっとしか載っていないものがある。
たとえば、comprehension(内包)、lambda expression(λ式)、string formatting などがそうだ。
その他、他の本で説明が見つからなかった場合に、結構見つかることがある。

でも、あまりにページ数が多く、丁寧過ぎるきらいがあり、延々と読み進めるには退屈に感じる。

ということで、傍らに置いておき、時々参照するのに適した本である。

もちろん、全部英語で書かれているのだが、平易な英語で、やさしく説明してあるので(本書は入門書)、本書を延々と読むことは、英語とPythonを両方同時に習得できて一石二鳥であろう。

このアーカイブについて

このページには、2017年3月に書かれたブログ記事が新しい順に公開されています。

前のアーカイブは2017年2月です。

次のアーカイブは2017年4月です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。