モデル選択の実験 (BIC を追加)

前回の記事 では AIC と AICc を比較した。今回はそれに BIC を追加してみた。BICはあまり使ったことがなかったが、個人的には結構おどろきの結果が得られた。

BIC は以下で定義される。n はデータ数、k はモデルのパラメータ数。

\begin{align}
BIC=-2\times\{\text{Maximum Likelihood}\}+(\log n)\times k
\end{align}

実際のデータ分析では当然、n は固定なので AIC とのちがいは k の前の係数が 2 という定数か、 \(\log n\) という定数か、の違いがあるが、これって同じようなもんでしょ、と思って BIC はわりとノーマークだったけど、今回実験してみて考えを改めることなった。
“モデル選択の実験 (BIC を追加)”の続きを読む

モデル選択の実験 (AIC vs AICc / R の AICcmodavg パッケージ)

前回の「モデル選択の実験 (AIC vs LOOCV)」の続きです。

小標本の場合は、AIC じゃなくて AICc を使うといいよとのことなので、今回は、前回同様の方法で AIC と AICc を比較してみた。

真のモデルやその他のモデルの設定などは前回と全く同様。図の見方も前回と同様です。LOOCV の結果も並べてみたかったが計算量の関係で断念。

いろんな意味で手抜き気味です。あしからず。

AICc

「正規ノイズの線形モデル」のケースでは以下で定義される情報量基準。

\begin{align}
AICc = AIC + \frac{2k(k+1)}{n-k-1}
\end{align}
k : パラメータ数、n : データ数。

GLM のケースではこの定式化は使えないらしい[1]。 そんなこんなで GLM の場合は R の AICcmodavg パッケージをつかう。

glm.fit <- glm(...)
AICc(glm.fit)

とすることで AICc の値を良きに計らって計算してくれる。
“モデル選択の実験 (AIC vs AICc / R の AICcmodavg パッケージ)”の続きを読む

  1. [1] でも実は間違えてこの定式化でやってみたけどそれなりに良い結果 (AIC よりもよい結果) が得られた。

モデル選択の実験 (AIC vs LOOCV)

最近読んだ「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (→ amazon)」に感化されて、モデル選択基準である、AIC、LOOCV (Leave One-Out Cross Validation) について実験してみた(参考:モデル選択の周辺の話の整理)。

LOOCV と AIC はある意味で両極端だ。LOOCV は利用するのに必要な仮定が少なく、直感的にも何をやっているのかわかりやすい(ただし計算量は非常に大きい)。一方で AIC は評価式の簡単さ(計算量の少なさ)とは裏腹に、背後には正則条件だとかネストだとかの直感的とは言いがたい仮定をしたりする。

まあ AIC とか LOOCV とかが常に正しいモデルを選択するわけじゃない、というのは直感的にはわかるけど、実際に試してみたことはなかったので、真のモデルを知っているという神の視点からみて、こいつらがどんな挙動を示すのかをシミュレーションで確かめてみた。

“モデル選択の実験 (AIC vs LOOCV)”の続きを読む

指数型分布族のメモ

モチベーション

推定や検定の理論を扱うときに、正規分布の場合~、ポアソン分布の場合~、みたいに分布ごとに性質を調べてたんじゃやってられない、ということで、さまざまな密度関数をある種のクラスとしてまとめて考えることで、一般論を展開したい。そのための枠組みの一つが指数型分布族。

実際、指数型分布族に該当する分布は推定や検定などで「極めていいかんじの」性質を持っている。一方で実用上よく使う分布の多くが、指数型分布族というクラスを「継承」している(ほとんどが、と言ってしまってもいいかもしれないが、機械学習系の確率モデルには当てはまらないことが多い)。

要するに、指数型分布族の理論について知っていれば、実用上よく使う分布たちが極めてよい性質を持っていることが統一的に理解できる。ちなみに、GLM(一般化線形モデル)も基本的には指数型分布族前提で定式化されている。

適当に書いて見ましたが、こんなところでしょうか。。。

