モデル選択の実験 (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)”の続きを読む

モデル選択の周辺の話の整理

モデルを選択したり、変数を選択したり、というようなことに関係しそうなネタを簡単に整理してみた。

情報量基準

AIC / BIC / DIC / TIC のような。データのあてはまりのよさとモデルの複雑度を天秤に図るタイプのやつ。

たくさん種類があるのは確率モデルに関する仮定と汎化誤差の近似の仕方の違いによるものだと理解している。

検定

回帰係数が0である、という帰無仮説を検定することである変数が貢献しているかどうかを定量化するタイプのやつ。棄却されなければ「えい!」と変数を削ってしまう。

ささっと分析してデータの雰囲気を掴みたい時に使うことはある。

L1正則化

寄与度の小さな(ある閾値より小さな)係数をゼロにしてしまう、というような感じのやつ。

事前分布としてラプラス分布を使うことに相当。単純に寄与度が低い変数は消してしまえ!というノリなのだろうか?もっと深遠な背景があるのだろうか?勉強不足でよくわからない。

ベイズモデル選択

複数のモデルに事前分布を設定して、「モデルの事後分布」を計算するたぐいのもの。事後分布が求まったあとはMAPなものを選んでくるか、事後分布で平均をとってしまうか。

たとえばディリクレ混合過程。これはGMMのような混合モデルの混合数の事後分布を求めることができる。

PRMLの変分ベイズのところで出てきた関連度自動決定もこのタイプだと思っていいのだろうか。これも勉強不足により不明。

***

いろいろと抜けがあるとは思うが、とりあえずすぐに思いついたのはこれくらい。場合によっては追記します。

ハイブリッドモンテカルロの実験

相関の強い二変量正規分布に対してハイブリッドモンテカルロを使ってみた。上から順に、サンプリング結果、x1の自己相関、x2の自己相関。

自己相関ほぼ完全になし、という結果になった。ギブスサンプラーだとこうはいかない。ただし、

  • 微分方程式を解く時間のスケールが小さすぎると自己相関が出たので良い感じのスケールをちょっとだけ探索した。
  • ステップ幅を固定にしたら怪しげな自己相関の挙動がでた。
  • 計算時間でギブスサンプラーと比較してどちらが有利かは今回は検討してません。

という点は追記しておきます。

“ハイブリッドモンテカルロの実験”の続きを読む

対数線形モデルとエントロピー最大化の関係

昔の勉強ノートを引っ張りだしてくるシリーズ.

機械学習の対数線形モデルが最大エントロピー法とも呼ばれる,みたいな記述は頻繁に目にするし,統計力学のボルツマン分布の話とか考慮すれば,なんとなくそうなってそうな気はするけど,実際どうなの?というのを (たんに好奇心を満たすために) 調べてみた.実用上は何の意味ないと思う.

対数尤度関数に L1 正則化項を加えるタイプの目的関数を使った場合,もはやエントロピーは最大化されない,とかそういうわりとどうでもいいことがわかったりするかもしれない.

概要

「言語処理のための機械学習入門 (→ amazon) 」などに出てくるタイプの対数線形モデルの係数の最尤推定量が,エントロピーを「ある制約条件下」で最大化した場合のラグランジュ未定乗数に対応することを説明する (クロス表の対数線形モデルとはたぶん別物).

ただし,記号が煩雑になるのを避けるため,対数線形モデルとほぼ同一の構造を持ち,記号が煩雑でない条件付きロジットモデルがエントロピー最大化と等価であることを見る.

本文の最後に対数線形モデルと等価なエントロピー最大化問題を示す.多少ややこしくなるが,同じ方針で証明可能.
“対数線形モデルとエントロピー最大化の関係”の続きを読む

逆Wishart分布を図示してみる

正規分布のベイズ推定で共分散行列の(共役)事前分布として使われる逆Wishart分布(→ wikipedia:en)というものがありますが、あまりイメージしにくかったので図にしてパラメータの変化がどのような違いをもたらすか調べてみました。

といっても逆Wishart分布は n x n 行列を値にとる確率分布なので、2×2 行列を考えても3次元なので(※対称行列)そのままプロットするとわかりにくい。そこで ellipse パッケージを使って可視化してみます。
“逆Wishart分布を図示してみる”の続きを読む

事後分布に収束していく様子をグラフ化してみる(ポアソン/ガンマ分布)

ポアソン分布のパラメータの事前分布としてガンマ分布を仮定すると事後分布もガンマ分布になる(過去の記事 → 共役事前分布について)。

