TIM Labs

2017年6月アーカイブ

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの四回目です。

この章ではMCMCやメトロポリス法について説明がされています。 その中で、メトロポリス法によるサンプリングと尤度が関係あることが解説されています。 なので、メトロポリス法によるサンプリングの分布と尤度を比べられるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

sampling <- function(n, init) {
  q <- init
  ql <- likelihood(q)
  samples <- vector(length = n)
  i <- 1
  while(i <= n) {
    p <- develop(q)   
    pl <- likelihood(p)
    r <- log(runif(1, 0, 1))
    if (ql < pl || r < pl - ql) {
      q <- p
      ql <- pl
    }
    samples[i] <- p
    i <- i + 1
  }
  samples
}

qs <- seq(0, 1, by = 0.01)
ls <- vector(length = length(qs))
for(i in 1:length(qs)) {
  ls[i] <- prod(dbinom(data, N, qs[i]))
}
plot(qs, ls, type = "l", xlab = "q", ylab = "Likelihood")

samples <- sampling(10000, 0.3)
hist(samples, breaks = seq(0, 1, by = 0.01), main = "10000 samples", xlab = "q")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 コードを実行すると次のような図がプロットされます。

metro_likelihood.png

この図は q を変化させた時の尤度のグラフです。 発生させたサンプルデータの偏り方によって多少ずれますが、だいたい 0.45 を中心とした山なりなグラフになります。

metro_hist.png

この図はサンプル数 10000 の場合のメトロポリス法によるサンプリングの結果のヒストグラムです。 確かに、さきほどプロットした尤度に比例した分布が得られていることが分かります。

ダメな統計学 (283x400).jpg

題名:ダメな統計学

悲惨なほど完全なる手引書

アレックス・ラインハート 著  西原 史暁 訳

A5、185項、本体2200円

2017年1月20日 発行

勁草書房

https://www.statisticsdonewrong.com/ (オリジナルサイト)

http://id.fnshr.info/2017/01/20/sdw-trans-publ/ (日本語サイト)


最近、世の中、統計学が流行っているようである。

確かに昔と違って、大量のデータもさくっと処理できるだけの計算機パワーがあり、様々な分析、予想、判断などに実際に利用できるようになった。

そして、ビッグデータ、人工知能などのためには、統計学は当然知っていなければいけない。


ということで、統計学を勉強と思ってではなく、世の中、どう考えてもオカシイ、誤った統計の利用、誤用、乱用、悪用が目に付いていたところでこの本が出版されたので読んでみた。


第一印象は、結構読みづらい本であった。

日本語訳があまりよくない。元々内容がねちねちしたところがあるのは当然なのだが、それゆえしっかりした、かつ読みやすい日本語にして欲しいところだが、その当たりはかなり改善して欲しいところだ。英語がすらすら読めるのならば、英語で読むべきかもしれない。

内容は、論文などでの統計の利用がいかに正しくないかが書かれている。
世界の主要な論文掲載誌に載る論文も、1/3程度は統計的におかしなことを書いており、結論が違っているとか。

とくに、医学、医療関連、つまり治療の有効性の判定などでの誤りの指摘が多い。
医療系の世界では、統計の誤解だけでなく、どうやら政治、ビジネス的な問題がからんで、世界中どこもドロドロなところがあると。

統計とくに検定についての言及が多いが、そのあたりを正しく把握しようと思うと、一通りの統計学の知識が必要になる。

いろいろ統計的に誤った実例が載っている。とくに、薬で、有効と認定されていたのが無効と判定されてしまったいきさつなどの紹介もある。

こういう研究は、直接的な成果に結びつかないので、関心を持たれることがないが、統計にきちんと基づいて論文の査読がされれば、誤った結論を出してしまうことが相当減少することは確実だ。

一部の論文誌では、査読に統計、検定の根拠を示すのを義務化するのも進み始めたみたいだ。

本書は、さくっと読んで、しっかり統計学を勉強するのが重要だ。


なお、本書の内容であるが、英語も、翻訳もネットでフリーで読むことができる。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの三回目です。

この章ではMCMCやメトロポリス法について説明がされています。 その中で、メトロポリス法によるサンプリングが解説されています。 なので、同じようにメトロポリス法でサンプリングしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

sampling <- function(n, init) {
  q <- init
  ql <- likelihood(q)
  samples <- vector(length = n)
  i <- 1
  while(i <= n) {
    p <- develop(q)   
    pl <- likelihood(p)
    r <- log(runif(1, 0, 1))
    if (ql < pl || r < pl - ql) {
      q <- p
      ql <- pl
    }
    samples[i] <- p
    i <- i + 1
  }
  samples
}

samples1 <- sampling(500, 0.3)
hist(samples1, breaks = seq(0, 1, by = 0.01), main = "500 samples", xlab = "q")

samples2 <- c(samples1, sampling(500, samples1[500]))
hist(samples2, breaks = seq(0, 1, by = 0.01), main = "1000 samples", xlab = "q")

samples3 <- c(samples2, sampling(1000, samples2[1000]))
hist(samples3, breaks = seq(0, 1, by = 0.01), main = "2000 samples", xlab = "q")

