ほくそ笑む

R言語と統計解析について

ベイズ統計の入門書が出版ラッシュなのでまとめてみた


【宣伝】2016/09/14

このページに来た方へ。あなたが求めている本はこれです。

まずこれを予約してから下記を読むといいです。
【宣伝終】

最近、ベイズ統計の入門書がたくさん出版されているので、ここで一旦まとめてみようと思います。

1. 基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門 (2015/6/25)

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

データ分析業界ではかなり有名な豊田秀樹先生の本です。
難易度としては易しいほうだと思います。
最新のハミルトニアンモンテカルロ法について容易に解説されているところが特徴です。
Stan ユーザでここらへん分かっていないのなら必読の本だと思います。
勉強会も開催されていて、スライドも公開されていますので、これらを見ながら読み進めていくといいと思います。

(追記:この本は有識者によると誤りが多いそうです)

2. 岩波データサイエンス Vol.1 (2015/10/08)

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

岩波から出版された「岩波データサイエンス」、シリーズの一冊目の特集が「ベイズ推論とMCMCフリーソフト」でした。
難易度としては簡単なものから難しいものまで多岐に渡っています。
ベイズ推定を行う際の実践的な内容が盛りだくさんというのが特徴です。
特筆すべきは、サポートページの充実ぶりです。購入を検討している人は、一度このサイトに目を通しておくといいと思います。

3. ベイズ計算統計学 (2015/10/10)

ベイズ計算統計学 (統計解析スタンダード)

ベイズ計算統計学 (統計解析スタンダード)

結構ガチな入門書。数式を使って入門したい人向け。
著者は下記の資料を書いた方なので、購入を検討している人は、一度見ておくといいと思います。

4. 完全独習 ベイズ統計学入門 (2015/11/20)

完全独習 ベイズ統計学入門

完全独習 ベイズ統計学入門

完全独習 統計学入門 を書いた人です。この本を読んで気に入った人ならぜひ読むべきでしょう。

著者による書籍紹介はこちら。

5. 見えないものをさぐる―それがベイズ (2015/11/28)

見えないものをさぐる―それがベイズ: ツールによる実践ベイズ統計

見えないものをさぐる―それがベイズ: ツールによる実践ベイズ統計

読んでませんが、Excelベイズをやるらしい。

解説する上で最小限必要とする数式は掲載しますが、ベイズ法で大きな障害となる「計算が難しい」という問題点をツール「Weka」や「Excel」を積極的に使用して簡略化し、データ分析の敷居を低くすることで、「理論より実践」を目指します。

6. ベイズ法の基礎と応用 (2016/1/12)

まだ発売されていませんが、東工大の間瀬先生も出されるようです。
間瀬先生といえば、R ヘルプの日本語訳で有名です。詳しくは下記記事。

7. まずはこの一冊から 意味がわかるベイズ統計学 (2016/2/24)

まずはこの一冊から 意味がわかるベイズ統計学 (BERET SCIENCE)

まずはこの一冊から 意味がわかるベイズ統計学 (BERET SCIENCE)

まだ発売されていませんが、こういうのも出るらしい。

具体的な例題も豊富な本書は、「専門書はとっつきにくいけれども数式のない入門書では物足りない」「これから本格的に学んでいこう」という読者にとって最適な「入門書」です。

8. 医薬データ解析のためのベイズ統計学 (2016/2/25)

医薬データ解析のためのベイズ統計学

医薬データ解析のためのベイズ統計学

“Bayesian Biostatistics”の翻訳
WinBUGSやSASなどの統計ソフトを用いたデータ解析を詳しく紹介
医薬統計分野の実務家にとっても参考となる

9. 身につく ベイズ統計学

身につく ベイズ統計学 (ファーストブックSTEP)

身につく ベイズ統計学 (ファーストブックSTEP)

ベイズ統計の基本的な考え方や活用について、初学者向けテキストとして、ていねいに解説していきます。
ベイズ統計学をこれから勉強したい大学生、ビジネスで実際に活用したい社会人に最適です。

10. はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

独自の解釈が多く、初心者が読むと混乱するためお勧めしない。

既刊

ベイズ統計の入門書はすでにたくさん出ていて、いいものがありますのでついでに紹介します。