平均3のポアソン分布から100個サンプルをとってきて、オンライン学習的に次々にデータを食わせて事後分布がどのように変化していくかを図示してみた。事前分布としては \(\alpha=0.1, \beta=0.1\) のガンマ分布を使用。

データを増やしていくとなんとなくいいかんじに分布が収束していくことがわかる。こういう図をうまく見せるのって難しい。以下、Rスクリプト。

library(RColorBrewer)
x <- rpois(100,3)
lam <- seq(0,8,0.01)
hc=brewer.pal(11,"RdYlGn")
plot(lam,dgamma(lam,0.1,0.1),type="n",
     ylim=c(0,3),xlab="Lambda",ylab="Density",col="gray")
lines(lam,dgamma(lam,0.1,0.1),col="gray",lw=2)
i <- 1
for( n in c(seq(1,10,2),seq(20,100,20)) ){
  lines(lam,dgamma(lam,0.1+sum(x[1:n]),0.1+n),
        col=hc[i],lw=2)
  i=i+1
}
title("Posterior of Poisson distribution parameter with Gamma(0.1,0.1) prior")
lines(c(3,3),c(0,5),lw=2,lt=2,col="darkblue")
legend("topright",
       legend=c("Prior (Gamma(0.1,0.1))",
                paste("n =",c(seq(1,10,2),seq(20,100,20)))),
       lw=2,col=c("gray",hc))
dev2bitmap(file="PoissonPosterior.png",
           gaa=4,taa=4,width=9,height=6)

正定値行列と共分散とフィッシャー情報量(と機械学習との関連について少々)


確率統計ででてくる基本的な2つの行列(共分散行列、フィッシャー行列)の半正定値性はほとんど同じ方法で証明できる、ということについて。最後の方は機械学習と関連したすこしマニアックな話題。

関連1 → Cramer-Rao’s theorem
関連2 → 正定値行列と共分散(に関する駄文)

共分散行列の半正定値性

そもそも、共分散行列が半正定値になる理由は、任意のベクトル \(u\) に対して
\begin{align}
u^T\rm{Cov}[\boldsymbol X]u
&=\int u^T(\boldsymbol x-\bar{\boldsymbol x})
(\boldsymbol x-\bar{\boldsymbol x})^TudF(x)\\
&=\int |u\cdot(\boldsymbol x-\bar{\boldsymbol x})|^2dF(x)\ge0
\end{align}
となることから明らか。

フィッシャー情報量の半正定値性

同様に、フィッシャー情報量が半正定値行列になる理由は \(l(x|\theta)=\log f(x|\theta)\) として
\begin{align}
u^TI(\theta)u
&=u^TE[\nabla_\theta l(X|\theta)\cdot\nabla_\theta^Tl(X|\theta)]u\\
&=\int |u\cdot\nabla_\theta l(x|\theta)|^2dF(x|\theta)\ge0
\end{align}
より明らか。よって逆行列が存在すれば(つまりフィッシャー行列が退化していなければ) \(I(\theta)^{-1}\) も正定値。

さらに、任意の行列 \(A\) に対して \(A^TI(\theta)^{-1}A\) も半正定値(\(A\) が退化することもあるので「半」がくっついてくる)。これも上と全く同じ論法で終了。なのでクラメール・ラオのエントリーで出てきた次の式
\begin{align}
\mathrm{Cov}_{\boldsymbol{\theta}}
\left(\boldsymbol{W}(\boldsymbol X)\right)\geq \frac {\partial \boldsymbol{\tau}
\left(\boldsymbol{\theta}\right)} {\partial \boldsymbol{\theta}}
[I\left(\boldsymbol{\theta}\right)]^{-1}
\left( \frac {\partial
\boldsymbol{\tau}\left(\boldsymbol{\theta}\right)}
{\partial \boldsymbol{\theta}}\right)^T
\end{align}
の右辺も左辺も半正定値になることがわかる。

機械学習とフィッシャー情報量

このくらいまで書くとフィッシャー情報行列が退化する、というワクワクする状況についても書きたくなってくる。フィッシャー情報量 \(I(\theta)\) がモデルの真のパラメータのところで退化するケース。ゼロ固有値を持ち、その逆行列は存在しない(形式的には固有値ゼロで割るので逆行列が無限大に発散する)。こういう状況は機械学習でしばしば生じる。