samples4 <- c(samples3, sampling(8000, samples3[2000]))
hist(samples4, breaks = seq(0, 1, by = 0.01), main = "10000 samples", xlab = "q")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 なので、このデータに対して q のサンプリングをすると、だいたい 0.45 を中心とした分布が得られるはずです。

コードではサンプル数が 5001000200010000 の時にヒストグラムをプロットしています。 コードを実行すると次のような図がプロットされます。

metro_500.png

この図はサンプル数 500 の場合です。 まだヒストグラムが凸凹しているのが分かります。

metro_10000.png

この図はサンプル数 10000 の場合です。 だいぶヒストグラムがなめらかになって、分布が収束してきているのが分かります。

PythonをAIなどで使おうと思ったら、かならず使うことになるのがNumPyである。 様々な計算処理をしようとすると、はやり配列を使って高速に計算する必要がでてくるので、配列をサポートしているNumPyは避けて通れない。
 NumPyにあるtutorilaなどを参考に勉強したり、ネット上にウヨウヨ存在する説明をみて勉強するのもあるのだが、もうちょっと面白い勉強方法、あるいは腕試しには、100 numpy exercises が良いかも。

題名どおり、100の小問で構成されている。
難易度が3段階になっていて、やさしい、基本的なものから、次第に難しい問題になっているようだ。
最初は1、2行で十分な簡単な問題ばかりだが、次第に長いプログラムを書かないといけなくなり、20行程度のプログラムになることもある。

問題と解答があり、どうしてそう書けば良いのかは自分で考えよう、というスタンスのようである。
もちろん、世の中には延々と解いて、解説ページを公開している人もいる。

解説ページを読むより、解答が理解できなかったら、ネットで色々探して、理解できるようになろう。
さらに、出題範囲以外のことを色々試すことをお薦めする。

あ、書き忘れていたことがある。
もちろん、全部英語だけれど、やさしいプログラミング英語なので、問題ないだろう。
どうせ、プログラムをいっぱいやるようになったら、英語しかない世界にすぐなるのだし、プログラミング英語に慣れる良い機会だ。

まず、最初の問題。こんなに簡単で、分からない人はいないはずだ。

1. Import the numpy package under the name np (★☆☆)

import numpy as np

そして最終問題は、こんな感じだ。

100. Compute bootstrapped 95% confidence intervals for the mean of a 1D array X (i.e., resample the elements of an array with replacement N times, compute the mean of each sample, and then compute percentiles over the means). (★★★)

# Author: Jessica B. Hamrick

X = np.random.randn(100) # random 1D array N = 1000 # number of bootstrap samples idx = np.random.randint(0, X.size, (N, X.size)) means = X[idx].mean(axis=1) confint = np.percentile(means, [2.5, 97.5]) print(confint)


こんな問題に対するプログラムをすらすら書けるようになったら、もちろんNumPyのお勉強は優秀な成績で卒業だ。
NumPyの勉強は、通過すべき最初のゲートみたいなものだ。どんどん先へ行こう。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの二回目です。

この章ではMCMCやメトロポリス法について説明がされています。 なので、実際にメトロポリス法を動かしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)
summary(data)
cat("Analytically, q =", sum(data) / (dataSize * N), "\n\n")

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

q <- round(0.3 + runif(1, 0, 0.3), 2)
ql <- likelihood(q)
i <- 0
cat("start q =", q, "\n")
while(i < 30) {
  p <- develop(q)   
  pl <- likelihood(p)
  r <- log(runif(1, 0, 1))
  if (ql < pl || r < pl - ql) {
    q <- p
    ql <- pl
  }
  cat("q =", q, "likelihood", ql, "\n")
  i <- i + 1
}
cat("last q =", q, "\n")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 このデータに対して、事象の発生確率 q が分かっていないものとして推定していきます。

メトロポリス法ではランダムに q を変化させた後、尤度が低くなる場合でも、変化の前後で尤度の比を計算して、その確率で q を更新します。 今の場合、計算の途中で使っているのは対数尤度なので、比ではなくて差になります。 コードで言うと以下の部分です。

  r <- log(runif(1, 0, 1))
  if (ql < pl || r < pl - ql) {
    q <- p
    ql <- pl
  }

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

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    2.0     3.0     4.0     3.7     4.0     6.0
Analytically, q = 0.4625

start q = 0.31
q = 0.31 likelihood -40.44572
q = 0.32 likelihood -39.35181
q = 0.32 likelihood -39.35181
q = 0.32 likelihood -39.35181
q = 0.33 likelihood -38.34881
q = 0.33 likelihood -38.34881
q = 0.34 likelihood -37.43294
q = 0.35 likelihood -36.60087
q = 0.36 likelihood -35.84958
q = 0.37 likelihood -35.17642
q = 0.38 likelihood -34.579
q = 0.39 likelihood -34.05522
q = 0.4 likelihood -33.60322
q = 0.39 likelihood -34.05522
q = 0.4 likelihood -33.60322
q = 0.41 likelihood -33.22138
q = 0.42 likelihood -32.90828
q = 0.42 likelihood -32.90828
q = 0.43 likelihood -32.66271
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.45 likelihood -32.37025
q = 0.46 likelihood -32.32184
q = 0.45 likelihood -32.37025
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
last q = 0.43

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。 メトロポリス法は確率的に尤度が低くなる変化もするので、q は一直線に最尤推定値に向かうのではなく、時々戻ったりするのが分かります。

