TIM Labs

2017年4月アーカイブ

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

この章では様々なタイプのGLMについて説明がされています。 その中でガンマ回帰について説明がされています。 なので、いろいろなパラメータについてガンマ分布をプロットしてくれるコードを用意しました。 コードはRで書きました。

rs <- seq(0.5, 4.5, by = 1.0)
ys <- seq(0.5, 19.5, by = 1.0)
xl <- "y"
yl <- "Probability"
legends <- paste0("r = ", rs)
for (s in 1:20) {
  title <- paste0("Gamma distribution, p(y|shape = ", s, ", rate = r)")
  plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
  for (i in 1:5) {
    lines(ys, dgamma(ys, s, rs[i]), type = "l")
    points(ys, dgamma(ys, s, rs[i]), pch = i)
  }
  legend("topright", legend = legends, pch = 1:5)
}

実行するとグラフが20枚プロットされます。 例えば下図のようなグラフがプロットされます。

gamma_4.png

ガンマ分布はパラメータの組によっては、指数分布のような分布や、ポアソン分布のような山なりの分布を取れることが分かります。

MNISTは、グレースケールの手書き数字のデータ・セットだった。
でも、もう飽きたので、写真などでテストすることにしよう。

といっても、自分でディープラーニング用の画像データを集めるのは気が遠くなるほど面倒だ。

それで、まず、ディープラーニング用のデータ・セットでどのようなものがあるか探ってみた。

http://deeplearning.net/datasets/

そうすると、こちらが探そうとしていることが、そのままURLになっているのが見つかった(笑)

These datasets can be used for benchmarking deep learning algorithms:

という文の下に、自由に使えそうなデータセットが並んでいるのだ。
MNIST関連は飛ばして、色々な写真などのデータセットらしいのを探そう。

CIFAR-10
これは、元トロント大、現在GoogleのAlex Krizhevsky氏が配布しているものだ。
サイトは、今もトロント大にあり、データ収集は
Alex Krizhevsky, Vinod Nair, and Geoffrey Hinton の3名が行ったとある。
最後の Geoffrey Hinton
は、the godfather of deep learning として知られ、Google/DeepMind の AlphaGo にも関係している人なのだ。
画像は32x32とかなりコンパクトで、学習用に50000枚、テスト用に10000枚用意されている。
これは外せないな。

Caltech 101
こちらは、その名の通り、カルテック、カリフォルニア工科大学が用意してくれている画像データセットだ。
101種類の画像で、サイズが300x200程度という。各種類40から500枚くらいで、かなりバラバラである。

その他にも色々あるので、あとは自分で見てみよう。

さて、どれを利用してみるか、それが問題だ。

mnistのデータセットは、こんな感じで読み込んだのだった。

train, test = chainer.datasets.get_mnist()
次のデータセットも、同じように読み込めれば楽である。
ということで、こんな感じに読み込めそうな、つまりchainerがサポートしているデータセットを探した。

Docs ≫ Chainer Reference Manual ≫ Dataset examples

この中が、General datasets と Concrete datasets に別れており、Concrete datasets の中にmnistがあったのだ。
mnist以外では、CIFAR10/100 と Pen Tree Bank があったが、後者は英文に関するデータセットで今考えているものとは違う。

ということで、すんなり、CIFARを選ぶことに決定した。
次回から、実際にデータを読み込んで、あれこれやってみよう。

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

この章では様々なタイプのGLMについて説明がされています。 その中でオフセット項わざという手法について説明がされているので、サンプルデータに対してこの手法を試してみるコードを用意しました。 コードはRで書きました。

dataSize <- 100
createData <- function(n, createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 5)
  A <- runif(dataSize, n, 10 * n)
  y <- createY(x, A)
  data.frame(x, A, y)
}

linkFunction <- function(x, A) rpois(dataSize, exp(0.1 + 0.3 * x + log(A)))
dataGenerator1 <- function() createData(1, linkFunction)
dataGenerator2 <- function() createData(100, linkFunction)

doGlm <- function(dataGen) {
  dataGen <- match.fun(dataGen)
  data <- dataGen()
  fit <- glm(y ~ x, offset = log(A), data = data, family = poisson)
  c(fit$coefficients[1], fit$coefficients[2])
}

showResult <- function(dataGen) {
  dataGen <- match.fun(dataGen)
  m <- replicate(100, doGlm(dataGen))
  meanm <- apply(m, 1, mean)
  sdm <- apply(m, 1, sd)
  cat("mean b1", meanm[1], "\n")
  cat("mean b2", meanm[2], "\n")
  cat("sd b1", sdm[1], "\n")
  cat("sd b2", sdm[2], "\n")
}

