StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

MCMCサンプルを{dplyr}で操る

RからStanやJAGSを実行して得られるMCMCサンプルは、一般的に iterationの数×chainの数×パラメータの次元 のようなオブジェクトとなっており、凝った操作をしようとするとかなりややこしいです。

StanとRでベイズ統計モデリング (Wonderful R)』のなかでは、複雑なデータ加工部分は場合によりけりなので深入りしないで、GitHub上でソースコードを提供しています。そこでは、ユーザが新しく覚えることをなるべく少なくするため、Rの標準的な関数であるapply関数群を使っていろいろ算出しています。しかし、apply関数群は慣れていない人には習得しづらい欠点があります。

一方で、Rのデータ加工パッケージとして、%>%によるパイプ処理・{dplyr}パッケージ・{tidyr}パッケージがここ最近よく使われており、僕も重い腰を上げてやっと使い始めたのですが、これが凄く使いやすい。%>%selectfiltermutategroup_bysummarize*_joingatherspreadだけをまずは覚えればほとんど不自由しませんでした。これらがないともう他の言語に移れないレベルです。これらのパッケージの練習のおかげで、ややこしいMCMCサンプルの処理についても、こんな感じでやれば毎回ウンウン唸らずに統一的に操作できそうかなぁ、というところまで来ましたので簡単にメモします。

* * *

手始めに以下の図を描いてみます。

この図はパラメータごとにMCMCサンプルの中央値と95%CIを表示した図です。{ggmcmc}パッケージや{bayesplot}パッケージに含まれる関数を使うと一撃で描くこともできます。しかし、練習のため自分で算出して作図します。

library(rstan)
library(ggmcmc)
library(dplyr)

data <- list(J=8, y=c(28,  8, -3,  7, -1,  1, 18, 12), sigma=c(15, 10, 16, 11,  9, 11, 10, 18))
model_code <- readr::read_file(url('https://raw.githubusercontent.com/wiki/stan-dev/rstan/8schools.stan'))
fit <- stan(model_code=model_code, data=data, seed=1234)

d_mcmc <- ggs(fit)
d_qua <- d_mcmc %>%
  filter(grepl('^theta\\[\\d+\\]$', Parameter)) %>%
  group_by(Parameter) %>%
  summarize(`2.5%` = quantile(value, probs=.025),
            `50%`  = quantile(value, probs=.5),
            `97.5%`= quantile(value, probs=.975))

p <- ggplot() +
  geom_pointrange(data=d_qua, mapping=aes(x=forcats::fct_rev(Parameter), y=`50%`, ymin=`2.5%`, ymax=`97.5%`)) +
  coord_flip() +
  labs(x='Parameter', y='Value')
ggsave(p, file='fig1.png', dpi=300, w=4, h=3)
  • 6行目:Web上から8schools.stanを読み込んで文字列としています。RStanの公式ページの例題で使われているモデルファイルです。
  • 7行目:stan関数はmodel_code引数でモデルを書いた文字列も指定できます(本来はファイル名を直接指定できればよかったのですがよくわかりませんでした)。
  • 9行目:{ggmcmc}パッケージのggs関数でtidyなデータにしておきます。tidyなデータについては西原さんの記事を参照。{ggmcmc}パッケージに含まれる関数を使うとd_mcmcからいろいろな図が描けます。詳しくは『StanとRでベイズ統計モデリング (Wonderful R)』の4章に書きましたので読んでいただけるとうれしいです。

なお、d_mcmcは以下のようなデータフレームになります。

> d_mcmc
# A tibble: 72,000 × 4
   Iteration Chain Parameter   value
       <dbl> <int>    <fctr>   <dbl>