定義

n次元の確率変数 \(X\) が d-パラメータの指数型分布族に従うとは、あるパラメトリゼーション \(\boldsymbol{\theta}\) と d 個の関数 \(T_i:\mathbb{R^n}\to\mathbb{R}\hspace{5pt}(i=1\cdots d)\) が存在して、\(X\) の密度関数が
\begin{align*}
f(x|\boldsymbol{\theta}):=h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)
\end{align*}
と書けること。\(h, A\) はそれぞれスカラー値をとる関数。\(T_i(\cdot)\) は十分統計量と呼ばれる。

指数型分布族の分布の例

正規分布、ポアソン分布、指数分布、ガンマ分布、ディリクレ分布、負の二項分布などなど。数え出したらきりがない。

指数型分布族でない分布の例

一様分布、コーシー分布など。混合ガウス分布とか隠れマルコフモデルとか、その手のいわゆる「特異モデル」は該当しない。

Natural exponential family

分布族を
\begin{align*}
f(x|\boldsymbol{\theta}):=h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot x-A(\boldsymbol\theta)\big)
\end{align*}
に限定したものが natural exponential family。

平均と共分散の計算

指数分布族に従う d 次元確率変数 \(T(X)=(T_1(X),\cdots,T_d(X))\) の平均と共分散は
\begin{align*}
E(T_i(X))&=\frac{\partial A}{\partial\theta_i}\\
Cov(T_i(X),T_j(X))&=\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}
\end{align*}
となる。なぜなら、正規化条件
\begin{align*}
1=\int h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx
\end{align*}
の両辺を \(\theta_i\) で微分すれば
\begin{align*}
0=E_\theta(T_i(X))-\frac{\partial A}{\partial\theta_i}\cdot1
\end{align*}
となり、もう一度、今度は \(\theta_j\) で微分すれば
\begin{align*}
0=&\int T_i(x)T_j(x)\cdot h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx\\
&-\frac{\partial A}{\partial \theta_j}\int T_i(x)\cdot h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx
-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}\\
=&E(T_i(X)T_j(X))-E(T_i(X))E(T_j(X))-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}\\
=&Cov(T_i(X),T_j(X))-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}
\end{align*}
ということでOK(積分と極限の交換の条件とかは無視してしまいましたが)。

JAGSを使ってギブスサンプリングを試してみた (1)

前から気になってた JAGS を試してみました。Just Another Gibbs Sampler を略して JAGS。Bugs ファミリーの中では現在、最もアクティブに開発がすすめられているようです。

ここでは、とにかく何も知らない状態から R+JAGS で簡単なモデル(ガウスノイズな線形回帰)の推定までできるように最低限のことをメモしました。

準備

JAGS本体のインストール

JAGS は R とは分離された独立なソフトウェアなので R の CRAN 経由で一発インストールというわけにはいきません。ですが、Ubuntu の場合は apt で一発です。他のOSは試していないので、インストールしていない人は別途調べてください。

sudo apt-get install jags

rjagsのインストール

JAGSをRから呼び出すためには rjags というパッケージをインストールする必要があります。apt 使うなら

sudo apt-get install r-cran-rjags

もしくはRのプロンプトで

install.pacakges("rjags")

とすればOK。

テスト

テスト問題として

\begin{align*}
y_i&\sim \text{Normal}(\mu_i,\sigma^2)\\
\mu_i&=\alpha+\beta x_i
\end{align*}
というデータの発生メカニズムを仮定します。これに対応した次のような「モデルファイル」をつくって、適当に名前をつけておきます(lm.jagsとする)。

model {
  for(i in 1:N){
    mu[i] <- alpha + beta*x[i]
    y[i] ~ dnorm(mu[i],1/(tau^2))
  }
  alpha ~ dnorm(0,0.0001)
  beta ~ dnorm(0,0.0001)
  tau ~ dgamma(0.1,0.1)
}