cat("case1\n")
showResult(dataGenerator1)
cat("\ncase2\n")
showResult(dataGenerator2)

サンプルデータに含まれる変数は xAy です。 説明変数は xA です。 応答変数は y です。 ポアソン回帰のリンク関数は λ = exp(b1 + b2 * x + log(A)) という形です。 なので、log(A) 部分がオフセット項になります。 本の例で言うと A は面積のような、y との割算で密度を計算できるような値です。

本の中では1000打数300安打の打者と10打数3安打の打者を同じ三割打者と見なしてよいかという問題に触れられています。 サンプルデータの変数名で言うと、A = 1000y = 300 の場合と、A = 10y = 3 の場合とを同じと見なしてよいかという問題です。 感覚的には、10回中の3回打つ打者より、1000回中の300回打つ打者の方がより正確に三割打者であると言えます。 なので、これらの状況には何らかの差が出てきて欲しいと思うところです。

この問題の状況を見るために、ポアソン回帰を2つの場合について行いました。 ケース1においては A1 から 10 までのランダムな値を取ります。 ケース2においては A100 から 1000 までのランダムな値を取ります。 x はどちらの場合も 0 から 5 までのランダムな値です。 yλ = exp(0.1 + 0.3 * x + log(A)) つまり λ = A * exp(0.1 + 0.3 * x) というリンク関数で得られる λ を平均に持つポアソン分布からサンプリングされます。 なので、A が100倍になれば、y の平均も100倍になります。

コードを実行すると次のような結果が得られます。

case1
mean b1 0.09576121 
mean b2 0.3013721 
sd b1 0.05874369 
sd b2 0.01892569 

case2
mean b1 0.1002007 
mean b2 0.299965 
sd b1 0.006423389 
sd b2 0.00181942 

推定されるパラメータの平均値はどちらも大差がありませんが、その標準偏差は10倍近い差があります。 確かに A の大きさの違いが推定値のばらつきに影響を与えているのが分かります。

ChainerでGPUを使ってみようで、MNISTのデータセットを使ったディープラーニングのサンプルプログラムである train_mnist.py のコマンドで、-gpu 0 でGPUが使えることを示した。
つまり、操作説明をしたに過ぎなかったが、今回は、プログラムレベルでどう対応しているか見てみよう。
プログラムは,https://github.com/pfnet/chainer/tree/master/examples/mnist である。

    parser.add_argument('--gpu', '-g', type=int, default=-1,
                        help='GPU ID (negative value indicates CPU)')
ここで、 -gpu の処理が行われ、args.gpu に整数値が入る。
デフォルト値は -1 で、GPUを使用しないことを示す。

    if args.gpu >= 0:
        chainer.cuda.get_device(args.gpu).use()  # Make a specified GPU current
        model.to_gpu()  # Copy the model to the GPU
GPUを使用するときだけ、指定したGPUを確保し、作ったモデルをGPUにコピーするようだ。

    # Set up a trainer
    updater = training.StandardUpdater(train_iter, optimizer, device=args.gpu)
    trainer = training.Trainer(updater, (args.epoch, 'epoch'), out=args.out)

    # Evaluate the model with the test dataset for each epoch
    trainer.extend(extensions.Evaluator(test_iter, model, device=args.gpu))

あとは、上記の2個所で、device指定があるので、args.gpu をそのまま与えている。
GPUのの使用に関係してくるのは、たったこれだけである。
 あとは、chainerがよきに計らってくれるので、何もすることはない。
フレームワークとは、こういうものである。

ということで、よくわからなくてもGPUを使うプログラムが書けてしまう。
便利というか、便利過ぎな気がする。

そのうち、Python から直接GPUいじるためにCuPyを紹介するかもしれない。 NumPy互換のインターフェイスでCUDA(NVIDIAのGPUライブラリ)を使えるらしい。
デフォルのままだと、中間層が2層で、それぞれ1000ユニットで、
784 ==> 1000 ==> 1000 ==> 10
として動かしていた。

中間層の2層のノード数を自由に変更できるようになったので、まず次のようにしてみた。
784 ==> 300 ==> 30 ==> 10