1          1     1        mu -1.3476
2          2     1        mu -0.9601
3          3     1        mu  7.0919
4          4     1        mu 15.0782
5          5     1        mu 20.0110
6          6     1        mu 20.4483
7          7     1        mu 13.0249
8          8     1        mu 11.8232
9          9     1        mu 15.9213
10        10     1        mu 17.9294
# ... with 71,990 more rows
  • 11行目:まずはモデルに含まれるtheta[数字]というパラメータだけ残しています。grepl関数でパラメータ名がマッチするか判定する際に正規表現を使う必要があります。ここが正規表現に慣れていない人は少し厳しいかもしれません。
  • 12~15行目:{dplyr}パッケージの典型的な使い方です。Parameter列ごとに要約量を算出します。列名が数字で始まる場合はバッククォートで囲む必要があります。1つずつ分位点を算出するのではなく、do関数で一行で算出する方法もあるのですが、分かりにくく、現状issueとして検討中のようです(ここここ)。
  • 18行目:ggplot2coord_flipすると下から上に向かってfactorが並びますので、{forcats}パッケージのfct_rev関数で逆順にしています。

* * *

次に以下の図を描いてみます。久保本11章の図に相当します。

library(rstan)
library(ggmcmc)
library(dplyr)

Y <- read.csv(url('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap12/input/data-kubo11a.txt'))$Y
I <- length(Y)
d <- data.frame(X=1:I, Y=Y)
data <- list(I=I, Y=Y)
model_code <- readr::read_file(url('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap12/model/model12-11.stan'))
fit <- stan(model_code=model_code, data=data, seed=1234)

d_mcmc <- ggs(fit)
d_qua <- d_mcmc %>%
  filter(grepl('^Y_mean\\[\\d+\\]$', Parameter)) %>%
  tidyr::separate(Parameter, into=c('Parameter', 'x'), sep='[\\[\\]]') %>%
  mutate(x=as.integer(x)) %>%
  group_by(Parameter, x) %>%
  summarize(`2.5%` = quantile(value, probs=.025),
            `10%`  = quantile(value, probs=.1),
            `50%`  = quantile(value, probs=.5),
            `90%`  = quantile(value, probs=.9),
            `97.5%`= quantile(value, probs=.975))

p <- ggplot() +
  geom_ribbon(data=d_qua, mapping=aes(x=x, ymin=`2.5%`, ymax=`97.5%`), alpha=1/6) +
  geom_ribbon(data=d_qua, mapping=aes(x=x, ymin=`10%`,  ymax=`90%`),   alpha=2/6) +
  geom_line(data=d_qua, mapping=aes(x=x, y=`50%`)) +
  geom_point(data=d, aes(x=X, y=Y), shape=1, size=2) +
  labs(x='i', y='Y[i]') +
  ylim(0, 22)
ggsave(p, file='fig2.png', dpi=300, w=4, h=3)
  • 15行目:{tidyr}パッケージのseparate関数を使って、Y_mean[20]Y_mean20という2つの列に分解しています。
  • 17行目:あとは集計の単位であるgroup_byの単位が場面によって多少変わるぐらいで、特に悩まずに色々な量が算出できます。

この記事ではMCMCChainごとに何かを算出することは取り上げませんでしたが、ggs関数で作ったd_mcmcChain列も含んでいますので自由自在です。

2次元以上のスポット検出を行う統計モデル

1次元の場合の変化点検出は以下の記事で扱いました。

変化点検出のポイントは状態空間モデルのシステムモデルに「ほとんどの場合は0の近くの値を生成するが,まれにとても大きな値を生成する」という性質を持つ分布を使うことでした。そこで、上記の例ではコーシー分布を使いました。しかし、コーシー分布はかなり厄介な分布なので逆関数法を利用したコーディングをしないとうまくサンプリングできません。そのため2次元以上の場合の同様のモデル(マルコフ場モデル)に拡張することができません。2次元の場合のマルコフ場モデルについては以下の記事で簡単に扱いました(少しコードが古いです)。詳しくは書籍「StanとRでベイズ統計モデリング」の12章を参照してください。

それでも2次元以上のスポット検出を行いたくて色々なモデルを試しました。コーシー分布+実装の工夫、線過程(ライン過程)、混合正規分布、ベルヌーイ分布と正規分布の混合分布、両側指数分布(ラプラス分布)を使ったモデルなどなど。しかし、少なくとも私の環境においてはこれらのモデルだと、推定がうまくいかない場合や、推定できてもスポットのようにシャープではなくボワッと境界があいまいに推定されてしまう場合が多かったです。その後も2年にわたり粘り続けて合計150個弱のモデルを試した結果、一つの解決策を見つけましたのでこの記事を書きます。

