TIM Labs

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

最近のコメント