Chainer$ python train_mnist1.py -g 0 -u 300 30
GPU: 0
# number: 60000
# units: [300, 30]
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           0.307806    0.150247              0.913467       0.954                     2.44043       
2           0.123374    0.112284              0.963534       0.9662                    4.77677       
3           0.0812434   0.0797795             0.975833       0.9758                    7.09853  
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・                  
18          0.00678134  0.0989654             0.997766       0.9799                    42.7799       
19          0.00394047  0.107429              0.998866       0.9782                    45.2334       
20          0.00806017  0.0969633             0.997466       0.9803                    47.6298       
Chainer$ 
かなり中間層を小さくしたのだが、テスト精度は98%あり、精度は落ちていない。
つまり、中間層を小さくしても影響は少ないようだ。

ということで、どんどん減らしてみよう。
中間層の第2層を出力層と同じ10まで下げてみた。
分類数と同じだが、どうなるだろうか?

784 ==> 50 ==> 10 ==> 10

Chainer$ python train_mnist1.py -g 0 -u 50 10
GPU: 0
# number: 60000
# units: [50, 10]
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           0.583669    0.312933              0.828284       0.9087                    2.45214       
2           0.26503     0.23423               0.925767       0.931                     4.88362       
3           0.202712    0.191492              0.942717       0.9448                    7.3304        
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・                  
18          0.0472231   0.132507              0.985798       0.9636                    42.9347       
19          0.0456997   0.120474              0.986715       0.9653                    45.3106       
20          0.0433583   0.123629              0.986148       0.9672                    47.6485       
Chainer$ 

かなり中間層を絞ったのだが、まだ96%の正解率である。
ということで、思いっきり絞ってみた。
784 ==> 40 ==> 2 ==> 10

Chainer$ python train_mnist1.py -g 0 -u 40 2
GPU: 0
# number: 60000
# units: [40, 2]
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           1.66372     1.36298               0.392867       0.5006                    2.45188       
2           1.20082     1.08511               0.562367       0.6294                    4.81112       
3           0.997415    0.951765              0.65775        0.6989                    7.13735       

18          0.3529      0.502587              0.915784       0.8996                    42.5103       
19          0.341643    0.48146               0.918735       0.9035                    44.8852       
20          0.329597    0.491561              0.920851       0.9042                    47.2294       
Chainer$ 
という訳で、中間層のユニット数をどんどん減らし、第2層のユニット数を2まで減らしても正解率90%を維持しているのは脅威だ。

といっても、手書き数字の読み取りでの場合なのだ。
少々手書き数字には飽きてきたので、別のデータセットなども考えてみることにする。

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

この章では様々なタイプのGLMについて説明がされています。 その中で、ロジスティック回帰についてAICを使ってモデル選択をする様子が説明されています。 なので、サンプルデータを用意して、いろいろなモデルについてAICを計算してくれるコードを用意しました。 コードはRで書きました。

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

dataGenerator <- function()
  createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.5 * x))))

data <- dataGenerator()

calcAIC <- function(formula) {
  fit <- glm(formula, data = data, family = binomial)
  fit$aic
}

cat("model\t\taic\n")
cat("null\t\t", calcAIC(cbind(y, N - y) ~ 1), "\n")
cat("x\t\t", calcAIC(cbind(y, N - y) ~ x), "\n")
cat("c\t\t", calcAIC(cbind(y, N - y) ~ c), "\n")
cat("x + c\t\t", calcAIC(cbind(y, N - y) ~ x + c), "\n")
cat("xc\t\t", calcAIC(cbind(y, N - y) ~ x : c), "\n")
cat("x + xc\t\t", calcAIC(cbind(y, N - y) ~ x + x : c), "\n")
cat("c + xc\t\t", calcAIC(cbind(y, N - y) ~ c + x : c), "\n")
cat("x + c + xc\t", calcAIC(cbind(y, N - y) ~ x * c), "\n")

サンプルデータは xcy を変数として含んでいます。 説明変数は x で応答変数は y です。 そこに y と相関のない変数として c を用意しました。 c は必要ない変数なので、モデルに組み込むとバイアスが大きくなるはずです。 コードを実行すると以下のような結果が出力されます。

model       aic
null         885.6348 
x            264.6336 
c            887.2804 
x + c        265.7918 
xc           667.2655 
x + xc       265.8728 
c + xc       381.9791 
x + c + xc   267.7918

確かに x だけがモデルに使われている時に一番AICが低くなるのが分かります。