Deep Learning していると、学習にとても時間がかかるのがネックになる。
時間だけでなく、メモリもいっぱい必要になる。
今はGPUなるものがあり、これを使うと計算処理が一気に早くなると言われている。
GPUの計算性能は、CPUの50倍以上とか言われているのだが、実際にはそこまで高速にDeep Learningできるようにはならない。
GPUを使う場合には、GPUにデータを送ったり、GPUからデータを受け取ったりなど、CPUとGPUの間でのデータのやり取りが発生し、この量が半端ではない。このため、実際の速度向上は10倍程度にしかならない。

でも、もっと高速にしたいと思ったらどうすればよいだろうか。
多数のGPUを同時に使えば、どんどん高速になるのではないだろうか。
分散して学習する場合に問題になるのは、別々のGPUで学習した結果をまとめて、それをまた全GPUに分配し直すという作業が必要である。これをいかに上手にするかで性能が殆んど決まってしまう。

これは人間にたとえると、多人数で勉強をするにの、それぞれ別の参考書を読んで学習し、読み終わったら全員で学習結果を寄せ集めることでより賢くし、学習結果をまた全員に分配することで、全員が全参考書を読んで賢くなったようにしてしまう。そして、またそれぞれが別の参考書で学習すれば、また賢くなっていく。これを繰り返せば、相当な速度で賢くなるはず。
一人でやるより、二人でやれば2倍近く賢くなり、4人、8人、16人、、、、と増やしていくと、賢くなる速度が倍々で速くなるはずだ。

でも、そんなにうまくいくだろうか。

ということで調べていたら、こんな記事にぶち当たった。
ChainerMN による分散深層学習の性能について

そして、この中には、以下のグラフが含まれていた。

dlsummit_01_throughput.png
GPUを128個使うと、100倍高速になったとある。
つまり、ほぼGPUの数に比例して高速になるという夢のような効率である。
本当にそこまで高速になるのか、どうしても疑いたくなってしまう。

さらに、他のフレームワークとのスケーラビリティの比較があった。

dlsummit_05_framework-throughput.png
ChainerはPythonなので、単純比較してしまうと高速ではないはずだ。
しかし、GPUをたくさん使った場合には、他のDeep Learning よりも遥かに高速になるらしい。
以上は、Chainerを作っているPrefered NetworksのPreferred Researchにあったものの引用である。
もっと詳しく知りたい場合は、ぜひオリジナルを読んでみよう。

このGPUのマルチノード対応パッケージは、既に公開されているので、GPUがゴロゴロしている人、GPUをふんだんに使える環境にある場合は、ぜひ実験してみよう。
詳細は、分散深層学習パッケージ ChainerMN 公開 で確認のこと。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの一回目です。

この章ではMCMCやメトロポリス法について説明がされています。 本の中ではMCMCの前に、まずはランダムにパラメータを更新していく手法が紹介されています。 なので、本の中で紹介されているような推定法を実際に動かしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)
summary(data)
cat("Analytically, q =", sum(data) / (dataSize * N), "\n\n")

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

q <- round(0.3 + runif(1, 0, 0.3), 2)
ql <- likelihood(q)
i <- 0
cat("start q =", q, "\n")
while(i < 30) {
  p <- develop(q)   
  pl <- likelihood(p)
  if (ql < pl) {
    q <- p
    ql <- pl
  }
  cat("q =", q, "likelihood", ql, "\n")
  i <- i + 1
}
cat("last q =", q, "\n")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 このデータに対して、事象の発生確率 q が分かっていないものとして推定していきます。

推定方法は q をランダムに変化させて、対数尤度を計算し、尤度が大きくなっていたら変化させた先の q に更新する、という方法です。 develop 関数の中で q0.01 足すか引くかして q を変化させています。 コードを実行すると次のような出力が得られます。

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      1       3       4       4       5       6
Analyticaly, q = 0.5

start q = 0.35
q = 0.36 likelihood -41.69397
q = 0.37 likelihood -40.76192
q = 0.38 likelihood -39.90848
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.4 likelihood -38.42821
q = 0.4 likelihood -38.42821
q = 0.41 likelihood -37.79737
q = 0.42 likelihood -37.23712
q = 0.42 likelihood -37.23712
q = 0.43 likelihood -36.74602
q = 0.44 likelihood -36.32282
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.46 likelihood -35.67609
q = 0.47 likelihood -35.45097
q = 0.48 likelihood -35.29055
q = 0.49 likelihood -35.19445
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
last q = 0.5

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。 q をランダムに変化させていくと、最終的にはその値に近づいていくのが分かります。

Chainerはどういう環境で動くであろうか。
Installation Guideを見ると、

  •     Ubuntu 14.04/16.04 LTS 64bit
  •     CentOS 7 64bit