1. 史上最強図解 これならわかる!ベイズ統計学 (2012/2/21)

史上最強図解 これならわかる!ベイズ統計学

史上最強図解 これならわかる!ベイズ統計学

最初に読む本の定番になっているのではないでしょうか。わかりやすいです。
最後に出てくるこの図が個人的にはお気に入りです。


2. 図解・ベイズ統計「超」入門 あいまいなデータから未来を予測する技術 (2013/12/18)

上記の本と内容は基本的に変わりませんが、こちらは読み物要素が入っているので、文系向きかも。

3. 異端の統計学 ベイズ (2013/10/23)

異端の統計学 ベイズ

異端の統計学 ベイズ

こちらは完全に読み物。長いですが、文章を読むのが好きな人にはいいかも。私は途中で挫折しました。

4. Think Bayes ―プログラマのためのベイズ統計入門 (2014/9/6)

Think Bayes ―プログラマのためのベイズ統計入門

Think Bayes ―プログラマのためのベイズ統計入門

未読ですが、Think Stats が面白かったのでこちらも良いかなと。
プログラマ向けの入門書で、Python で書きながら学ぶ感じです。

(追記:個人的感想に過ぎませんが、ちょっといまいちでした)

ということで、英語版の PDF が公式ホームページから無料でダウンロードできます。

5. 入門ベイズ統計―意思決定の理論と発展 (2008/6)

入門ベイズ統計―意思決定の理論と発展

入門ベイズ統計―意思決定の理論と発展

私が最初に学んだ本なので思い入れが深いですが、分かりやすさでは上で紹介した入門書の方が良いかもしれません。
ただ、非常に面白い本なので、ある程度学んだあとでいいので、ぜひ手に取ってほしいです。

6. 計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (2005/10/27)

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

Amazon レビューでも書かれていますが、伊庭先生の一章は必読です。

R言語でWebアプリケーションを作るためのチュートリアルを翻訳しました

RStudio社が開発した Shiny パッケージは、R言語で簡単に Web アプリケーションを作るためのフレームワークを提供します。

この Shiny による Web アプリケーションの開発方法を学ぶには、公式のチュートリアルを読むのが一番です。

しかし、公式は英語なので、読むのがしんどいです。
そこで、チュートリアル全文を日本語に翻訳しました。
訳文は Qiita で公開しています。このページは目次として活用していただければと思います。

Shiny チュートリアル目次

7 つのレッスンからなるこのチュートリアルは、R プログラマを Shiny 開発者へと導きます。
1 つのレッスンは 20 分ほどで終了し、各レッスンごとに新しい Shiny スキルを 1 つ学ぶことができます。
すべてのレッスンを終えたとき、あなたは Shiny アプリを構築し、公開することができるようになっています。


各レッスンにはエクササイズ(練習問題)が用意されています。
はやる気持ちは分かりますが、このエクササイズを読み飛ばしてはいけません。
エクササイズを行うことによって重要な学びを得ることができます。
我々はエクササイズが重要な意味を持つように、このチュートリアルを作成しました。


それでは、レッスン1 をクリックして、Shiny の世界へ飛び込みましょう!

【速報】神ハドリーの名著 "Advanced R" の翻訳が出ます

ビッグニュースです。

統計言語 R の業界において神とあがめられているハドリー・ウィッカム
彼は ggplot2 や dplyr などの数々の超便利パッケージを作成したことで知られていますが、その R に対する深い知識をもとに書かれた "Advanced R" という名著を昨年秋に出版しました。

Advanced R (Chapman & Hall/CRC The R Series)

Advanced R (Chapman & Hall/CRC The R Series)

本書は、これまで複雑でわかりにくいとされてきた R 言語を、エッセンスにしぼり、体系的にここまで綺麗にまとめることができたのかと、思わずうなるほどの出来で、アメリAmazon で星 4.8 という超高評価をつけられている名著です。

この本の日本語訳が、ついに出版されます!

R言語徹底解説

R言語徹底解説

訳者はこれまで数々の R 関連書を翻訳してきた石田基広先生と、R 界隈で知らない人はいないお三方。
翻訳の品質については折り紙付きです。
値段はややお高いですが、はっきりいって、この本一冊あればもう他の R 解説本はいらないんじゃないかと思えるほどの内容なので、むしろお得です。