サンプルプログラム train_mnist.py のネットワークの定義は以下のようになっていた。

    def __init__(self, n_units, n_out):
        super(MLP, self).__init__(
            # the size of the inputs to each layer will be inferred
            l1=L.Linear(None, n_units),  # n_in -> n_units
            l2=L.Linear(None, n_units),  # n_units -> n_units
            l3=L.Linear(None, n_out),  # n_units -> n_out
        )
このプログラムでは、入力層のサイズは指定がなく、データにより自動決定される。
2つの中間層はいずれも n_units でノード数が与えられるので、同じサイズになる。
デフォルトでは1000である。
出力層は、0から9を判別するためにn_out に10がパラメータとして与えられる。

中間層の2層が同じサイズというのは、ちょっと面白くない。
それぞれサイズを変えることができれば、中間層のノード数の変化と学習の関係をより詳しく調べられる。

デフォルトのままだと 784 ==> 1000 ==> 1000 ==> 10 であるのを、 784 ==> 150 ==> 50 ==> 10 という感じの指定ができるようにしたい。

コマンドラインパーサーの説明は、先日行った。
その中間層のサイズの解析部分は以下のようになっている。

    parser.add_argument('--unit', '-u', type=int, default=1000,
                        help='Number of units')
これをどう変更すればよいだろうか?

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

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、サンプルデータを用意して、ロジスティック回帰をやってみるコードを用意しました。 コードはRで書きました。

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

dataGenerator <- function()
  createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.05 * x))))

data <- dataGenerator()
fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

linkFunction <- function(x) 1 / (1 + exp(-b1 - b2 * x))

xs <- 1:200
ys <- 0:N
title <- "Logistic regression"
xl <- "x"
yl <- "The most probable y"
maxY <- numeric(length = 200)
for (x in xs) {
  q <- linkFunction(x)
  maxY[x] <- which.max(dbinom(ys, N, q)) - 1
}
plot(0, 0, type = "n", xlim = c(0, 200), ylim = c(0, N), main = title, xlab = xl, ylab = yl)
lines(xs, maxY, type="l")

前回もサンプルデータに対してロジスティック回帰をやってみて結果をプロットしました。 今回は説明変数 x に対して、どの応答変数 y が一番高い確率で得られるかをプロットしました。 コードを実行すると下図のような図がプロットされます。

logistic_regresson_2.png

ロジスティック関数のような曲線が得られました。 もちろん y の値はこの曲線の周囲に分布しているわけですが、一番取りやすい値をつなげると、このような曲線になるということです。

実装DL (285x400).jpg

題名:実装ディープラーニング

株式会社フォワードネットワーク 監修

藤田一弥+高橋歩 共著

A5、272項、本体3200円

2016/11/30 発行

オーム社

ISBN 978-4-274-21999-3


本書は、題名どおり「実装」について説明した本である。
実装といっても、ゼロから実装ではなくて、世界的にも有名なデータや、ニューラルネットワークを自分のマシンで走らせる方法についての手引書といったところ。
ディープラーニングに関する基本的な知識はざっとは書いているが、ちゃんと勉強するには他書を読んだほうが良い。

使っているフレームワークなどは、適当に色々使っている感じ。

Theano, Keras, Chainer などをインストールした上で行う。
もちろん、Numpyを初めとしたPythonの人工知能系でよく使われる各種ライブラリもインストールする。
そして、そうなると問題になるのが各ライブラリのバージョンの整合性の問題に悩むことになる。
本書では、使用した各ライブラリについてバージョンが示されている。
NumPyについては、ダウングレードして難を逃れたと書かれていた。
現状で、PythonでAIをやろうと思ったら、こういうことに慣れるしかない。

また、一部はCのプログラムも紹介されている。

それよりも、GPUの利用について、本書の最初に説明がされている。
ただし、GPUは変化が激しい世界なので、かならずネットで調べるなり、すでにGPUを利用している人に聞くなりするべきだ。

10層、あるいはそれ以上のモデルでの実験を行うので、GPUを前提とした内容と思ったほうが良く、まだGPUを入手していない場合は、本書の最初のGPUの解説や、ネットの情報、NVIDIAの情報などを参考に、NVIDIAのGPUを入手しよう。

テーマは、ほとんどが画像処理ばかりという感じの本だが、最後の方に、3目並べを強化学習で強くする方法について解説しているようだ。

使用する学習データもサイズが大きくなり、どれも数百MBくらいはある。
まとめてダウンロードしようとすると何GBにもなるので、回線が遅い場合は何らかの手段を考えよう。