とあるのだが、日本で特に多いWindowsが載っていない。
Windowsをサポートしていなくても、クラウド上のChainerを使って、手元のPCはWindowsというのはできるようだが、私はやったことがない。

それが、ついに、ChainerがMicrosoft Azure および Windows にも対応するようになるようだ。
詳細は、ASCII.jp の 「Chainer」をMicrosoft Azureに最適化、Windowsもサポート を参照のこと。

Microsoftサイトでの紹介:Preferred Networks とマイクロソフト、ディープラーニングソリューション分野で戦略的協業

ところで、Microsoft自体は、CNTK(Microsoft Cognitive Toolkit)を持っている。
CNTK(Microsoft Cognitive Toolkit)のホームページ(英語)があり、GitHubにてMITライセンスにてオープンソースとして公開されている。

マイクロソフトというと、まだ頭の固い連中は、オープンソースと対極にある会社と思っている者だらけだと思うが、ずいぶん前からOSSの活動を行っている。まだ、ビル・ゲイツが現役の最後の頃に日本にやってきて、「マイクロソフトはこんなにOSSをやっていて....」と説明したものの、そもそも相手がOSSなるものが何かも分からない連中で、全然マイクロソフトのOSS活動を伝えられなくて、ガックリとして帰国したことがあった。

マイクロソフトは、クラウドビジネスの拡大を目指しており、自社のソフトという縛りをつけていたらビジネスにならないので、積極的にOSS活動をしている。たとえば、 マイクロソフトによるオープン ソースへの取り組みが、ビジネス チャンスを拡大を見れば、その一端が分かるだろう。


2017年6月現在、Deep Learningのツール、フレームワークはいったいいくつあるのだろうか?
まだ、これからも増えそうな気がする。
また、各フレームワークは、機能、性能、開発環境などをどんどん向上させてきているので、どれが良いといっても、それは判断した瞬間だけしか有効でない。
GPUの性能向上も著しく、うまくGPUを活かせるのはどれかも重要な判断材料になる。
数年後に、Deep Learningのフレームワークの競争はどのようになっているか、今使われているもののうちどれが優勢になるかの判断は困難だ。
2020年頃には、消えてなくなる、弱体化しているのが多いだろうし、それはITの世界の運命だ。
でも、主流になるものを今から使っておきたい.....と誰もが思うだろう。

だからといって、Deep Learning を避けて通るというのは難しいというか、取り残されてしまう可能性もある。

でも、使う側にとって、フレームワークなので、別のに乗り換えるのに、スクラッチで書き上げたソフトを別の言語に移植するような困難は通常は無いので安心である。
とりあえず、自分の都合の良いのを使い込むべきだろう。
そして、いつ状況が変わってもよいように準備しておくべく、あまりにも使用フレームワークに依存したのを書かないのが良いだろう。

人工知能、機械学習、Deep Learning などは、大学でも普通に教えている情報科学、情報工学の一般教養になりつつある(かな?)


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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、サンプルデータを用意して、実際に一般化線形混合モデルを最尤推定してみるコードを用意しました。 コードはRで書きました。

library('glmmML')

dataSize <- 100
N <- 20

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

createData <- function() {
  x <- runif(dataSize, 0, 5)
  y <- vector(length = dataSize)
  for (i in 1:dataSize) {
    r <- rnorm(1, mean = 0, sd = 2)
    l <- linkFunction(-2, 1, x[i])
    y[i] <- rbinom(1, N, l(r))
  }
  id <- 1:dataSize
  data.frame(x, y, id)
}

data <- createData()
glmmML(cbind(y, N - y) ~ x, data = data, family = binomial, cluster = id)

コードを実行するには、glmmMLパッケージをインストールする必要があります。

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。

サンプルデータは、0 から 5 までの一様分布から x を生成し、標準偏差 s = 2 の正規分布から r を生成し、b1 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、最尤推定すると、b1 = -2b2 = 1s = 2 に近い値を推定してくれるはずです。

コードを実行すると以下のような出力が出力されます。

Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = data,      cluster = id)


              coef se(coef)      z Pr(>|z|)
(Intercept) -2.163   0.4570 -4.732 2.22e-06
x            1.210   0.1604  7.539 4.73e-14

Scale parameter in mixing distribution:  2.144 gaussian
Std. Error:                              0.1654

        LR p-value for H_0: sigma = 0:  5.114e-131

Residual deviance: 311.6 on 97 degrees of freedom       AIC: 317.6

期待通りにパラメータが推定されているのが分かります。

自然言語処理と深層学習,C言語による (424x600).jpg

題名:自然言語処理と深層学習

《C言語によるシミュレーション》

小高 知宏 著

A5、224項、本体2500円

平成29年3月25日 発行

オーム社

ISBN 978-4-274-22033-3



深層学習(Deep Learning)のブームが続いているが、多くの場合は画像系が話題の中心になっている。
深層学習は、何も画像処理だけでなく、自然言語処理にも使える、他にも、たぶん色々なことに使えるであろう。
ということで、少々画像関連ばかりをやるのに飽きてきたというのもあって、自然言語処理の本を取り上げることにした。


