MCMC

logsumexpを使って乱数生成 + X^2検定

空港で暇にしていたので書いてみる(出発時間が遅れた)。takanori-i君がベイジアンHMMを作っているらしく、相談に乗る。確率の積が入ってきて、数値計算で死んでしまうときがあるとのことだったので、logsumexpについて教える。logsumexpについては高村本が分…

マルコフブランケット

HMMのGibbs Samplingを考えるとき、マルコフブランケットの変数のみを考えればよいのだがなんでだっけとかアホなことを考えたりしたのでメモ。なんでマルコフブランケットだけ考えればいいかについてはパターン認識と機械学習 下 - ベイズ理論による統計的予…

とあるモデルのMCMC

飽きもせずにGibbs Samplingとかばっかりやってますが、久しぶりにはまった。離散確率分布とかからサンプリングするのにrandとかを使ってたんだけど、これがとてもとてもとてもいけなかった。CとかC++のrandは線形合同法で実装されているとかで周期性が問題…

SIR(重点サンプリング)を簡単な例で

木曜だと思っていたゼミが明日だということにさっき気がついてあたふたと準備をしています。。。担当している箇所はパターン認識と機械学習 下 - ベイズ理論による統計的予測の11.1.4の重点サンプリングと11.1.5のSIRです。重点サンプリングのところは去年研…

マルコフ連鎖の収束判定に関連するトピック

以前までのMCMCのシミュレーションでは、不変分布に収束したかと見なせるかは、時系列plotをして、変な動きをしていなければ収束しているんじゃね?という大雑把なものだった。ということで、今回はマルコフ連鎖がどういう状況ならば収束していると判定してい…

俺的MCMCまとめ

12月くらいからMCMCの勉強しだして、いくつか代表的なアルゴリズムによるサンプリングをやったのでまとめておく。 Example of Rejection Sampling - yasuhisa's blog Example of importance sampling - yasuhisa's blog Example of Metropolis Hastings Algo…

Gibbs Sampler Algorithmによって多変量正規分布からのサンプル抽出を行なう

もちろんRで。 2変数正規分布 2変数の正規分布からのサンプリングはここのをそのまままねすると B <- 2000 x1 <- rep(NA,B+1) x2 <- rep(NA,B+1) x1[1] <- -2 x2[1] <- 1 for (i in 1:B){ x1[i+1] <- rnorm(1,1+0.7*(x2[i]-2),sqrt(1-0.7^2)) x2[i+1] <- rno…

Metropolis Hastings Algorithmの続き

前回のところではガンベル分布からの乱数が生成できないとして、独立連鎖において正規分布を使って乱数生成をやった。しかし、上の記事の通り、全区間に渡って0以上となるような分布ならとりあえずできるらしいということが分かった。なので今回は独立連鎖の…

Gibbs Sampler and Metropolis Hastings Algorithm

先輩と二人でExample of Metropolis Hastings Algorithm - yasuhisa's blog付近の続き。メトロポリスとかがどういうアルゴリズムかとかは分かった。実装もして欲しい乱数列を取得することもできた。が、なんでこれでうまく行くのかが分からない。よし、調べ…

Example of Metropolis Hastings Algorithm

ということで初めてのMCMCだよ。パターン認識と機械学習 下 - ベイズ理論による統計的予測のP252-253、マルコフ連鎖モンテカルロ法 (統計ライブラリー)のP9-14とここの資料を参考にしながら作りました。What to doガンベル分布(二重指数分布)からのサンプリ…

Example of Rejection Sampling

をやってみるてすと。PRMLの第11章サンプリング法のところのP242とかです。Suppose we can not make the random samples which come from gamma distribution. Let us consider making them from the cauchy distribution. So proposal distribution is the …

Example of importance sampling

マルコフ連鎖モンテカルロ法 (統計ライブラリー)の1.2.3を参考にしつつ。rejection samplingでは、target distributionからの乱数を間接的に取ってきているわけですが、今度はxの関数h(x)の期待値に興味があるとしましょう。で、積分が難しいので、サンプリ…