ディープラーニングにおいて、極度に単純な入門例の解説を読み終えて、もうちょっと本格的なものに手を出してみたい人には丁度お手ごろな本だろう。

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

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、サンプルデータを用意して、ロジスティック回帰をやってみるコードを用意しました。 コードはRで書きました。

N <- 10
dataSize <- 100

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

dataGenerator <- function()
  createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.5 * x))))

data <- dataGenerator()

fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

linkFunction <- function(x) 1 / (1 + exp(-b1 - b2 * x))

xs <- seq(5, 20, by = 3)
ys <- 0:N
title <- "Logistic regression"
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
for (x in xs) {
  q <- linkFunction(x)
  y <- dbinom(ys, N, q)
  lines(ys, y, type="l")
  points(ys, y, pch = x)
}
legend("topleft", legend = legends, pch = xs)

作ったサンプルデータにロジスティック回帰をやって、リンク関数のパラメータを推定しているのは以下の部分です。

fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

こうして得られた b1b2 を使って以下のリンク関数を作れば、説明変数 x に対して、事象の発生確率 q が得られます。

linkFunction <- function(x) 1 / (1 + exp(-b1 - b2 * x))

これで事象を N 回繰り返した時の応答変数 y の分布が求まるので、あとはそれをプロットします。 コードを実行すると、下図のようなグラフがプロットされます。

logistic-regression.png

推定される b2 は正の数なので、x が大きくなるにしたがって y の分布も大きな値を取るようになっています。

学習の進行状況を把握するにはグラフがなんと言ってもわかりやすい。
仕方がないから、ちょっとプログラムして、グラフを作るしかないかな。

と思ったら、resultというフォルダができていて、その中に何やら色々なものができていた。
その中に、accuracy.png と loss.png があり、こんな感じになっていた。

accuracy.png      -n 60000
MNIST6000accuracy.png 青線が学習精度曲線であり、オレンジ線がテスト精度曲線である。
学習曲線は、1.00にかなり近いところまでいっている。
テスト曲線(検証曲線)との差はグラフからはかなりあるように見えるが、0.015(1.5%)程度である。

もう1つのグラフがロス関数(誤差関数)というものであり、上の図をちょうど上下反転したような形になっている。

loss.png      -n 60000
MNIST60000loss.png
さて、データ数が少なくなったら学習曲線がどんな風に変わったか、n=500の場合の精度の変化の様子を見てみよう。

accuracy.png      -n 500
MNIST500accuracy.png今度は、学習精度が早い段階で1.00になってしまっている。
これだけ見ると、非常に学習が順調に進んだように見える。

しかし、テスト精度の方は、0.85あたりでふらふらしているので、学習とテストでは精度が0.15(15%)も違い、エポックが進んでも差は縮まらないようだ。
学習データが60000個の時の10倍もの誤差になってしまった。

こういう風に、学習精度だけ上がって、テスト精度との差が縮まらない状態を過学習という。
つまり、学習のデータを学習し過ぎて、学習データは完璧に判定できるようになったものの、本番データを入れたら精度が悪いということである。
これを避けるためには、大量のデータを使って学習しないとダメということだ。


それにしても、こんなグラフまで何もしないで作ってくれるchainerのサンプルプログラムは素晴らしい。

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

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、いろいろなパラメータについてロジスティック回帰のモデルをプロットしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

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

distribution <- function(x, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dbinom(y, N, l(x))
  }
}

xs <- seq(0, 10, by = 2)
ys <- 0:N
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
for (b1 in seq(-0.5, 0.5, by = 0.2)) {
  for (b2 in seq(-0.5, 0.5, by = 0.2)) {
    l <- linkFunction(b1, b2)
    title <- paste0("q = 1 / (1 + exp(-(", b1, " + ", b2, "x)))")
    plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
    for (x in xs) {
      d <- distribution(x, l)
      lines(ys, d(ys), type = "l")
      points(ys, d(ys), pch = x)
    }
    legend("topright", legend = legends, pch = xs)
  }
}

リンク関数はロジスティック関数です。 コードでは以下の部分です。

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

これが二項分布のパラメータの一つである事象の発生確率 q を与えます。 もう一つのパラメータは試行回数 N ですが、これは 20 にしました。 ということで、リンク関数のパラメータ b1b2 さえ与えれば、説明変数 x に対して応答変数 y の分布がどのように変化するかが分かります。

コードを実行すると様々な b1b2 の値に対してモデルをプロットします。 例えば下図のような図がプロットされます。