このニュースに、R 界隈の著名人たちが一斉に反応しています。

これを読まずして R を語ることは許されません。
今後 10 年、この本が R ユーザーのバイブルとなることはほぼ間違いないでしょう。

というわけで、みなさん、今すぐ予約しましょう! ぼーっとしてたらすぐに売り切れてしまいますよ!

R言語徹底解説

R言語徹底解説

可視化で理解する「負の二項分布」

みどりぼんでカウントデータの過分散対策のために使われると書かれている負の二項分布ですが、Wikipediaの説明を読んでもよく分かりません。

そこでおススメなのが、このスライドです。

ようするに、負の二項分布は、\lambda がガンマ分布に従うようなポアソン分布だと思えばだいたい OK みたいです。

今日はこれを可視化してみます。

負の二項分布(Negative Binomial Distribution)

負の二項分布はパラメータを 2つ持ちます。成功回数を表す r と成功確率を表す p です。

統計言語 R には負の二項分布に従う乱数を生成する関数 rnbinom() があり、これらのパラメータはそれぞれ引数 size と prob に対応しています。

したがって、r=4, p=0.3 の負の二項分布は次のようにして描画できます。

negative_binom_values <- rnbinom(10000, size = 4, prob = 0.3)
barplot(table(negative_binom_values)[1:21])

これに対応するのは、\lambda がガンマ分布  {\rm Gamma}(r, p/(1-p)) に従うポアソン分布です。
このガンマ分布から \lambda を生成し、各 \lambda に従うポアソン分布から乱数を生成という手順でやってみましょう。

gamma_values <- rgamma(10000, 4, 0.3/(1-0.3))
poisson_values <- sapply(gamma_values, function(lambda) rpois(1, lambda))
barplot(table(poisson_values)[1:21], main="λがガンマ分布に従うポアソン分布")

見事に同じような分布の形状になりました。

少しわかりにくいので、この様子をアニメーションにしてみましょう。*1

最下段がガンマ分布であり、このガンマ分布からサンプリングされた値を赤い線で示しています。
この値をパラメータ  \lambda として持つポアソン分布が中段のグラフです。
このポアソン分布からサンプリングされた値を赤い線で示し、この値からなるヒストグラムを最上段に描いています。
毎回異なるポアソン分布から値がサンプリングされており、それらの値が負の二項分布を形成していく様子がわかります。

混合ポアソン分布

同じ状況を別の角度から見てみましょう。
ガンマ分布からサンプリングされた値をパラメータ \lambda として持つポアソン分布を複数作成します。

library(ggplot2)

df <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:16]))
ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + facet_wrap(~ lambda, nrow=4)

これらの分布を合成したものが負の二項分布です。

ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none")

上図は 16個のポアソン分布を合成したものですが、1000個ほど合成すればほとんど負の二項分布と一致することがわかります。

library(gridExtra)

