TIM Labs

2017年2月アーカイブ

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

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

createHistogram <- function(n) {
  lambdas <- vector(length = 1000)
  for (i in 1:1000) {
    y <- rpois(n, lambda = 4)
    data <- data.frame(y)

    fit <- glm(formula = y ~ 1, family = poisson, data = data)
    lambdas[i] <- exp(fit$coefficients[1])
  }
  title <- paste("Number of data =", n)
  hist(lambdas[lambdas <= 10], breaks = seq(0, 10, by = 0.01), main = title, xlab = "Estimeted lambda")
}

for (i in c(1, 5, 10, 100, 200)) {
  createHistogram(i)
}

実行すると平均 4 のポアソン分布から n 個のサンプルを取って、ポアソン回帰をするということを1000回繰り返して、推定結果のヒストグラムを描きます。 リンク関数は定数です。 n10 の時は下図のようになります。

estimation_number_of_data10.png

n100 の時は下図のようになります。

estimation_number_of_data100.png

推定結果が真の値である 4 の周囲に分布しています。 データの数が多くなれば分布のピークは 4 に集中していきますが、データが 100 個あっても、やはりずれる時はずれるのが分かります。

前回セーブしたものをロードして、ちゃんと動くか、脳の中身はちゃんと蘇るか確認してみよう。

saveの時と同じで、最初のおまじないのimportにつづいて、DigitsChainもimportしておく。
続いて、学習済みの 'digitslearnt.pkl' をモデルに読み込む。
データも 'digitstestdata.pkl' から読み込む。

このあたりは簡単過ぎて、説明することもないだろう。

from digitschain import DigitsChain

# 学習結果とTestデータをファイルから読み込む

import pickle
of = open('digitslearnt.pkl','rb')
model = pickle.load(of)
of.close()
of = open('digitstestdata.pkl','rb')
xtest,yans = pickle.load(of)
of.close()
ここまでで、前回セーブした時点の状態に戻ったはずである。

なので、ここから先はChainer:学習進行状況を確認しようで説明したのと同じコードで動作確認ができる。

# Testデータで学習成果を確認する
def check(data,answer,model):

    xt = Variable(data, volatile='on')	# 学習不要
    yy = model.fwd(xt)

    ans = yy.data
    nrow, ncol = ans.shape
    ok = sum([np.argmax(ans[i,:])==answer[i] for i in range(nrow)])
    return (ok * 1.0)/nrow	# 正解率
上記の内容を "digits0l.py" にまとめて書いて走らせると、こんな結果が得られた。
Chainer$ python digits0l.py
558 / 599  =  0.931552587646
Chainer$

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

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

legends <- c("ideal", "regression")
ltys <- c("dashed", "solid")
for (i in 1:10) {
  y <- rpois(5, lambda = 4)
  data <- data.frame(y)

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

  title <- paste("Poisson regression, estimated lambda = ", l)
  plot(1:10, dpois(1:10, lambda = 4), ylim = c(0.0, 0.25), type = "l", lty = "dashed", main = title, xlab = "", ylab = "")
  lines(1:10, dpois(1:10, lambda = l), type = "l")
  legend("topright", legend = legends, lty = ltys)
}

平均 4 のポアソン分布から5つのサンプルをサンプリングして、そのデータに対してポアソン回帰をしています。 リンク関数は定数です。 正確に推定されるなら、推定されるパラメータは 4 に近い値になるはずです。 実際に実行してみると、5個しかサンプルが無いわりには結構正確に推定されるのが分かります。 しかし、だいたい10回も推定すると、1回くらいは目に見えておかしな推定をするのが分かります。 例えば下図のような推定をします。

poisson_regresion_estimated4.png

真の値は 4 なのに、推定結果は 5.4 です。 このように、与えられたデータに対して最尤推定するしかない以上、推定は推定でしかなくどうしても真の値からずれてしまいます。

今回は、学習して、その結果をファイルにセーブする。

さて、学習結果なのだが、実際には model というDigitsChainクラスのオブジェクトの中に存在するので、このmodelをファイルにセーブすれば良いはずだ。
オブジェクトのセーブ、それもあれこれいっぱいデータが入っているはずのオブジェクトはどうやってセーブすれば良いだろうか?
Pythonには、オブジェクトをsave/loadするときに、オブジェクトを直列化してsaveするのが流儀であり、そのための仕組みが用意されている。