binomial_model.png

これは b1 = -0.3b2 = 0.3 の場合です。 x0 ならq = 1 / (1 + exp(0.3)) = 0.425... です。 x10 ならq = 1 / (1 + exp(-2.7)) = 0.937... です。 b2 が正なので x の増加に従って exp の引数部分がどんどん小さくなって、q 全体としては 1 に近づくわけです。 事象の発生確率が高くなれば、N 回の試行でその事象が起きる回数は大きくなります。 図を見ると、確かに x の増加に従って y の分布が大きな値を取るようになっているのが分かります。

学習データ数を減らして実験をしてみよう。
オリジナルは、MNISTの学習用データ数の60000であるが、10000, 2000, 500と減らしながら実験してみた。
もちろん、時間がもったいないので GPU を利用した。

学習データ 60000個
Chainer$ python train_mnist0.py -g 0
GPU: 0
# number: 60000
# unit: 1000
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           0.193926    0.0882014             0.941183       0.9735                    2.43387       
2           0.075335    0.0752461             0.976834       0.976                     4.83063       
3           0.0480114   0.074878              0.984965       0.978                     7.23542       
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・             
18          0.00998753  0.104438              0.996899       0.9828                    43.1969       
19          0.00818152  0.0978973             0.997833       0.9821                    45.6029       
20          0.00540837  0.0948488             0.998149       0.9817                    48.0003       
Chainer$
学習データ 10000個
Chainer$ python train_mnist0.py -g 0 -n 10000
GPU: 0
# number: 10000
# unit: 1000
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           0.4188      0.217062              0.8749         0.9334                    0.925599      
2           0.153328    0.196418              0.9534         0.9377                    1.74652       
3           0.0852114   0.139428              0.9756         0.9561                    2.58325  
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・             
18          7.43423e-05  0.158279              1              0.9674                    14.9663       
19          6.48571e-05  0.159226              1              0.9674                    15.7797       
20          5.74103e-05  0.160547              1              0.9673                    16.6219       
Chainer$
学習データ 2000個
Chainer$ python train_mnist0.py -g 0 -n 2000
GPU: 0
# number: 2000
# unit: 1000
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           0.949915    0.49678               0.7335         0.8402                    0.648741      
2           0.315252    0.36667               0.9065         0.8894                    1.21073       
3           0.16682     0.316182              0.957          0.9049                    1.80692       
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・             
18          0.00063135  0.327181              1              0.9257                    10.3965       
19          0.000569121  0.328971              1              0.925                     10.9622       
20          0.000512893  0.331511              1              0.9251                    11.5304       
Chainer$
学習データ 500個
Chainer$ python train_mnist0.py -g 0 -n 500
GPU: 0
# number: 500
# unit: 1000
# Minibatch-size: 100
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           1.814       1.16017               0.456          0.7062                    0.596051      
2           0.66647     0.718226              0.85           0.7652                    1.12284       
3           0.345737    0.635446              0.888          0.7999                    1.65856       
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・             
18          0.00110594  0.596009              1              0.863                     9.70698       
19          0.00101897  0.600444              1              0.8626                    10.2394       
20          0.000941102  0.6048                1              0.8626                    10.7713       
Chainer$
学習データが60000個から500個に減っても、速度は4倍程度しか高速になっていない。
やはり、GPUを使う場合、色々オーバーヘッドが多いということであろうか。

いや、それよりも、テスト時の正解率の減少の方が問題だ。
98%を超えていた正解率が、86%で終わっている。

そうだ、正解率など、学習経過を延々と数字で見るのではなく、パッと一瞬で分かる方法が欲しいものだ。
どうすれば良いだろうか?

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

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 ロジスティック回帰で使われる分布は二項分布なので、いろいろなパラメータについて二項分布をプロットしてくれるコードを用意しました。 コードはRで書きました。

qs <- seq(0.1, 0.9, by = 0.2)
xl <- "y"
yl <- "Probability"
legends <- paste0("q = ", qs)
for (N in 1:20) {
  title <- paste0("Binomial distribution, p(y|N = ", N, ", q)")
  plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
  for (i in 1:5) {
    lines(0:N, dbinom(0:N, N, qs[i]), type = "l")
    points(0:N, dbinom(0:N, N, qs[i]), pch = i)
  }
  legend("topright", legend = legends, pch = 1:5)
}

実行するとグラフが20枚プロットされます。 例えば下図のようなグラフがプロットされます。