df2 <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:1000]))
g1 <- ggplot(df2, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none") + xlab("混合ポアソン分布(n=1000)") + ylab("")
g2 <- qplot(x=0:20, y=dnbinom(0:20, 4, 0.3)) + geom_bar(stat="identity") + xlab("負の二項分布") + ylab("")
grid.arrange(g1, g2, ncol=1)

この様子も無駄にアニメーションにしてみましょう。


まとめ

みどりぼんでは、個体差を持つカウントデータに対処するために GLMM(一般化線形混合モデル)を導入しています。

この本では、GLMM の例として、R の glmmML パッケージを使って「パラメータが正規分布に従う二項分布」の推定を行っています。
これに対して、負の二項分布は「パラメータがガンマ分布に従うポアソン分布」だと考えれば同じ枠組みでとらえることができると思います。

R で負の二項分布を用いた GLMM は、MASS パッケージの glm.nb() 関数によってできます。

負の二項分布がパラメータによってどのような形状をとるかの感覚を知りたい場合は、次のサイトが便利です。

みなさまの理解の一助になれば幸いです。

追記

ホクソエムの切り込み隊長 teramonagi さんが glm.nb() の使い方について詳しく書いてます。

*1:アニメーションの作成コード:https://gist.github.com/hoxo-m/3060e2fe269ff1fdf4a3

Stan でパラメータに大小関係の制約をつける

EC サイトなどでは、1 ページの読み込み時間の平均値が短いユーザほどコンバージョン率が高くなるという話があります。
この原因については、サクサク読み込みが終わる速いサイトの方がユーザは気前よく購入する、というわけではなく、単に速い回線を持つほどお金持ちなので購入しやすいという身も蓋もない話もあります。
今回はこの原因についてではなく、制約条件が付いた時のパラメータ推定の話をします。
とりあえず擬似データを作成しましょう。

true_cvr <- c(0.209, 0.126, 0.096, 0.093, 0.086, 0.077, 0.067, 0.057)

load_time <- c("0-1", "1-3", "3-7", "7-13", "13-21", "21-35", "35-60", "60+")
session <- c(1000, 6000, 4000, 1500, 700, 500, 200, 150)
set.seed(71)
cv <- unlist(Map(function(n, p) rbinom(1, n, p), session, true_cvr))

data.frame(load_time, cv, session)
  load_time  cv session
1       0-1 208    1000
2       1-3 769    6000
3       3-7 366    4000
4      7-13 142    1500
5     13-21  54     700
6     21-35  35     500
7     35-60  17     200
8       60+  11     150

ロード時間ごとに真のコンバージョン率(true_cvr)を設定して、コンバージョン数(cv)を二項分布からの乱数(rbinom)として生成しています。

このときのコンバージョン率を Stan で推定してみましょう。

// load_time.stan
data {
  int<lower=0> N;
  int<lower=0> cv[N];
  int<lower=0> session[N];
}
parameters {
  real<lower=0, upper=1> cvr[N];
}
model {
  cv ~ binomial(session, cvr);
}
library(rstan)
model <- stan_model("load_time.stan")
standata <- list(N = length(cv), cv = cv, session = session)
result <- sampling(model, data = standata, chains = 3)

この結果を見やすく可視化してみます。

library(tidyr)
df <- rstan::extract(result, pars="cvr")$cvr %>%
  data.frame %>%
  setNames(load_time) %>%
  gather(load_time, cvr)

cvr_hat <- colMeans(rstan::extract(result)$cvr)
point_df <- data.frame(load_time, true_cvr, cvr_hat) %>%
  gather(load_time, cvr) %>%
  setNames(c("load_time", "cvr_type", "cvr"))

library(ggplot2)
ggplot(df, aes(x = load_time, y = cvr)) + 
  geom_violin() + ylim(0, NA) +
  geom_point(data = point_df, aes(x=load_time, y=cvr, color=cvr_type))

パラメータの推定値(cvr_hat)と真値(true_cvr)には開きがあります。

制約条件の追加

ここからが本題です。
今、「読み込み時間が短いユーザほどコンバージョン率が高くなる」という条件が分かっています。
この条件を考慮してパラメータを推定することはできないでしょうか?
Stan では、パラメータを ordered 型で定義することによって、パラメータに順序関係の制約を加えることができます。
このような制約条件を付けた場合の Stan コードは次のようになります。

// load_time_ordered.stan
data {
  int<lower=0> N;
  int<lower=0> cv[N];
  int<lower=0> session[N];
}
parameters {
  ordered[N] cvr_rev;
}
transformed parameters {
  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- inv_logit(cvr_rev[N - i + 1]);
  }
}
model {
  cv ~ binomial(session, cvr);
}
library(rstan)
model <- stan_model("load_time_orderd.stan")
standata <- list(N = length(cv), cv = cv, session = session)
result <- sampling(model, data = standata, chains = 3)

非常にきれいにパラメータを推定することができました。

注意点

まず、parameters ブロックにおいて、ordered 型のパラメータを作成しています。

  ordered[N] cvr_rev;

この際、ordered 型は、「小さい順」という制約であり、求めたいパラメータは大きい順なので、これを逆順にする必要があります。
この逆順にするという操作を行っているのが transformed parameters ブロックです。
ここで、新しいパラメータ cvr に cvr_rev を逆からひとつひとつ代入する、つまり、

  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- cvr_rev[N - i + 1];
  }