今回は、pickleを利用する。要するに、漬物(=pickle)にしてしまおうという訳だ。
詳しい説明はここではしないので、リンクを先を見るなどして自分で学習しよう。

import pickle
of = open('digitslearnt.pkl','wb')
pickle.dump(model,of)
of.close()
セーブは、たったこれだけで可能になる。 直列化してセーブしたデータの拡張子は".pkl"とする。
ファイル名とモード(バイナリの書き込みモード)でオープンし、そのファイルオブジェクトをofとしている。
次に、pickle.dumpで、セーブしたいオブジェクトと、ファイルオブジェクトを指定することで、ファイルに書き出してくれる。
最後にcloseする。

これだけだ。簡単過ぎる。

次回、別のプログラムで、この'digitslearnt.pkl'をロードしてテストしてみよう。
なので、テストデータもセーブしておこう。

of = open('digitstestdata.pkl','wb')
pickle.dump([xtest,yans],of)
of.close()
こんどは、2つのオブジェクトをセーブしないといけないので、リストにすることで、ひとつのオブジェクトに見えるようにした。 これで ".pkl" ファイルが2つできた。
Chainer$ ls -l *.pkl
-rw-rw-r-- 1 fuji fuji  20440  2月 15 15:56 digitslearnt.pkl
-rw-rw-r-- 1 fuji fuji 158385  2月 15 15:56 digitstestdata.pkl
今回は、ここまで。 次回、これらを読み込んで、テスト可能かどうか確かめよう。

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

この章ではポアソン回帰について説明がされています。 ポアソン回帰によって得られるのは、説明変数と応答変数の分布の関係です。 変数が変化すると別の変数の分布が変化するという、結構複雑な関係が得られるわけです。 分布そのものにも興味はありますが、やはり一番興味があるのは、分布の中でどこが一番確率が高くなるかです。 そこで、本の中で使われているサンプルデータを使ってポアソン回帰をやってみて、得られた分布の中の一番確率が高い場所をプロットしてみました。 コードはRで書きました。

d <- read.csv("data3a.csv")

fit <- glm(y ~ x, data = d, family = poisson)

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

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

fitted <- function(x, y) {
  dpois(y, lambda = linkFunction(x))
}

cat("p(y | exp(a + bx)) where p is Poisson distribution\n")
cat("a =", b1, "b =", b2, "\n\n")

y <- 1:20
max <- numeric(length = 20)
for (x in 1:20) {
  at <- fitted(x, y)
  cat("when x = ", x, ", y = ", which.max(at), " is most probable, probability is ", max(at), "\n")
  max[x] <- which.max(at)
}
title <- paste0("p(y|exp(", b1, " + ", b2, "x))")
xl <- "x"
yl <- "the most probable y"
plot(1:20, max, type = "b", main = title, xlab = xl, ylab = yl)

実行すると以下のようなグラフが得られます。

most_probable.png

サンプルデータの説明変数 x と応答変数 y には正の相関があるので、 x の増加に従って y の分布は大きな値を取るようになります。 グラフを見ると一番高い確率で得られる y の値が x の増加に従って増えているのが分かります。 正の相関の一面が見られますね。

リンク関数は指数関数なので、もっときれいに指数関数的増加が見られないかとプロットする x の範囲を広げてみました。

most_probable2.png

細かい増加は分かりにくくなりましたが、全体的な指数関数的増加がよく分かるようになりました。

Deep Learningを色々いじってきて、なんとなく学習出来ているのは分かったのだが、学習には結構時間が掛かる。そのため、一度学習した結果をファイルに出力しておき、必要になったときにファイルを読み込むことで学習することなく賢くして、使いたいものである。

ということで、学習結果(脳)のファイルへのsave/loadについて検討してみよう。

そのために、いままで利用してきたプログラムを使おう。

いままでのプログラムは、学習部分も、学習結果の判定のためのテスト部分も1つのファイルになっていたので、まずは2つの部分に分けようと思う。

このとき、学習結果が入っているオブジェクト(脳)は同じでなければならず、プログラムではそれを表すクラスが使われていた。
手書きデータの場合には、次のクラスがそれだ。