本書の最大の特徴は、C言語を使って実装されているということだ。

深層学習というと、今は圧倒的にPythonが多い。それも、多くはフレームワークを使っていて、深層学習がどのように実装されているのとかは全く分からないけれど、とりあえず動かすことができ、何かAI的な動きになっているのを感じられるという本が多いように思う。

この本は、そういう流れとはまったく違う。
C言語で、何もライブラリを利用せず、もろにゼロから作り上げている。

最初の方は準備として基本的な説明がされており、その後、畳み込みニューラルネットワーク(CNN)とリカーレントニューラルネットワーク(RNN)が、それぞれ1つの章として書かれている。
RNNについては、説明している書籍はかなり少ないようなので、参考になるかも知れない。

Cでゼロから作っているし、200ページ余のかなり薄い本なので、扱っている内容は基本的なことだけであるが、実行結果やその説明はかなり丁寧に書かれている。

Cのソースコードおよびデータは、オーム社のダウロードサービスのページから取得可能である。
書籍の実行結果を見ると、Windows上での実行結果が載っていたが、C言語の標準的な部分しか使用していないので、C言語が動くコンピュータなら何でも動作するのではないかと思う。
私はUbuntu/LinuxのGNUのCでコンパイル・実行してみたが、特に問題は無かった。

難点があるのは、Cのインデントである。
インデントの単位が1カラムになっており、読むのにちょっと困る。
印刷する場合、8tabだと都合が悪いかもしれないので、4tabで字下げをして欲しかった。
今では、C Beautifier(Cプログラム整形ソフト) という便利なツールが色々あるので、読みやすいように自動整形してから利用するのをお薦めする。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中で二項分布とと正規分布を混ぜて使う一般化線形混合モデルについて説明されています。 なので、サンプルデータを用意して、様々なパラメータにおける対数尤度の値をプロットしてくれるコードを用意しました。 コードはRで書きました。

dataSize <- 50
N <- 20

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