のようにすれば良いように思いますが、それではうまくいきません。
ordered 型には lower や upper のような最小値・最大値の制約がつけられないため、−∞〜+∞までの値を取ります。
しかし、cvr には 0〜1 という制約が付いているで、cvr_rev を適当にサンプリングすると、ほとんどがはみ出してしまいます。
このような場合に便利なのが inv_logit() です。
これは、logit 関数の逆関数(logistic 関数!)で、−∞〜+∞ の値を 0〜1 にマッピングしてくれます。
したがって、

  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- inv_logit(cvr_rev[N - i + 1]);
  }

と書けば、cvr_rev がどのような値を取ろうが、cvr は 0〜1 の値となり、cvr の制約条件に違反しなくなります。

まとめ

Stan でパラメータに大小関係の制約を付ける方法を紹介しました。
このように、パラメータ間の制約を柔軟に加えることができるのも、ベイズ統計モデリング(Stan)の魅力かなと思います。
Enjoy!

協調フィルタリングについて発表しました

某所で協調フィルタリングについて発表してきました。資料を公開します。
機械学習とかの知識を全く持たないエンジニアさん向けの導入資料です。


協調フィルタリングは、レコメンド(推薦)に使われるアルゴリズムの一つです。

レコメンドは適用できる場面が多く、効果も分かりやすいので、実務では広く使われています。

レコメンドシステムの概要については、神嶌先生の論文が非常に分かりやすく、参考にさせていただきました。

また、協調フィルタリングについては、下記 2 つの記事を参考にさせていただきました。

みなさまのご理解の一助となれば幸いです。

Shiny アプリを収集するサイトを作った

RStudio の作成した Shiny パッケージは、R 言語で Web アプリを簡単に作成するためのフレームワークを提供します。

また、RStudio は shinyapps.io というホスティングサービスを提供しており、Shiny で作成した Web アプリを簡単にアップロードして動かすことができます。

現在、shinyapps.io には様々な Shiny アプリがアップロードされています。
例えば、次のアプリは、入力中の文章の次の単語を予測するアプリです。

テキストエリアに "to be or not to be that is the" まで入れてみると、"question" と出てきます。

また、id:ksmzn さんが作った、次のアプリも Shiny 製です。

いろいろな確率分布の形状を、パラメータを変更しながら確認できるため、非常に便利です。


このように、shinyapps.io には、日々様々なアプリが作成されてはアップロードされています。
しかし、shinyapps.io のサイトでは、自分のアプリを管理することはできても、他の人がどんなアプリを作ったのかを知ることはできません。
せっかく作られたアプリケーションも、誰にも知られずにウェブの藻屑となって消えてしまっている状況です。
また、気に入ったアプリがあっても、あとで探すときに記憶があいまいで苦労するということもあります。

  • 他の人がどんなアプリを作っているか知ることはできないでしょうか?
  • shinyapps.io に新たに追加されたアプリを知る方法はないでしょうか?
  • 気に入った Shiny アプリを記録してまとめておくことはできないでしょうか?


というわけで、作りました。

shinyapps.io に登録されたアプリを収集するサイトです。

原理的には、Google 検索にひっかかった Shiny アプリを片っ端から登録してます。


トップページはアプリが見つかった順番で並んでおり、新しいものほど上に表示されます。
画像をクリックすると、アプリの詳細ページに行き、そこのリンクからアプリにアクセスすることができます。

また、右上のログインボタンから Google アカウントでログインすると、アプリをお気に入りに入れることができるようになります。
お気に入りに入れるには、詳細ページの星マークを押します。

お気に入りに入れたアプリはマイページ(ログインしたら右上にボタンが出ます)で一覧表示できます。


これにより、

  • shinyapps.io に新しく登録されたアプリを知ることができる
  • 作者名をクリックすれば、その人が作ったアプリの一覧が見れる
  • 気に入ったアプリを登録して自分だけの一覧を作れる

ということができるようになりました。
やったね!


また、このサイトに登録された新規 Shiny アプリの情報を、Twitter で流しています。

ぜひフォローして最新情報を追ってください!


というわけで、Shiny アプリを収集するサイト、Shinyapps Gallery を作りました。ぜひ使ってみてください。

それではよき Shiny ライフを!