class DigitsChain(Chain):
    def __init__(self):
        super(DigitsChain, self).__init__(	
            l1=L.Linear(64,32),		# 1-2層
            l2=L.Linear(32,10),		# 2-3層
        )
        
    def __call__(self,x,y):
        return F.mean_squared_error(self.fwd(x), y)

    def fwd(self,x):
         h1 = F.sigmoid(self.l1(x))
         h2 = self.l2(h1)
         return h2
今までのファイルは、「学習&セーブ」と「ロード&テスト」の2つのファイルに別れる。
このとき、上記クラスは共通だから、これを"digitschain.py"という名前のファイルにする。

Pythonのモジュールは異常に簡単だ。 上のように、何でも良いから、.pyがついたファイルにしてしまえば、それだけでモジュールになってしまう。 他の言語のように、モジュールにするにはあれこれ宣言しないとダメとか全然ない。 そして、そのモジュールをインポートすれば使えるようになる。 記述フォーマットはいくつか種類があるが、モジュール名(=ファイル名)と、インポートする物を記述するのが次である。
from digitschain import DigitsChain
ということで、実際にimportしてみた。
>>> from digitschain import DigitsChain
Traceback (most recent call last):
  File "", line 1, in 
  File "/home/fuji/Study/Python/Chainer/digitschain.py", line 11, in 
    class DigitsChain(Chain):
NameError: name 'Chain' is not defined
>>> 
クラスだけのファイルはダメだった。
Chainというのは親クラスなので、
from chainer import Chain
の1行だけを加えたら、一応エラーは出なくなった。
L.やF.が使われているのにエラーにならないということは、メソッドを定義だけでは中身は検査されず、エラーにならないようだ。
でも、一応、おまじないを最初に全部加えておいた。
 
"digitschain.py"
# Define model

import numpy as np
import chainer
from chainer import cuda, Function, gradient_check, Variable 
from chainer import optimizers, serializers, utils
from chainer import Link, Chain, ChainList
import chainer.functions as F
import chainer.links as L

class DigitsChain(Chain):
    def __init__(self):
        super(DigitsChain, self).__init__(	
            l1=L.Linear(64,32),		# 1-2層
            l2=L.Linear(32,10),		# 2-3層
        )
        
    def __call__(self,x,y):
        return F.mean_squared_error(self.fwd(x), y)

    def fwd(self,x):
         h1 = F.sigmoid(self.l1(x))
         h2 = self.l2(h1)
         return h2
さて、次回は、これをimport して、学習内容をファイルに書き出すことをやってみよう。
前回、学習の進行状況をプリント出力できるようにしたのだが、数字の羅列を見るのは面倒なので、グラフとして見えるようにしよう。 そのために、学習ループの中で、リストresultに、[周回数,学習正解率,テスト正解率]を毎回追加するようにした。
# 学習ループ
result = []

for i in range(100000):
    x = Variable(xtrain)
    y = Variable(ytrain)
    model.zerograds()
    loss = model(x,y)		# lossを求める (forward)
    loss.backward()		# 微分(backward)
    optimizer.update()		# 調整
    if i % 100 == 0:
        lck = check(xtrain,trans,model)
        tck = check(xtest,yans,model)
        result.append([i,lck,tck])
        print("%6d   %6.4f   %6.4f" % (i,lck,tck))

lck = check(xtrain,trans,model)
tck = check(xtest,yans,model)
result.append([i,lck,tck])
print("%6d   %6.4f   %6.4f" % (i,lck,tck))
以上で集めたデータを、matplotlibで表示しやすいように分解し、表示するようにした。 グラフの内容をそのままpngのファイルとして出力もした。
# グラフ
import matplotlib.pyplot as plt
ilist = [d[0] for d in result]
lck   = [d[1] for d in result]
tck   = [d[2] for d in result]
plt.plot(ilist,lck,color="blue")
plt.plot(ilist,tck,color="red")
plt.savefig("digits0graph.png")
plt.show()

そして出力された画像がこれである。
digits0graph.png青が学習データでの成功率で、赤がテストデータでの成功率である。
最初、一気に .90位まで上昇し、その後、上昇ペースがどんどんゆっくりになっている。
青と赤の差が乖離誤差と呼ばれ、この差が激しくなると、過学習に陥っているという。
過学習とは、模擬試験には強いのだが本番はダメとか、練習の王者で本番ではメダルが取れないという感じとでも言えばよいだろうか。

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

