読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

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

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

Stan 理論

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

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

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

UBN分布

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

 

 

この分布は [a,b]の範囲の値を生成する、正規分布に一様分布のゲタをはかせたような分布です。ここで、 a,bは一様分布の範囲を決めるパラメータ、 uは一様分布のゲタの高さを決めるパラメータ、 \mu,\sigma正規分布の平均と標準偏差 Cは正規化項です。この記事ではこの分布を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 Stan

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

統計学入門

色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布とヒストグラム、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換、検定、散布図と箱ひげ図、回帰、相関あたりを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によるグラフ作成のレシピ集

情報量規準LOOCVとWAICの比較

R Stan

この記事はStan Advent Calendar 2016およびR Advent Calendar 2016の12月7日の記事です。StanコードとRコードは記事の最後にあります。

背景は以下です。

  • [1] Aki Vehtari, Andrew Gelman, Jonah Gabry (2015). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. arXiv:1507.04544. (url)
  • [2] 渡辺澄夫. 広く使える情報量規準(WAIC)の続き (注4)【WAICとクロスバリデーションの違いについて】 (url)
  • [3] Sumio Watanabe. Comparison of PSIS Cross Validation with WAIC. (url)

leave-one-outクロスバリデーション(以下LOOCV)およびWAICは予測のよさをベースにしたモデル選択に用いられる情報量規準であり、ともに汎化誤差(Generalization Error、以下GE)の近似です。それにもかかわらず[1]では、本来性能評価では必須と思われる汎化誤差との比較がありません。実データ(真のモデルが未知の状況)で用いているためと思いますが、これではいけないように思います。この記事では僕が日常的に使用するような5つの基本的なモデルを使い、真のモデルが既知の状況でGE・LOOCV・WAICの性能比較を行いました。

具体的には[2]に以下のようなコメントがあります。なお[1]ではPareto Smoothed Importance SamplingでLOOCVを算出しており、PSISCVとも呼ばれるようです。

  • (0) まず同じデータに対してマルコフ連鎖モンテカルロ法を何度も行ったときの値の揺らぎを調べてみましょう.WAICの分散はISCVおよびPSISCVの分散よりも小さくなります。つまり、WAICはISCVおよびPSISCVよりもマルコフ連鎖揺らぎに対して強いということができます.
  • (2) しかしながら,CVもWAICも汎化誤差を推定することが本来の目的です.CVとWAICの両方の厳密値が計算できたとして,(つまりMCMC法で無限にサンプルが取れたとき),どちらの方が汎化誤差の推定として優れているのでしょうか。(中略) 我々の実験では,GEを汎化誤差とするとき,ほぼ,いつでも E[|PSISCV-GE|] > E[|WAIC-GE|] が成り立つのですが・・・。このページをご覧の皆様にはぜひ実験してみていただければと思います。なお、E[ ]は学習用データのでかたについての平均を表しています。

そこでこれら2点について検証してみました。先に「(2)のコメント」について説明します。

GE・LOOCV・WAICの比較

使用した5つのモデルとシミュレーションの手順を説明します。

重回帰

真のモデルは以下です。

  Y \sim Normal(1.3 - 3.1 X_1 + 0.7 X_2, 2.5)

あてはめたモデルは以下です。

  Y \sim Normal(b_1 + b_2 X_1 + b_3 X_2, \sigma)

データ点の数Nについては10,30,100,300を試しました。例としてN = 10の場合を説明します。まず乱数でデータX(すなわち X_1, X_2)を生成します。次にそのXの値を使って1000通りのYを生成します(学習用データのでかたの平均をとるため)。各YについてStanでiter=11000, warmup=1000, chains=4で実行して合計40000個のMCMCサンプルを得てGE・LOOCV・WAICと「LOOCV - GE」と「WAIC - GE」を求めました。その後、Nごとに「LOOCV - GE」と「WAIC - GE」のboxplotを描きました。

ロジスティック回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

  Y \sim Bernoulli(inv\_logit(0.8 - 1.1 X_1 + 0.1 X_2))

あてはめたモデルは以下です。

  Y \sim Bernoulli(inv\_logit(b_1 + b_2 X_1 + b_3 X_2))

  b_1,b_2,b_3 \sim Student\_t(4,0,1)

非線形回帰 ミカエリス・メンテン型

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

  Y \sim Normal(10.0 X / (2.0 + X), 0.8)

あてはめたモデルは以下です。

  Y \sim Normal(m X / (k + X), \sigma)

  k \sim Uniform(0, 12)

  m \sim Uniform(0, 20)

真のモデルが含まれない場合

あてはめたモデルが以下の場合も試しました。

  Y \sim Normal(a + b X, \sigma)