binomial_15.png

二項分布は確率 q で起きる事象が N 回中 y 回起きる確率の分布です。 この図は N = 15 の場合です。 q が小さいうちは事象は少ない回数起きる確率が高く、 q が大きくなるにしたがって多くの回数起きる確率が高くなるのが分かります。

MNISTの学習データは6万画像もあって、非常に多い。
もっと少ない枚数でも大丈夫かも。
ということで、学習データ画像数とテストの正解率の関係を知りたいと思ったのだが、コマンド引数にはデータ数を指定する項目がない。

それで、サンプルソース(train_mnist.py)プログラムを見たら、コマンド引数の処理をargparseで行っているようだ。
main()の最初の部分を以下に示す。

import argparse

def main():
    parser = argparse.ArgumentParser(description='Chainer example: MNIST')
    parser.add_argument('--batchsize', '-b', type=int, default=100,
                        help='Number of images in each mini-batch')
    parser.add_argument('--epoch', '-e', type=int, default=20,
                        help='Number of sweeps over the dataset to train')
    parser.add_argument('--frequency', '-f', type=int, default=-1,
                        help='Frequency of taking a snapshot')
    parser.add_argument('--gpu', '-g', type=int, default=-1,
                        help='GPU ID (negative value indicates CPU)')
    parser.add_argument('--out', '-o', default='result',
                        help='Directory to output the result')
    parser.add_argument('--resume', '-r', default='',
                        help='Resume the training from snapshot')
    parser.add_argument('--unit', '-u', type=int, default=1000,
                        help='Number of units')
    args = parser.parse_args()

    print('GPU: {}'.format(args.gpu))
    print('# unit: {}'.format(args.unit))
    print('# Minibatch-size: {}'.format(args.batchsize))
    print('# epoch: {}'.format(args.epoch))
    print('')
ということで、Pythonのマニュアルを調べたら、説明があった。
16.4. argparse -- コマンドラインオプション、引数、サブコマンドのパーサー

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

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 対立仮説が正しい場合でも、尤度比検定は必ずしも帰無仮説を棄却できるわけではありません。 そこで、パラメトリックブートストラップ法で帰無仮説を棄却できない場合がどれくらい現れるかを見るためのコードを用意しました。 コードはRで書きました。

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

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
}

test <- function(n, b2) {
  expData <- createData(n, 0.1, b2)
  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(100, bootstrap(fit1$coefficient[1]) > dev))
  P > 0.05
}

xl <- "Data size"
yl <- "The ratio of failing to reject"
b2 <- c(0.05, 0.07, 0.1, 0.12)
legends <- paste0("b2 = ", b2)
plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0, 1), xlab = xl, ylab = yl)
n <- 1:20
for(i in 1:4) {
  r <- vector(length = 20)
  for (n in 1:20) {
    r[n] <- mean(replicate(100, test(n, b2[i])))
  }
  lines(1:20, r, type = "b", pch = i) 
}
legend("topright", legend = legends, pch = 1:4)

実行するとグラフを1枚プロットします。 最尤推定を 202 * 100 * 20 * 4 回やるので結構時間がかかります。

有意水準としては5%を設定しました。 なので、パラメトリックブートストラップ法を使ってP値を計算し、得られたP値が 0.05 より大きいと帰無仮説を棄却できません。 test は検定に使うデータサイズ n と、データをサンプリングする分布の中の説明変数 x の係数 b2 を受け取り、P値を計算して帰無仮説が棄却されなかった場合 true を返します。 あとは、各 nb2 について、それを100回繰り返して、true が帰ってくる確率を求め、プロットします。 プロットされるのは下図のような図です。

failing_to_reject.png

傾向として、データの数が多くなって x との相関が強いほど、帰無仮説は棄却されやすくなるのが分かります。

せっかくGPUが使えるようになったので、その効果ができるだけ出るようにしてみたいものだ。
ということで、とりあえずミニバッチの数をデフォルトの100から1000に増やしてみた。

Chainer$ python train_mnist.py -g 0 -b 1000
GPU: 0
# unit: 1000
# Minibatch-size: 1000
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           0.398479    0.165815              0.8895         0.9507                    0.972913      
2           0.130305    0.103915              0.9622         0.9661                    1.82951       
3           0.0807383   0.0763431             0.976533       0.9764                    2.68751       
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・   
18          0.000367914  0.0654017             1              0.9838                    15.6677       
19          0.000316801  0.0666274             1              0.984                     16.5294       
20          0.000283742  0.0675182             1              0.9839                    17.3905       
Chainer$  
実行時間は、17.39秒であり、b=100のとき(49秒)と比べて約2.8倍高速になった。 まとめて処理する単位を増やすと、全体の計算量は変わらないはずなのだが、かなり効率的になるようだ。
学習データ自体を使ったテスト(main/accuracy)は正解率100%になっているが、テストデータでは98.4%程度であり、デフォルトの場合との違いはほとんどない。