この章にはポアソン回帰の例が載っています。 本の中で使われているサンプルデータがここで公開されているので、解説されている通りにポアソン回帰をやってみました。 コードはRで書きました。

d <- read.csv("data3a.csv")
cat("summary of input data\n")
summary(d)

fit <- glm(y ~ x, data = d, family = poisson)

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

cat("\nlink function\n")
cat("exp(", b1, " + ", b2, "x)\n")

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

ps <- 0:20
xl <- "y"
yl <- "probability"
legends <- c("x = 0","x = 5","x = 10", "x = 15", "x = 20")
pchs <- c(0, 5, 10, 15, 20)
title <- paste("lambda = exp(", b1, "+", b2, "x)")
plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 0.25), main = title, xlab = xl, ylab = yl)
for (x in seq(0, 20, by = 5)) {
  d <- function(y) {
    dpois(y, lambda = linkFunction(x))
  }
  lines(ps, d(ps), type = "l")
  points(ps, d(ps), pch = x)
}
legend("topright", legend = legends, pch = pchs)

実行するとポアソン回帰を行い、得られる分布をプロットしてくれます。 実際にポアソン回帰をやっているコードは以下の部分です。

fit <- glm(y ~ x, data = d, family = poisson)

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

b1b2 が最尤推定されたパラメータです。 推定された結果は b1 = 1.291721b2 = 0.07566191 でした。 これらのパラメータが得られればリンク関数 λ = exp(b1 + b2 x) が得られるので、説明変数 x に対して y がどのような平均 λ でポアソン分布するかが分かります。 その分布は下図のようになります。

glm_result.png

グラフを見ると x が大きくなるにしたがって y の分布が大きな値をとるように移動しているのが分かります。 これは、推定されたパラメータのうち、 b2 が正だからです。 サンプルデータの xy には正の相関があるのですね。

とりあえず、前回で学習ができているらしいことが確認できた。
しかし、それは、延々と学習した後で、テストデータとして温存しておいたデータにより確認しただけだ。
途中の学習状態、つまり、成功率がどのように変化していくのかを見たい!

ということで、最後にあった成功率を求める部分を、関数としてまとめてみた。

# Testデータで学習成果を確認する
def check(data,answer,model):

    xt = Variable(data, volatile='on')	# 学習不要
    yy = model.fwd(xt)

    ans = yy.data
    nrow, ncol = ans.shape
    ok = sum([np.argmax(ans[i,:])==answer[i] for i in range(nrow)])
    return (ok * 1.0)/nrow	# 正解率
データ、正解、モデルを与えると、このモデルを使ってデータを評価し正解との一致率を返す。

この関数(クラス内でないので、メソッドとは言わない)は、テストデータによる評価に使えるのだが、この形式にデータを合わせさえすれば、学習データを用いた場合の正解率も計算できる。

ということであるが、学習データ(教師データ)の正解の形式が違っているので、リストtransに教師データの答えを入れた。
# 学習データ(xtrain,ytrain)とテストデータ(xtest,yans)に分ける
index  = np.arange(N)
xtrain = X[index[index % 3 != 0],:]
ytrain = Y2[index[index % 3 != 0],:]
trans  = Y[index[index % 3 != 0]]	# 追加
xtest  = X[index[index % 3 == 0],:]
yans   = Y[index[index % 3 == 0]]
必要な追加は済んだので、最後の部分を書き換えた。
# 学習ループ

for i in range(10000):
    x = Variable(xtrain)
    y = Variable(ytrain)
    model.zerograds()
    loss = model(x,y)		# lossを求める (forward)
    loss.backward()		# 微分(backward)
    optimizer.update()		# 調整
    if i % 1000 == 0:
        print("%6d   %6.4f   %6.4f" % (i,check(xtrain,trans,model),check(xtest,yans,model)))

print("%6d   %6.4f   %6.4f" % (i,check(xtrain,trans,model),check(xtest,yans,model)))

forループの制御変数iを、forループを抜けてからも使っている。
C言語ならエラーになるところだが、Pythonではそんなことはないのだ。

これで、こんな感じに動いた。
Chainer$ python3 digits0ck.py
     0   0.0751   0.0634
  1000   0.5601   0.5309
  2000   0.6820   0.6477
  3000   0.7821   0.7379
  4000   0.8297   0.7863
  5000   0.8706   0.8280
  6000   0.8890   0.8715
  7000   0.9040   0.8765
  8000   0.9124   0.8865
  9000   0.9249   0.8948
  9999   0.9316   0.9032