UBN分布

イデアは「ほとんどの場合は0の近くの値を生成するが,まれにとても大きな値を生成する」という性質を工夫して表現することです。そこで、状態空間モデルのシステムモデルに相当する部分に以下の分布を使いました。

 

 

この分布は [a,b]の範囲の値を生成する、切断正規分布に一様分布のゲタをはかせたような分布です。ここで、 a,bは一様分布の範囲を決めるパラメータ、 uは一様分布のゲタの高さを決めるパラメータ、 \mu,\sigma正規分布の平均と標準偏差 Cは正規化項です。この分布は一様分布と切断正規分布の混合分布なのですが、ベルヌーイ分布とポアソン分布の混合分布をZero-inflated Poisson分布と分かりやすく呼ぶのになぞらえて、この記事では上記の分布をUniformly Boosted Normal Distributionと呼び、以下では略してUBN分布と呼びます。思いつけばコロンブスの卵ですが、UBN分布はコーシー分布の利点である裾の長さを保ったまま、サンプリングの非効率さを解消します。例えば、 a=-1, b=9, u=0.1, \mu=3, \sigma=1の場合には Cはほぼ2となり、以下のような確率密度になります。

なお、Stanの計算に必要な対数尤度は以下のように式変形して表現することで高速化できます。

 

ここでlog1p_exp関数とnormal_lpdf関数は計算の高速化するためのStanの便利関数です。

2次元のスポット検出例 その1

30x30の大きさで中央部分にシグナルがあり、ノイズを加えたデータYを考えます。図にすると以下です。

上下左右のマスとつながるマルコフ場モデルを、UBN分布を利用して実装したStanコードは以下になります。

  • 5行目: 上記で説明したUBN分布の uを固定値Uで与えています。
  • 9行目: muにノイズがのってデータYになります。設定した下限と上限がUBN分布の a bにそれぞれ対応しています。
  • 17,20行目: 17行目がマルコフ場モデルの左右のつながり、20行目が上下のつながりに対応します。ここではUを固定値として与えていること、またs_muの値が高々1.5ぐらいと考えられることから、UBN分布の u/Cはほぼ一定の定数です。Stanでは対数確率の偏微分が重要であり、微分で消える定数項を無視して構いません。そのため u/Cを無視してコーディングしています。
  • 21行目: ノイズをのせている部分です。状態空間モデルの観測モデルに相当します。
  • 22,23行目: なくても構いませんが、変分ベイズが安定になることがあるのでs_mus_Yに弱情報事前分布を設定しています。

このモデルは局所最適解が多くMCMCの一種であるNUTSではうまく推定できませんが、自動変分ベイズだと高速に推定できるモデルとなっています。データを生成し、推定を実行するRコードは以下になります。

  • 4~9行目:データ生成部分です。
  • 12行目:ここではU0.1としました。少し値を変えて実行してみて結果が頑健かチェックするとよいでしょう。

変分ベイズが高速なため、計算時間はSurface Pro 3(core i5)でも10秒以下です。推定結果は以下です。以下の図はmuの乱数サンプルの中央値を2次元でプロットした図です。

以下の図はRowの真ん中(ここでは15)に固定してスライスし、横軸をColumn・縦軸をmuにしてプロットした図です。薄い灰帯は乱数サンプルから計算した95%ベイズ信頼区間、濃い灰帯は80%ベイズ信頼区間、青線は中央値、赤線はデータです。

シグナルのスポットを検出できているのがわかるかと思います。上記のRコードと結果はシグナルの大きさがノイズのSDの3倍となっている場合ですが、シグナルの大きさがSDの2.5倍となっている場合の結果は以下になります。徐々にカドから崩れていきます。

このシグナルが正方形の場合、少なくともシグナルがノイズのSDの2.3倍ぐらいないとスポットとして綺麗に検出されないようでした。

2次元のスポット検出例 その2

素数が300x600の場合も試してみました。データ例を図にすると以下です。