バッチサイズの変更で、GPUのメモリ使用量は、こんな感じで増加していった。
-b 100 (デフォルト)
|    0     16937    C   python                                         209MiB |

-b 1000
|    0     16917    C   python                                         317MiB |

-b 2000
|    0     16998    C   python                                         339MiB |

-b 5000
|    0     17031    C   python                                         493MiB |

-b 10000
|    0     17343    C   python                                         639MiB |
ミニバッチのサイズを大きくするとより高速になっていくと言いたいところだが、それは100を1000にした時には効果的だったのだが、10000にしても、速度はほとんど上昇しない。
Chainer$ python train_mnist.py -g 0 -b 10000
GPU: 0
# unit: 1000
# Minibatch-size: 10000
# epoch: 20

epoch       main/loss   validation/main/loss  main/accuracy  validation/main/accuracy  elapsed_time
1           1.37287     0.512147              0.660883       0.8614                    0.905861      
2           0.437023    0.361621              0.87135        0.901                     1.69545       
3           0.335709    0.294568              0.905283       0.9191                    2.50435       
・・・・・・・・・・・・・・・      中  略      ・・・・・・・・・・・・・・・   
18          0.0516563   0.0811957             0.986567       0.9746                    14.6888       
19          0.0475925   0.0765162             0.987367       0.9759                    15.4938       
20          0.0428052   0.074131              0.988917       0.9778                    16.29         
Chainer$ 
それよりも、accuracyが少し落ちてきたのがわかる。
ミニバッチ毎に学習が分けて行われるので、どうやら適度なサイズで行わないと、accuracyが落ちるようだ。
速度と精度のバランスを考えてミニバッチのサイズは決定しないといけないらしい。

ところで、MNISTデータセットは、学習用データが60000、テストデータが10000ある。
学習用データの数を変更してみようと思うが、どうやらコマンドパラメータからはできないらしい。
さて、どうしよう。
GPU(GEFORCE)がやってきて、CUDAが動くことが確認できた。
次にやることは、当然、ChainerからGPUを呼び出して高速に処理できるかどうか確認だ。

でも、どうすれば良いだろうか。
ネットをググるといろいろ見つかる。

Chainer, Cupy入門
Using GPU(s) in Chainer

でも、ちゃんと理解して読むのは面倒だ。
もっと横着して、とりあえず高速化できたことを実感したいだけなのだ。

と、思いながら眺めていると、sampleプログラムが用意されているらしい。
ということで、Example: Multi-layer Perceptron on MNIST にたどり着いた。
プログラムの説明がされているが、ここは無視して先へ行こう。

プログラムは,https://github.com/pfnet/chainer/tree/master/examples/mnist
にあり、その中の train_mnist.py を使おう。

落としてきたのを早速使おうと思ったが、引数などがわからない。
ということで、こうしてみた。

Chainer$ python train_mnist.py -h
usage: train_mnist.py [-h] [--batchsize BATCHSIZE] [--epoch EPOCH]
                      [--frequency FREQUENCY] [--gpu GPU] [--out OUT]
                      [--resume RESUME] [--unit UNIT]

Chainer example: MNIST

optional arguments:
  -h, --help            show this help message and exit
  --batchsize BATCHSIZE, -b BATCHSIZE
                        Number of images in each mini-batch
  --epoch EPOCH, -e EPOCH
                        Number of sweeps over the dataset to train
  --frequency FREQUENCY, -f FREQUENCY
                        Frequency of taking a snapshot
  --gpu GPU, -g GPU     GPU ID (negative value indicates CPU)
  --out OUT, -o OUT     Directory to output the result
  --resume RESUME, -r RESUME
                        Resume the training from snapshot
  --unit UNIT, -u UNIT  Number of units
Chainer$ 
このhelpから、GPUを使ったり、使わなかったりできるようである。
また、MNISTのデータセットによる学習の方法を色々変更できるようである。

ということで、とりあえず GPUを使ってみる。

このアーカイブについて

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

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

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