dmixture <- function(b1, b2, s) {
  function(y, x) {
    l <- linkFunction(b1, b2, x)
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

likelihood <- function(b1, b2, s) {
  d <- dmixture(b1, b2, s)
  function (y, x) {
    d(y, x)
  }
}

logLikelihood <- function(data, b1, b2) {
  function (s) {
    sl <- 0
    ll <- likelihood(b1, b2, s)
    for (i in 1:dataSize) {
      sl <- sl + log(ll(data$y[i], data$x[i]))
    }
    sl
  }
}

createData <- function() {
  x <- runif(dataSize, 0, 5)
  y <- vector(length = dataSize)
  for (i in 1:dataSize) {
    r <- rnorm(1, mean = 0, sd = 1)
    l <- linkFunction(-2, 1, x[i])
    y[i] <- rbinom(1, N, l(r))
  }
  data.frame(x, y)
}

data <- createData()
xl <- "s"
yl <- "Lilelihood"
ss <- seq(0.2, 2, by = 0.2) 
b2s <- seq(-2, 2, by = 1)
legends <- paste0("b2 = ", b2s)
for (b1 in seq(-2, 2, by = 1)) {
  title <- paste0("logit(q) = ", b1, " + b2 * x + r, r ~ N(0, s)")
  plot(0, 0, type = "n", xlim = c(0, 2), ylim = c(-1000, 0), main = title, xlab = xl, ylab = yl)
  for (i in 1:5) {
    ls <- logLikelihood(data, b1, b2s[i])(ss)
    lines(ss, ls, type = "l")
    points(ss, ls, pch = i)
  }
  legend("topleft", legend = legends, pch = 1:5)
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。

サンプルデータを生成するには、まず 0 から 5 までの一様分布から x を、標準偏差 s = 1 の正規分布から r を生成します。 その後、b1 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、様々なパラメータについてモデルを試すと、b1 = -2b2 = 1s = 1 の時に最もデータへの当てはまりがよくなるはずです。 つまり尤度が最大になるはずです。

コードを実行すると5枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_likelihood.png

図を見ると、期待通りに b1 = -2b2 = 1s = 1 の場合に尤度が最大になっているのが分かります。 実際の問題では、これらのパラメータを最尤推定すればいいわけです。

出力結果は非常に単純なものである。
座標値(位置)に対する変位電流など電磁気学的な量を計算するだけである。
要するに、横軸が位置xで、縦軸が f(x) というだけの簡単なグラフを描きたいのだが、1980年初期にはグラフを描くのはなかなか大変なので、単に数値だけがだらだらと出力ファイルになって出てくるものである。
もちろん、今は、計算結果をExcelに読み込んで、グラフを描いてシミュレーション結果を評価しているようだった。

このプログラム、1980年代のコンピュータを想定しているので、配列などのサイズが非常に小さく制限されている。
当時は、大型コンピュータといっても、せいぜい1MB程度しかメモリが使えないのが普通であった。
そういう制限の下で、ごちゃごちゃと計算していたのだ。
せっかくWindowsになるのだから、グラフをもっと細かく精密に出したいという当然の要望があった。

実際のプログラムでは、サイズの上限の数値がプログラム中に散在していたのだが、それを1つの定数変数として宣言することで、1箇所直せばサイズを変更できるようにした。
グラフの細かさは100倍くらい楽勝と思われたが、何とか現状の2倍にして欲しいというので、10倍に設定した。

さて、この機能強化により、より細かい精密な評価ができるようになる..........はずだった。
x軸上の測定間隔を半分にしたら、出発点が同じなら、1つおきに元の値と完全一致する筈だ。
でも、一致しなかったのだ。
まあ、計算誤差というのもあるから、ある程度誤差が出るのは仕方がないのであるが、明らかにかなり違う。

想定外の問題が出てしまった。どうしよう。
ということで、デバッガを2つ立ち上げて、計算間隔を丁度倍違うようにして、同じx値に対する計算が始まる箇所で止めを、変数の値を全部調べ比較したのだ。
どこで何が違うのかを調べた。
違っても問題ないのかもしれないし、プログラムの意味を知らずに調べているので苦しい。
それでも、デバッガでステップ動作で調べていると、徐々に分かってくるものである。

やっと分かったのは、やはり初期化の問題だった。
内部ループの最初で行うべき初期化が、なぜか省略されていた。
電磁波の反射などの影響を延々と足しこんでいるらしい配列だったのだが、1回前の測定点の値にそのまま足しこんでいるようだった。
恐る恐るループの最初で疑わしい複素配列を初期化すると、結果がまるで変わった。

数字の羅列では良く分からないし、でもExecelに結果を読み込むのも面倒なので、Pythonで比較グラフを表示するツールを作った。
青線が元の計算結果、赤が内部ループにゼロクリアを入れた場合の計算結果である。
計算は、右から左に向かって進んでいく。
振動が消えてしまったが、これで良いのだろうか、悪化したのか、情報がないので分からないので、お客さんに聞いてみた。
compare.png 昔は、赤い線のようなグラフになっていたそうだ。
このグラフが出るようになってから、急に資料の発掘が進み、30年以上前の資料でグラフが載っているものがあり、振動がないことがわかり、昔の状態に復元できたようで、魔法使いのように扱われ、とても感謝された。

ということで、一件落着になったのだが、こんなに上手くいくことは珍しい。
非常に古いプログラムで、資料が無い、作った人が居ない、その他問題が多いプログラムを扱うのは拒絶するのが正しい。
それでも諸般の事情、大人の事情などで作業しないといけないこともあるが、その場合は細心の注意で望まれたい。

このような作業をおこなうと、非常に力が身に付くことも確かである。
ハッカー養成ギブスとしては良いかも知れない。

(完)

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中でポアソン分布と正規分布を混ぜた分布について説明されています。 なので、ポアソン分布と正規分布を混ぜ合わせた混合分布をグラフにしてくれるコードを用意しました。 コードはRで書きました。

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

distribution <- function(s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dpois(y, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

xl <- "y"
yl <- "Probability"
xs <- 1:5
ys <- 0:10
b1 <- 0.1
legends = paste0("x = ", xs)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    title <- paste0("lambda = exp(", b1, " + ", b2, " * x + r, r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 0.5), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      l <- linkFunction(b1, b2, xs[i])
      d <- distribution(s, l)
      lines(ys, d(ys), type = "l")
      points(ys, d(ys), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。 r は平均 0 、標準偏差 s の正規分布から生成されます。 r を生成する正規分布の標準偏差が0に近い数の時は、これはポアソン回帰のリンク関数とほとんど同じです。 混合分布はポアソン分布と正規分布をかけて r で積分したものです。

コードを実行すると20枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_poisson.png

この図は b1 = 0.1b2 = 0.5s = 1 の場合です。 ポアソン分布よりも広く分布した分布が得られているのが分かります。

プログラムの内容の詳細は書く訳にいかないのだが、技術計算で、電磁気学に関するシミュレーションである。
数学的には、複素関数、複素級数、偏微分(マックスウェルの方程式)などが使われている。
たぶん今だったら、MatLab(Octave)でちょっと書けば済みそうな内容であるが、困ったことに内容が分からない。

日々科学技術計算をやっている訳ではないので、どうなるか分からないが、まあ調べれば何とかなるのではと思った次第である。
科学技術計算とかいってずいぶん難しいことをやっていると予想して勉強してから行くと、ほとんど科学技術とは無縁のプログラムであることが多いので、たぶんこれも大丈夫と予想したのであった。
それでも、一応、いざというときのために、関連情報などをネットで軽くは調べておいた。

このくらい古くなると、作ったときの資料も無ければ、作った人も居ない。既に彼方の世界の人になっている。
でも、ソースプログラムがあるから恵まれているほうである。
ソースプログラムが無い移植などの仕事もあるのだが、そういうのに比較すれば易しい。(比較がオカシイかも)
かなり面倒な計算をしているのだが、詳細な計算内容はとても追いかけるのは無理で、核心部分はブラックボックスのままであるが、それでも移植は可能であり、全く同じ動作をするかを確認する手段を考えて、引き受けた。

以上は、元のプログラムに何も問題がなければ、全く同じ動作をすれば移植のテスト完了となるのであるが、世の中そんなに甘くないのである。

プログラムは9本に分かれていて、技術計算をいっぱいしているのは2本だけである。
ある特殊な装置の電磁波の影響を計算するプログラムで、水平方向と垂直方向で別々のプログラムに分かれている。
垂直方向が5本、水平方向が4本の構成である。
計算以外のプログラムは、入力データを作るためなのだが、昔々のプログラムなので全てCUIになっている。

まずはコンパイルエラーを取るのだが、ワーニングがかなり出てうっとうしい。
コンパイル時にワーニングを出なくすることは可能だが、問題があったときに大事なワーニング情報を見失う可能性がある。
それで、ワーニングが出なくなるようにするFORTRANのソースコード・コンバーターをちょっとPythonで書いた。

実際のコンパイルオプションはこんな感じである。
gfortran -Wall -pedantic -std=f95 -fbounds-check -fall-intrinsics -Wuninitialized 
	-ffptrap=invalid,zero,overflow -o HOGE.exe HOGE.for
FORTRANのコンパイラオプションを決めるのに、東京大学のFortranデバッグオプションのページを参考にした。
まだFORTRANをどんどん使っている研究室があるんだ。

コンパイルが済んだら、とりあえず走らせてみるのだが、特殊な技術計算用のプログラムなので、データを用意するのが困難なので、データをファイルでもらって、それを読み込んで走らせることにした。

最初は、実行時エラーが出てしまった。
変数に値をセットする前に参照していると叱られた。
つまり、配列の初期化がされていない。
これがやたらに多かったのだが、昔のFORTRANの一部には、配列宣言をしたら初期化も自動で行われるもの(コンパイラのオプション指定)もあったためだろう。
ということで、初期化されていない配列の初期化をかなり行った。

一応、初期化を済ませたら、計算処理プログラムが何とか走り出した。
コンパイル時のオプションもどんどん減らして、最終的にはソースとオブジェクトの指定だけまで減らした。
やれやれである。

入力データを、まだ旧ソフトが動くコンピュータ上で実行してもらい、入力データと出力データをもらって、移植したプログラムでも実行し、2つの出力ファイルを比較した。
完全一致したら、移植が完璧ということになる。
出力ファイルはテキストファイルなので、ファイル比較を行い、一文字も差がないことを確認すればよい。
もし、旧プログラムにバグがあれば、バグも完全に再現できているという素晴らしい移植ということになる。

技術計算なので、どうしても数値計算誤差がでてくる。
計算結果は固定小数点数で、有効数字も3,4桁程度であるので、少々の計算誤差があっても、結果は同じになる。
それでも、2つ差異が出てしまった。

1つは、0.00と-0.00の差である。
0に非常に近い値のとき、計算誤差により正になったり負になったりする。
これは値としては同じと考えられるので、OKとした。

計算結果の最後の1行が出たり出なかったりした。
これは、リミットの計算を計算誤差を無視して、単に不等号で判断していたのが原因で、ちょっと余裕をもたせるようにして対処した。

これで、同じ入力データを新旧のプログラムで走らせた結果の2つの出力ファイルが一致したことになる。
とりあえず移植完了。

移植以外に、いくつかの機能UPの要望があったのだ。
どれも比較的簡単なものと思って気楽に引き受けたのだが、思わぬ展開になるのであった。

次回につづく。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中でポアソン分布と正規分布を混ぜた分布について説明されています。 なので、ポアソン分布と正規分布を掛けた値をグラフにしてくれるコードを用意しました。 コードはRで書きました。

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

distribution <- function(y, s, linkFunction) {
  l <- match.fun(linkFunction)
  function(r) {
    dpois(y, l(r)) * dnorm(r, mean = 0, sd = s)
  }
}

xl <- "r"
yl <- "Probability"
ys <- seq(0, 8, by = 2)
rs <- seq(-5, 5, by = 0.5)
b1 <- 0.1
legends = paste0("y = ", ys)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    l <- linkFunction(b1, b2, 2)
    title <- paste0("lambda  = exp(", b1, " + ", b2, " * 2 +  r), r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.2), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      d <- distribution(ys[i], s, l)
      lines(rs, d(rs), type = "l")
      points(rs, d(rs), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。 r は平均 0、標準偏差 s の正規分布から生成される変数です。 あとは、正規分布から r が生成される確率と、平均 λ のポアソン分布を掛け合わせれば、欲しい値になります。

コードを実行すると20枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixed_poission.png

y によってグラフの山の大きさが変わってくるのが分かります。 混合分布はこれを r で積分したものです。 なので、グラフ上で面積の大きな y が確率の高い分布になります。

今回は、ちょっと別の昔話、ではなくて今年の話。

色々なことをやっていると、突然変な仕事が降ってくることがある。
FORTRANという、たぶん最近の人は見たことも聞いたこともないプログラムの話だ。

昔々、COBOLとFORTRANというプログラミング言語が主に使われていた。
COBOLが事務処理用で、FORTRANが技術計算用である。
当時は、もちろんネットもないし、WEBなど全然ない時代である。
以下、コンピュータ考古学の話である。

さて、今回、ずいぶん昔にFORTRANで書かれたプログラムを、Windows10で動くようにして欲しいというので持ち込まれた。
いや、泣きつかれたというほうが正しい。
平成ではなく、昭和の時代に書かれたもので、80年代に開発が行われ、当時、大企業のコンピュータセンターの大型コンピュータで動いていたものだ。
といっても、当時のコンピュータは呆れるほど非力で、今のパソコンはもちろん、Raspberry Pi にも劣る能力である。
でも、そこで計算していたものは、今でも高価な機械に関するシミュレーションプログラムで、30年以上使い続けられていて、今後も使う予定のものだ。

まあ、諸般の事情で逃げられないので引き受けたのだが、FORTRANを最後に使ったのは確か1982年だと思う。
その後は、C, C++, Java, Python などなど言語放浪を続けている。
なので、覚えているか全然自信がなかったのだが.....

それで、まずはWindows上で動くFORTRANを探した。
持ち込まれたFORTRANのプログラムが、非常に古い書き方をされたもので、まだパンチカード時代を偲ばせるものだ。
要するに古文書である。
FORTRANも長い間に仕様が色々変遷してきたのだが、1970年代の FORTRAN JIS 7000 という仕様で書き始められ、手を入れた部分だけがやや最近の書き方になっていた。
こういうプログラムでもちゃんとコンパイルできるものを探さないと修正作業が面倒になる。
とくに古いプログラムを扱うとき、まるで国宝を扱うように、可能な限り原型のままで動くようにしないといけない。
ちょっとでも手を入れると、何が起きるか怖いのである。
そのため、古い書き方でも対応してくれ、Windows10も対応しているFORTRANを探した。


その結果見つけたのが、GNUのFORTRANであり、gfortran(gfortran -- the GNU Fortran compiler, part of GCC)と呼ばれるものだ。
実際にWindows環境にインストールする方法はいくつかあるようだが、MinGWを使った。
これだと、古い書き方から最新の書き方まで混在しても大丈夫なようだった。
とりあえずコンパイルできることは分かったが、それだけでは動かないようだった。

また、デバッグにgdbを使うので、Cのデバッグと同じ感覚でデバッグでき、便利であった。
このデバッガが1980年頃に使えたらどんなに便利だったかと思う。

さて、実際の移植作業がどうだったかについては、次回紹介しよう。
プーリングのサイズを大きくするとすぐにダメになるのではと思ったが、そんなことはないらしいのが前回わかった。
それで、プーリングサイズをどんどん大きくし、テスト精度が最も良くなるのはどのあたりかを調べてみた。

グラフの右には、プーリングのサイズ、そのサイズでテスト成功率が最大になったepochとテストの最大精度を示している。

accuracy-CNN2s.pngプーリングなし  3 0.6667

accuracy-CNN5x5s.png プーリング 5x5  6 0.7354

accuracy-CNNp6x6s.png プーリング  6x6  7 0.7412

accuracy-CNNp7x7s.png プーリング 7x7  8 0.7436

accuracy-CNNp8x8s.png プーリング 8x8  6 0.7402

accuracy-CNNp10x10s.png プーリング 10x10  6 0.7307


この実験では、7x7のプーリングで、プーリングなしと比較して7.69%(66.67→74.36)テスト精度が向上した。
プーリングだけでこれだけ向上するということは、プーリングは何らかの形で使うのが良さそうだ。

上のグラフを眺めると、もう1つの変化が見えてくる。
緑がテスト精度であるが、青の学習精度の方に注目してみよう。
プーリングなしだと、10epochあたりで98%に達しているのだが、プーリングサイズを上げると、学習段階での精度の上昇ペースがかなり鈍るようである。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、二項分布と正規分布を混ぜ合わせた混合分布をグラフにしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

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

distribution <- function(s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

xl <- "y"
yl <- "Probability"
xs <- c(-5,-2,0,2,5) 
ys <- 0:20
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)) {
    for (s in c(0.1, 0.5, 1, 1.5, 2)) {
      title <- paste0("logit(q) = ", b1, " + ", b2, " * x + r, r ~ N(0, ", s, ")")
      plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 0.5), main = title, xlab = xl, ylab = yl)
      for (i in 1:5) {
        l <- linkFunction(b1, b2, xs[i])
        d <- distribution(s, l)
        lines(ys, d(ys), type = "l")
        points(ys, d(ys), pch = i)
      }
      legend("topright", legend = legends, pch = 1:5)
    }
  }
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。 r を生成する正規分布の標準偏差が0に近い数の時は、これはロジスティック回帰のリンク関数とほとんど同じです。

混合分布は二項分布と正規分布をかけて r で積分したものです。 今は単純に r-100 から 100 まで 0.1 刻みで変化させて、面積の合計を求めておきました。 コードで言うと以下の部分です。

    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }

これで様々な b1b2s について混合分布のグラフをプロットできます。 コードを実行すると180枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_distribution.png

この図は b1 = 0.5b2 = 0.5s = 1 の場合です。 二項分布とはたいぶ違った分布が得られているのが分かります。

このアーカイブについて

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

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

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

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