Chainer$ 

数字は、「繰り返し回数、教師データでの正解率、テストデータでの正解率」である。

まだまだ成功率上昇中なので、もっともっとループを繰り返してみた。

226000   0.9958   0.9666
227000   0.9958   0.9666
228000   0.9958   0.9699
229000   0.9958   0.9716

学習データでのテストは99.5%を超えて非常に高いが、テストデータでは 97%を超えるのがやっとであった。つまり、過学習状態になっているのが分かる。
これ以上学習させても、このモデルでは成績は微小な変動を繰り返すのみであった。

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

この章では主にポアソン回帰について解説されています。 ポアソン回帰というのは、ポアソン分布の平均がリンク関数という関数を通してデータと関連するモデルを使った回帰です。 モデルが複雑になったので、この回帰で得られる物が何なのかをグラフに描くコードを用意しました。 コードはRで書きました。

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

distribution <- function(x, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dpois(y, lambda = l(x))
  }
}

ps <- 0:10
xl <- "y"
yl <- "probability"
legends <- c("x = 0","x = 5","x = 10", "x = 15", "x = 20")
pchs <- c(0, 5, 10, 15, 20)
for (b1 in seq(-0.1, 0.1, by  =  0.05)) {
  for (b2 in seq(-0.1, 0.1, by = 0.05)) {
    l <- linkFunction(b1, b2)
    title <- paste("lambda = exp(", b1, "+", b2, "x)")
    plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
    for (x in seq(0, 20, by = 5)) {
      d <- distribution(x, l)
      lines(ps, d(ps), type = "l")
      points(ps, d(ps), pch = x)
    }
    legend("topright", legend = legends, pch = pchs)
  }
}

この例では、説明変数xに対してyが分布する場合について考えています。 リンク関数を作るコードは以下のようになります。 パラメータb1b2を与えると、そのパラメータを使ったリンク関数を返しています。

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

リンク関数は引数xを取って、パラメータと足したり掛けたりした後、累乗を計算する関数です。 ポアソン回帰ではリンク関数の戻り値を、ポアソン分布のパラメータとして使います。

リンク関数があれば、あるxに対するyの分布が得られます。 つまり、yを引数に取りyが起きる確率を返す関数が得られます。 以下の関数により得られる関数でyの起きる確率が計算されます。

distribution <- function(x, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dpois(y, lambda = l(x))
  }
}

これがポアソン回帰のモデルです。 後はリンク関数で使われるパラメータを最尤推定で求めれば良いわけです。

今はとりあえずモデルの具体例が見たいだけなので、パラメータは自分で与えて、様々な値についてグラフをプロットしてみました。 コードを実行すると、いろいろな場合についてグラフがプロットされます。 具体例として、観ていて面白いのはb1 = -0.1b2 = 0.1の場合です。 以下に図を貼ります。

glm_b1_-0.1b2_0.1.png

リンク関数はexp(-0.1 + 0.1x)になります。 xがある数だった場合、yはこの値を平均に持つポアソン分布で分布するということです。 例えばx = 1の場合は1です。 x = 11の場合はe = 2.718...です。

b2が正なのでxが大きくなるほどリンク関数の指数の肩は大きな数になります。 つまり、大きなxではyは大きな平均を持つポアソン分布をします。 図を見ると、確かにxが大きくなるにしたがってyの分布は大きな値を取るようになっています。

このように、ポアソン回帰で得られるのはxの変化によってyの分布がどのように変化するかという傾向です。 ただ変化すると言うだけでなく、どれくらいの分布を伴って変化するかも分かります。

前回で、データを読み込み、DLのモデルを DigitsChainクラスとして定義した。
今回は、このDigitsChainクラスよりモデルを作り、読み込んでいたデータを使って、学習とテストをやってみる。

実際のプログラムは以下のようになる。
# 学習モデルの初期化

model = DigitsChain()           # モデルの生成
optimizer = optimizers.SGD()
optimizer.setup(model)

# 学習ループ

for i in range(10000):
    x = Variable(xtrain)
    y = Variable(ytrain)
    model.zerograds()
    loss = model(x,y)		# lossを求める (forward)
    loss.backward()		# 微分(backward)
    optimizer.update()		# 調整