(誤差項に正規分布を仮定して確率モデルとして扱った場合の)多層パーセプトロンとかがこれに該当する。他にはHMMとか隠れノードのあるベイジアンネットとか。hidden state がある場合に退化することが多いっぽい。ちなみに SVM はそもそも確率モデルでないので対象外。

フィッシャー行列が退化すると何がまずいのか、というと、多層パーセプトロンはバックプロパゲーションするとローカルミニマムに落ちるからそもそも最尤推定は難しいのだけれども、直感的に言うと、たとえ最尤推定できたとしても、クラメール・ラオの下限が発散するので、最尤推定量のバラつきも無限大になる、ということを暗に意味している。

これは(悪い意味で)指数型分布族とは著しく異なる性質だ。指数型分布族的なものを正則モデル(regular model)というのに対して、特異モデル(singular model)というらしい。こいつらに対しては AIC などのモデル選択も理論的にはうまくいかない(実用的にはそれなりにうまくいくという話は周りの人たちからはよく聞く)。

そういうわけで、点推定はうまくいかないので、ベイズ推定をすることになる。ベイズ推定なら特異モデルに対しても漸近的な性質を導くことができる。くわしいことは特異モデルの機械学習のベイズ推定の理論の本 “Algebraic Geometry and Statistical Learning Theory” に書いてある(難易度:個人的な感想としては激高)。

Cramer-Rao’s theorem

Cramer-Rao’s theorem について.簡単に言うと確率モデルのパラメータを推定する際の精度の理論的限界を定める定理(もちろんこの定理を適用するためにはいくつかの仮定=前提条件があるのですが).1パラメータのケースは 『Casella & Berger 本』 に書いてあるけれど,パラメータが1次元より大きなときについては書いてなかったので考えてみた.

(自分のメモをブログに転載しただけなのであまりフレンドリーな書き方ではないかもしれませんが,いつか誰かの役に立つかもしれないので残しておきます.昔に書いた関連エントリはコチラ

考えている状況とか、記号の定義とか

\(\boldsymbol X\) を \(f(\boldsymbol X|\boldsymbol\theta)\) に従う確率変数とする.\(\boldsymbol{\theta}\) は分布の母数.\(\boldsymbol\tau(\boldsymbol\theta)\) は母数を一対一変換するような微分可能な関数.\(\boldsymbol W(\boldsymbol X)\) は \(\boldsymbol\tau(\boldsymbol\theta)\) の不偏推定量,つまり
\begin{align}
E_{X|\theta}(\boldsymbol W(\boldsymbol X))&=\boldsymbol\tau(\boldsymbol\theta)
\end{align}
が成立するとする.ようするに \(\boldsymbol W\) は観測されたデータをモデルパラメータに変換する写像(アルゴリズム)と考えればよい.\(\boldsymbol W(\boldsymbol X)\) は確率変数なので,バラつきがある.つまり推定されたモデルパラメータ自身が確率変数になっている.このバラつきが大きいと精度の良い推定は望むべくも無い.そこでパラメータのバラつきを評価するの理論的な限界(下限)を定めるための定理が Cramer-Rao.推定量といってもいろいろあって,「データがどんな値をとってもパラメータは0!」と決め打ちしてしまうアルゴリズムも考えられる.こういう推定量はバラつきはゼロだけど精度がいいとは言えない.不偏推定量という制約条件を課しているのはこのようなおかしな推定量を排除して「筋の良いもの」だけを相手にするための仕掛けと考えればいい.

この定理があってはじめて,「よい推定方法とは何か?」という問いに答えられるようになる.すなわちその答えは「不偏でかつバラつきの理論的限界を達成するような推定量」こそが最良である,となる.

定理;一般化されたクラメール・ラオの下限

\(I(\boldsymbol\theta)\) をフィッシャー情報量とするとき,
\begin{align}
\mathrm{Cov}_{\boldsymbol{\theta}}
\left(\boldsymbol{W}(\boldsymbol X)\right)\geq \frac {\partial \boldsymbol{\tau}
\left(\boldsymbol{\theta}\right)} {\partial \boldsymbol{\theta}}
[I\left(\boldsymbol{\theta}\right)]^{-1}
\left( \frac {\partial
\boldsymbol{\tau}\left(\boldsymbol{\theta}\right)}
{\partial \boldsymbol{\theta}}\right)^T
\end{align}
が成り立つ.ただし,\(A\ge B\) は \(A-B\) が半正定値行列であることを意味している.

以下、証明っぽいもの。 “Cramer-Rao’s theorem”の続きを読む