モデルは同じままで、実行するRコードもほぼ同じです(データ生成部分だけわずかに異なります)。推定に要する時間は5分ぐらいでした。シグナルがノイズのSDの3倍の場合の推定結果は以下になりました。図の凡例は上記の場合と同じです。丸いスポットの検出は比較的得意ですが、尖っている箇所はやはり難しいようです。

この場合、少なくともシグナルがノイズのSDの3倍ぐらいないとスポットとして綺麗に検出されないようでした。

試した範囲では、この方法によって地図上でも3次元でもスポット検出できそうです。Enjoy!

統計・R・Stan関連の本、用途別のオススメ10冊

年末年始向けに、比較的読みやすい本を中心にオススメします。

統計学入門

色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布とヒストグラム、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換、検定、散布図と箱ひげ図、回帰、相関あたりをRなどを使いながらシンプルに説明していく本があるといいと思うのですが、なかなかバランスのとれたいい本がありません。初歩の初歩しか説明してない、グラフが少ない、検定にページを割きすぎ、分厚い、ちょっと難しいなどの不満点があります。立ち読みして自分にあった本を選ぶのがいいと思います。ネットで検索して調べるのでもいいと思います。

統計学演習

統計の実力向上のためには、「紙と鉛筆で手を動かして計算する」という経験がどこかで必要になります。そこで、シグマと積分はすでに慣れ親しんで、統計の初歩は少し知っている人に以下の演習書をオススメします。

統計学演習

統計学演習

演習問題の難易度と内容が素晴らしいです。しかし、演習の前にある簡潔なまとめは説明がないので、復習として使うのがよいでしょう。

統計や検定を少し知っている人がRをはじめたい

「Rで楽しむ統計」がオススメです。例題がすごくいいです。また、重要であるのにもかかわらず一般的な本では説明を飛ばしているようなトピックを扱っていて勉強になります。

Rで楽しむ統計 (Wonderful R 1)

Rで楽しむ統計 (Wonderful R 1)

統計や検定を少し知っている人が統計モデリングをはじめたい

「データ解析のための統計モデリング入門」がオススメです。検定を乱用している人はぜひ読んでほしいです(乱用している自覚がある人は乱用しないかもですが)。

プログラミング少しできる人が統計モデリングをはじめたい

我田引水ですが「StanとRでベイズ統計モデリング」がオススメです。詳しくはこちらの記事を参照してください。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

Stanにもっと詳しくなりたい

Stanのマニュアルがオススメです。主要パートについては日本語訳もほとんど完成しており、以下のプロジェクトのReleasesタブに200ページ弱のpdfがあります。 github.com もともとの英語が難しく、訳に苦戦している箇所もあるので、プルリクエストお待ちしております。

因果推論や効果測定を詳しく知りたい

「岩波データサイエンス vol.3」がオススメです。統計を扱う上で非常に重要なのが因果関係ですが、この本を除いて読みやすい類書がほとんどありません。傾向スコアだけでなく、反実仮想を考えた回帰モデルやRCT(ランダム化比較試験)との対応がきちんと書いてあって勉強になりました。

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

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

解析で利用する立場からRをバランスよく詳しくなりたい

「改訂3版 R言語逆引きハンドブック」がオススメです。基本グラフィックスやファイル操作も詳しいです。見開きで読みやすいです。

改訂3版 R言語逆引きハンドブック

改訂3版 R言語逆引きハンドブック

統計に興味はないけど、プログラミング言語としてのRを知りたい

R言語徹底解説」がオススメです。エンジニア向けと思います。可視化などはあまり載っていません。

R言語徹底解説

R言語徹底解説

Rのデータ処理に詳しくなりたい

Rといえばデータ処理能力と作図能力がずば抜けています。データ処理については今のところ「Rデータ自由自在」がオススメです。

Rデータ自由自在

Rデータ自由自在

しかし、内容が少し古くなってきた感じがあります。現在は{dplyr}{tidyr}などのパッケージを使ってデータ処理するのが主流になりつつあります。2017年にWonderful Rのシリーズとして良い本が出るようなので待ちましょう。

Rの可視化(ggplot2)を使ってみたい

今のところ「Rグラフィックスクックブック」がオススメです。カラーページが多くびびります。

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集