# Testデータで学習成果を確認する

xt = Variable(xtest, volatile='on')	# 学習不要
yy = model.fwd(xt)

ans = yy.data
nrow, ncol = ans.shape
ok = sum([np.argmax(ans[i,:])==yans[i] for i in range(nrow)])
print( ok, "/", nrow, " = ", (ok * 1.0)/nrow )		# 正解率
これって、Irisのときと同一に見えるはずだ。 たった1行違うだけ。
model = DigitsChain()           # モデルの生成
つまり、粗い手書き文字用に作ったDigitsChainクラスのオブジェクトを作って、前と同じ変数modelに与えている。

今回の変更は1行、それ以前の変更も数行であり、全体で10行も変更していない。
というか、変更する必要がありそうなところが見つからない。

さて、とりあえず動かしてから考えよう。

Chainer$ python3 digits0.py
546 / 599  =  0.911519198664
Chainer$ 

テスト結果は91.15%の正解率で、ちゃんと動いているようだ。

IrisからDigitsへの変更は、同じクラス分けタイプのデータを使ってクラス分けしたので、僅かな変更で済んだ。
4種のデータからなるIrisと、8x8の画像データからなるDigitsが、ほぼ同じプログラムで動いてしまうのだ。
この適応性の高さは恐ろしいほどだ。

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

この章では主にポアソン分布と、ポアソン分布のパラメータの最尤推定について解説されています。 ポアソン分布とは、下図のような分布です。

pois5.png

ポアソン分布はパラメータを一つ持っていて、それが決まると分布が決まります。

最尤推定とは、あるデータが得られている時、そのデータが生成される確率が最大になるように、分布のパラメータを選ぶ推定方法です。 データは何でもいいですが、現実の世界で測定された数値です。 その数値をモデルを使って説明しようというのですから、実際に得られているデータに最もよく適合するモデルが欲しいと思うのは当然です。 最もよく適合するというのは、確率の言葉で言えば最大の確率を与えるということなので、このようにパラメータを推定するのが合理的ですね。

さて、問題が単純な場合には、解析的に確率を最大化するパラメータを求めることができるのですが、複雑な問題では普通は数値計算で求めます。 その様子を見やすくするデモを用意しました。 コードはRで書きました。

n <- runif(1, 1, 10)
data <- rpois(20, lambda=n[1])

cat("input data taken from Poisson distribution (lambda = ", n[1], ")\n\n")
data
summary(data)
cat("\n")

cat("likelihood of given data at\n")
res <- vector(length = 10)
for (i in 1:10) {
  prob <- function(x) {
    dpois(x, lambda=i)
  }
  res[i] <- sum(log(prob(data)))
  cat("lambda = ", i, " ", res[i], "\n")
}

cat("\nmaximum likelihood = ", max(res), "when lambda = ", which.max(res), "\n")

このコードが何をやっているか詳しく解説します。 まず最尤推定に使うデータを用意しています。 以下のコードで 1 から 10 までの乱数を一つ用意して、その乱数を平均に持つポアソン分布からサンプルを 20 用意しています。

n <- runif(1, 1, 10)
data <- rpois(20, lambda=n[1])

後は 1 から 10 までの整数 i について、このデータが平均 i を持つポアソン分布から得られる確率を計算しています。 確率を計算しているのは以下のコードです。

prob <- function(x) {
  dpois(x, lambda=i)
}
res[i] <- sum(log(prob(data)))

dpois(x, lambda=i)λ=i のポアソン分布から x が得られる確率です。 sum(log(prob(data)))data に含まれる x の全てについて、 dpois(x, lambda=i) を計算し対数を取って和を計算しています。確率の積をそのまま計算すると値が小さくなりすぎて計算できないので対数を取って和を計算しました。ちなみにこの確率には尤度という名前がついています。

後は結果をコンソールに出力しています。 実行結果は次のようになります。

input data taken from Poisson distribution (lambda =  2.923534 )

 [1] 5 3 4 6 5 4 0 2 6 3 6 0 3 4 7 2 1 1 2 2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     2.0     3.0     3.3     5.0     7.0 