ノイズがt分布に従う重回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

  Y \sim Student\_t(4, 1.3 - 3.1 X_1 + 0.7 X_2, 2.5)

あてはめたモデルは以下です。

  Y \sim Student\_t(4, b_1 + b_2 X_1 + b_3 X_2, \sigma)

階層モデル

手順は重回帰の場合とおおよそ同じですが、データ点の数とモデルが異なります。グループの数を10に固定し、データ点の数Nについては20,50,130,400を試しました(それぞれ各グループで2,5,13,40人)。また真のモデルは以下です。

  \mu\left[g\right] \sim Normal(0, 10)   g = 1,2,\dots,10

  Y\left[n\right] \sim Normal(\mu\left[g\left[n\right]\right], 2.0)   n = 1,2,\dots,N

あてはめたモデルは以下です。

  \mu\left[g\right] \sim Normal(\mu_0, \sigma_0)   g = 1,2,\dots,10

  Y\left[n\right] \sim Normal(\mu\left[g\left[n\right]\right], \sigma)   n = 1,2,\dots,N

  \mu_0,\sigma_0 \sim Student\_t(4, 0, 10)

階層モデルにおいては何を予測したいのか(どういう状況の汎化誤差を知りたいのか)を注意深く考える必要があります。以下の記事を参照。

statmodeling.hatenablog.com

ここでは以下の2つの場合を計算しました。

  • 既存の各グループに、新しく1人ずつ加わる場合
  • 別の新しいクラスができて、新しく1人が加わる場合

結果

重回帰

大きな差はありませんでした。Nが小さい場合にWAICの方がわずかにGEに近くなる傾向があるようです。

ロジスティック回帰

大きな差はありませんでした。Nが小さい場合にLOOCVの方がわずかにGEに近くなる傾向があるようです。

非線形回帰 ミカエリス・メンテン型

大きな差はありませんでした。Nが小さい場合にWAICの方がわずかにGEに近くなる傾向があるようです。

  • 真のモデルが含まれない場合

この場合はNを増やしても汎化誤差引くエントロピー(=予測分布と真の分布のKL情報量)は0に近づかず、0.89ほどで下げ止まります。そして、Nが小さい場合にWAICの方がGEに近くなる傾向があるようです。

ノイズがt分布に従う重回帰

Nが小さい場合にWAICの方がGEに近くなる傾向があるようです。なおboxからはみ出るoutlierの値が大きく、そのままプロットすると見づらくなるので、図の縦軸を制限しました。

階層モデル

  • 既存の各グループに、新しく1人ずつ加わる場合

Nが小さい場合にWAICの方がGEに近くなる傾向があるようです。

  • 別の新しいクラスができて、新しく1人が加わる場合

グループ数が少ない場合、グループあたりの人数を増やしてもLOOCV - GEおよびWAIC - GEの中央値は0に近づいていきません。グループ数が少ないと、グループあたりの人数を増やしてもグループを生成する平均と標準偏差のパラメータは精度よく求められないことを反映しているのだと思います。全体的にWAICの方がわずかにGEに近くなる傾向があるようです。

階層モデル その2 2016.12.14追記

伊庭先生から以下のような要望がありました。

そこでモデルは同じものを用いて、グループあたりの人数を固定し(2,5,13のうちいずれか)、グループの数が10,30,100の各場合で実行してみました。

  • 既存の各グループに、新しく1人ずつ加わる場合

グループあたりの人数が少ない場合、グループ数を増やすとLOOCV - GEおよびWAIC - GEのバラツキは少なくなるもの、それらの中央値は0に近づいていきません。グループ数を増やしても各グループには2人しかいないため、グループごとの予測はあたらないままということを反映しているのだと思います。グループあたりの人数が小さい場合にWAICの方がGEに近くなる傾向があるようです。

  • 別の新しいクラスができて、新しく1人が加わる場合

Gが小さい場合にWAICの方がわずかにGEに近くなる傾向があるようです。

* * *

5つのモデルを通して見ると、Nが小さい場合、すなわち1サンプルの重みが大きい場合にはWAICのほうがLOOCVよりもよい汎化誤差の近似になっているようです。[2]の「PSISCVとWAIC:実験例 追加」によると、影響の大きい(重みの大きい)サンプルが存在する場合にも同じような結果になるようです。以下の伊庭先生のツイートはこのような状況を指していると思われます。

WAICはLOOCVよりMCMCの揺らぎに強いか?

次に「(0)のコメント」について検証しました。前述の5つのモデルに対し、N10または100とし、各Nについてデータを5通り生成しました。さらに各データセットに対して、MCMCのシードを1000通り試し、1000個のLOOCVとWAICを求め、それぞれの分散と平均と変動係数を求めました。