見た感じ、Rのコードにかなり似ているけど、実際には細かい点で異なります。dnormが正規分布なのは同じだけど、第一引数は存在しないし、JAGSでは分散ではなくて精度(precision;分散の逆数)を指定することになっているので注意。Rを知っていればなんとなくわかってしまう構文は素晴らしいですが、上で挙げたような注意の他に、分布の名前自体がちがっていたりすることもあるのでハマる可能性もあります。ちなみに、for ブロックの後ろの3行は未知パラメータの事前分布を設定しています。

データの生成

\(\alpha=2, \beta=3\) としてデータを発生させます。

# determine the number of samples
N <- 50
# data generation
x <- rnorm(N) 
y <- 2 + 3*x + rnorm(N,0.1)

JAGSをRから呼び出す

モデルファイルを保存したら、Rで以下のスクリプトを走らせます。すると下のようなグラフが表示されるはずです。coda[1]でプロットさせているのでMCMCglmmなど他のRのベイズパッケージを使った時と同じ図になります。

library(rjags)
library(coda)

# compiling and initializing
lm.mcmc <- jags.model("lm.jags",
                       data=list('x'=x, 'y'=y, 'N'=N),
                       n.chains=4,  # the number of parallel chains
                       n.adapt=100)

# mcmc sampling
update(lm.mcmc,1000) # burn-in
posterior <- coda.samples(lm.mcmc, c('alpha','beta'), 2000)
# plotting the result
plot(posterior) # plot.coda is called


こんな感じで推定出来ました。すこしスクリプトについて書いておくと

  • jags.model 関数の第一引数では先程のモデルファイルを指定する。
  • 第二引数 data では JAGS に渡したいRオブジェクトをリスト形式で与える。クオーテーションする必要があるので注意。
  • n.chains で並列パスの数を指定
  • update で burn-in している
  • coda.samples で指定した数のサンプルを coda 形式で取得する。

まとめといろいろと蛇足

階層ベイズでも何でもないつまらない例を使って、ほんとに最低限のJAGSの使い方について書きました。たぶん続編があるはずです。次は混合ポアソン回帰でも試してみようかな。

わりとどうでもいいことですが、’J’AGS という名前(←JRubyみたいな意味で。)とクロスプラットフォームを謳っていることから Java で書かれていると思っていましたが、c++で実装されているようです。この勘違いのせいでずいぶん長い間敬遠してきた、という経緯があります。

ちなみに、これまで自分でMCMC実装するときは

  • cで書く
  • Rでものすごく遅いやつを書く(ちかいうちにこのブログにも登場予定)
  • MCMCglmmみたいなRのパッケージを使う

という選択肢の中から選んでいたのですが、これにJAGSを追加して比較するとこんな感じでしょうか?(いろいろと異論はあるかもしれませんが)

自由度 実行スピード 実装の手軽さ
c言語 ×
R ×
R-packages × ○ (?)
JAGS ○ (?)

おしまい。

  1. [1] RでMCMCの結果を可視化するときのほぼ事実上の標準形式

[R] MCMCglmm でポアソン回帰


R ネタです。一般化線形モデル(Generalized linear mixed model : GLMM)をベイズ推定する MCMCglmm というパッケージを使ってポアソン回帰をするためのメモです。もしかすると MCMCglmm の入門的、あるいはチュートリアル的なものになっている、かもしれません。 “[R] MCMCglmm でポアソン回帰”の続きを読む

分布の混合による負の二項分布の生成

ポアソン分布 (poisson distribution) は平均値 \(\lambda>0\) を唯一のパラメータとする離散確率分布ですが、これをガンマ分布 (gamma distribution) で混合してみることを考えます。つまりポアソン分布を \(p(x|\lambda)\)、ガンマ分布を \(q(\lambda|a,b)\) として、\begin{align*}r(x|a,b):=\int_0^\infty p(x|\lambda) q(\lambda|a,b)d\lambda\end{align*}を計算すると何が起こるのか、という話です。
“分布の混合による負の二項分布の生成”の続きを読む