likelihood of given data at
lambda =  1   -75.51993 
lambda =  2   -49.77221 
lambda =  3   -43.01152 
lambda =  4   -44.0245 
lambda =  5   -49.29702 
lambda =  6   -57.2638 
lambda =  7   -67.08986 
lambda =  8   -78.27679 
lambda =  9   -90.5031 
lambda =  10   -103.5493 

maximum likelihood =  -43.01152 when lambda =  3

データを生成するのに使われたパラメータは 2.923534 です。 なので、尤度はそれに近い 3 の時に最大になっています。 何度か実行してみると、データを生成するのに使ったポアソン分布の平均に近い時に尤度が一番大きくなっていることが分かると思います。 データの本当の分布は最初に引いた乱数を平均のポアソン分布なのだから、その乱数に近い推定を行った時に、最も当てはまりがよくなるのは順当な結果です。 このように最尤推定を行えば、統計処理したいデータに対し、順当な統計分布を見つけることができるのですね。


今回は、読み込んで、Irisのときと同様に、画像のデータをX、教師データをYに読み込む。
前回は、出力が3つだったが、今回は数字の0〜9までで10個に変わるのを考えて、以下のように書き換えた。
 
# digitsデータの読み込み
from sklearn import datasets
digits = datasets.load_digits()
X = digits.data.astype(np.float32)
Y = digits.target.astype(np.int)
N = Y.size
Y2 = np.zeros(10 * N).reshape(N,10).astype(np.float32)
for i in range(N):
    Y2[i,Y[i]] = 1.0
このデータを学習データとテストデータに分割する訳だが、データ数が150から1797に増加したので、2/3を学習データに、1/3をテストデータにしてみた。

# 学習データ(xtrain,ytrain)とテストデータ(xtest,yans)に分ける
index  = np.arange(N)
xtrain = X[index[index % 3 != 0],:]
ytrain = Y2[index[index % 3 != 0],:]
xtest  = X[index[index % 3 == 0],:]
yans   = Y[index[index % 3 == 0]]
次に、一番肝心なモデルについて変更する。
Irisのときのモデルのノード構成は、4-6-3 となっていた。中間層が6ノードである。

digitsのとき、入力ノードは8x8の画像なので、64になる。
出力ノードは、0〜9の各値に対してノードがあるので、10になる。
中間層のノードであるが、64と10から適当な値を考えるべきである。
しかし、ここではとりあえず、64の半分の32を与えてみる。

中間層のノード数が変わっても、プログラム中で変更するのは、たった2個所である。
中間層のノード数に名前を付けてしまえば、一箇所にすることも可能だ。

# Define model
class DigitsChain(Chain):
    def __init__(self):
        super(DigitsChain, self).__init__(	
            l1=L.Linear(64,32),		# 1-2層
            l2=L.Linear(32,10),		# 2-3層
        )
        
    def __call__(self,x,y):
        return F.mean_squared_error(self.fwd(x), y)

    def fwd(self,x):
         h1 = F.sigmoid(self.l1(x))
         h2 = self.l2(h1)
         return h2
これでほぼ準備は出来たはずなので、走らせるのは次回にしよう。
元々のプログラムはIrisのデータ用になっていて、そのデータに対して3層のニューラルネットワーク(NN)でDeep Learning をし、学習結果をテストするようになっていた。
このとき、次のような図をイメージするのだ。
この図は、
R for Deep Learning (I): Build Fully Connected Neural Network from Scratch
からのものだ。
あ、これはPythonではなく、RでNNをやる説明のようだ。まあいいか、言語なんて大差ない。

図にはデータの流れが示されていないが、データは左から入り、右に向いて流れていく。
3つの層を構成するノードはタテに並んだ丸で示され、左から、入力層、中間層、出力層になっている。

iris_network.png
この図を、次の『Chainerによる実践深層学習』のIrisのプログラムのうち、上の図に対応した部分だけを示す。

class IrisChain(Chain):
    def __init__(self):
        super(IrisChain, self).__init__(
            l1=L.Linear(4,6),
            l2=L.Linear(6,3),
        )
        
    def __call__(self,x,y):
        return F.mean_squared_error(self.fwd(x), y)

    def fwd(self,x):
         h1 = F.sigmoid(self.l1(x))
         h2 = self.l2(h1)
         return h2

このうち、今注目すべきは、以下の部分だけである。
            l1=L.Linear(4,6),              # 入力層→中間層
            l2=L.Linear(6,3),              # 中間層→出力層