結果

いずれの場合についても実用上の差が認められませんでした。N100の場合より10の場合の方がMCMCの揺らぎがありますが、それでも変動係数にして高々1%程度でした。なお、この結果はMCMCサンプルを求めるアルゴリズムの違いやMCMCサンプルの数にも依存すると思います。

まとめ

  • Nが大きい場合には、WAICとLOOCVにほぼ差がないと言えるでしょう。
  • Nが小さい場合には、WAICの方が汎化誤差のよい近似になっていると言えるでしょう。
  • さらに理論的な美しさや計算速度も含めて総合的に判断すると、WAICに軍配が上がると思います。

おまけ

データ Xとは異なる xにおける yの予測分布 y_{pred}(y|x)を考えたい場合があると思います。その場合は、一般によい情報量規準があるか未解明で、研究対象として興味深いようです(渡辺澄夫先生(私信))。例えば単回帰の場合には、真のモデルがあてはめたモデルに含まれており、かつモデルに依存する量を使うと情報量規準に準ずる量を構築できるようです。興味深いです。


ソースコード

「(2)のコメント」を検証するための、重回帰の場合と階層モデルの場合のStanコードとRコードを公開します。

重回帰

Stanコード

重回帰のモデルの部分は「StanとRでベイズ統計モデリング」の9.2.1項と同じです。異なる点はgenerated quantitiesブロックでGE・LOOCV・WAICを算出するためにlog_liky_predを算出している点です。

Rコード

例としてN100の場合を載せます。理解しやすさのため、並列化していないコードにしてありますが、実際には色々並列化して計算しています。

  • 16~17行目: 学習用データの出かたの平均をとりたいので、データYを乱数で生成しています。
  • 22~31行目: データごとのエントロピーと汎化誤差を計算しています。以前の記事参照。このコードのように十分なMCMCサンプルから予測分布の近似関数を求めて汎化誤差を算出する方法のほか、直接MCMCサンプルを使って求める方法もあるかと思います(渡辺先生はそうしています)。
  • 23~24行目: 予測分布は滑らかだと仮定し、40000個のMCMCサンプルから予測分布の密度関数を計算しています。Rのデフォルトのdensity関数よりも{KernSmooth}パッケージのbkde関数の方が優秀っぽいのでこちらを使っています(参考pdf)。
  • 25行目: 真の分布です。
  • 26行目: エントロピーの計算の際に積分される関数です。
  • 27行目: 汎化誤差の計算の際に積分される関数です。予測分布の密度推定した結果をapproxfunで関数に変換している関係上、f_predがごくまれに絶対値の小さな0以下の値を返します。それを避けるためにifelse関数をかませてあります。
  • 28~29行目: それぞれエントロピーと汎化誤差を計算しています。ここではf_true正規分布なので、-6SD~+6SDまで積分すれば十分よい近似となります。Rのデフォルトのintegrate関数はちょっと賢くてadaptiveに積分しているので計算は早いのですが、まれに不安定で計算ができない場合があります。そのため、少し遅いですが安定な数値積分の手法であるRombergの方法を使っています({pracma}パッケージのromberg関数または{Rmpfr}パッケージのintegrateR関数を使うことができます)。
  • 34行目: 汎化誤差(GE)のサンプルに関する平均を求めています。学習用データと同じNXの値を持つ新しいデータセットに対して予測を行い、1サンプルあたりの汎化誤差を求めていることに相当します。
  • 35~36行目: Stanのチームが開発している{loo}パッケージを用いてlooicwaicを求めています。彼らの情報量規準のスケールはAICやDICとあわせて 2N倍となっているので、1サンプルあたりにスケールをあわせる意味で2*Nで割っています。

階層モデル「既存の各グループに、新しく1人ずつ加わる場合」

階層モデルの場合のWAICの詳しい解説は以前の記事を参照してください。

Stanコード

Rコード

重回帰の場合と似ています。

  • 23~32行目: グループgに1人加わった場合の汎化誤差error_by_groupを求めています。重回帰の場合にサンプルごと(nごと)だったのが、グループごと(gごと)にインデックスが変わっただけで、内容は変わっていません。
  • 35行目: 「既存の各グループに、新しく1人ずつ加わる場合」なので、各グループの汎化誤差の和になります。
  • 36~37行目: 各グループごとにLOOCVまたはWAICを求めて和をとっています。

階層モデル「別の新しいクラスができて、新しく1人が加わる場合」

Stanコード

Rコード

  • 25~26行目: この場合は真の分布が積分を含んでいるので少し複雑です。

残りは「既存の各グループに、新しく1人ずつ加わる場合」とほぼ同じです。

なお、この記事は以下のツイートによってやってみようかなと思いました。