数字が、上の図ときちんと対応しているのが分かるだろう。

入力層の4は、データが4種の値から構成されているからである。

中間層の6は、まあ適当に決めて良いのだが、詳しい説明は面倒なので省略する。

出力層の3は、Irisでは3種類の種の分類(クラス分け)をしたいので、3になっている。


ここまで理解したら、前回説明した8x8のめちゃくちゃ粗い手書き数字データの分類に挑戦してみよう。
ということで、次回につづく。
DLの最初に使われるデータ例が、ほとんどの場合、iris(アヤメ)のデータだ。
3種類のアヤメについて、花びらの長さと幅、がくの長さと幅のデータが、それぞれのアヤメの種類に対して50組ずつデータがあり、全体で150組のデータがある。
これから学習データとテストデータに分けて、学習後、テストデータによりアヤメの種類を正しく分類できるかようになっているか調べるのであった。

でも、これには飽きたので、他のデータで試したい。
そのためには、まず、アヤメのデータがどこに存在し、どのように利用されているかを調べないとダメだ。

from sklearn import datasets
iris = datasets.load_iris()
X = iris.data.astype(np.float32)
Y = iris.target
N = Y.size
Y2 = np.zeros(3 * N).reshape(N,3).astype(np.float32)
for i in range(N):
    Y2[i,Y[i]] = 1.0

この中で肝心なのは最初に1行で、datasetsをインポートし、その中からirisのデータを読み込んでいることが分かる。
そして、scikit-learnのサイトにたどりついた。


scikit-learn_top.png
この中のドキュメントに、それらしい情報があった。

5.2. Toy datasets

scikit-learn comes with a few small standard datasets that do not require to download any file from some external website.

load_boston([return_X_y]) Load and return the boston house-prices dataset (regression).
load_iris([return_X_y]) Load and return the iris dataset (classification).
load_diabetes([return_X_y]) Load and return the diabetes dataset (regression).
load_digits([n_class, return_X_y]) Load and return the digits dataset (classification).
load_linnerud([return_X_y]) Load and return the linnerud dataset (multivariate regression).

These datasets are useful to quickly illustrate the behavior of the various algorithms implemented in the scikit. They are however often too small to be representative of real world machine learning tasks.

この中のload_irisが呼び出されていたのだ。

これを他のものに変えて、元のプログラムをちょっとだけ変更して、動くかどうか確かめよう。

上の表のloadメソッドをクリックすると、それぞれの説明が現れる。説明から、iris to digits がclassificationのためのデータらしいので、とりあえず同じタイプのdigitsを使うことにしよう。

とりあえず、load_digitsに載っているサンプルでデータのロード表示をしてみよう

from sklearn.datasets import load_digits
>>> digits = load_digits()
>>> digits.data.shape
(1797, 64)
>>> import matplotlib.pyplot as plt
>>> plt.gray()
>>> plt.matshow(digits.images[0])
<matplotlib.image.axesimage object="" at="" 0x7f5c88085a90>
>>> plt.show()
データが1797個あり、1つのデータは64個のデータからなる。 matplotlibで最初のデータを表示している。 1つのデータは64個の値からなるのだが、実際は8x8の画像データであることが分かる。 実際に最初のデータの内容を以下のようにして数値で見ると、画像との対応がよく分かる。 scikit-learn_zero.png
>>> X = digits.data.astype(np.float32)
>>> X[0].reshape(8,8)
array([[  0.,   0.,   5.,  13.,   9.,   1.,   0.,   0.],
       [  0.,   0.,  13.,  15.,  10.,  15.,   5.,   0.],
       [  0.,   3.,  15.,   2.,   0.,  11.,   8.,   0.],
       [  0.,   4.,  12.,   0.,   0.,   8.,   8.,   0.],
       [  0.,   5.,   8.,   0.,   0.,   9.,   8.,   0.],
       [  0.,   4.,  11.,   0.,   1.,  12.,   7.,   0.],
       [  0.,   2.,  14.,   5.,  10.,  12.,   0.,   0.],
       [  0.,   0.,   6.,  13.,  10.,   0.,   0.,   0.]], dtype=float32)
>
これから、digitsのデータは1797個の画像データであることが分かる。

とりあえず、データの取得が出来たので、今回はここまで。

このアーカイブについて

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

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

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

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