前回、グラフを描くのにつかったコードを関数化してみる。例によってエラーチェックなどは最低限しかしていない。 library(ggplot2) qplot.mcmc.list <- function(data, type = "trace", nrow = NULL, ncol = NULL, xlab = "Interactions", ylab = "Value", col = "black") { if (class(data) != "mcmc.list") { stop("class of data is not mcmc.list.") } if (type != "trace" & type != "density") { stop("type is not 'trace' nor 'density'.") } n.var <- length(varnames(data)) x

R2WinBUGSの使い方 このページは、R上からWinBUGSを操作することを可能にするパッケージ、R2WinBUGSの使い方を解説するページです。 モデル式の部分を加筆中 累積訪問者数（from 2008/12/17）： - 今日来てくれた人： - はじめに 計算機の進歩により、Bayes統計（というかMCMC計算）が可能になり、その利用に関する関心が急速に高まっています。しかし、実行方法の難しさが多くのユーザーへの拡大を妨げているといえるでしょう。そこでMCMC（Gibbsサンプリングのみですが.......）計算を行うためのWinBUGSが開発されました。WinBUGSは単体で結果の表示、作図などを行ってくれます。しかし、それら以外の結果が取り出したいケースもありますし、何よりRで使えたら論文などへの二次利用も容易となるでしょう。 そこで開発されたのがR2WinBUGSです。名前は

今さら感はあるが、MCMC (Markov Chain Monte Carlo; マルコフ連鎖モンテカルロ)を使えるようになろうと、まずは簡単な例から試してみた。 手始めに、正規乱数から生成した標本の平均と標準偏差を推定してみる。 やはりRを使用。MCMCpackパッケージを あらかじめインストールしておいて、呼び出す。MCMCpack中のMCMCmetrop1R()関数を利用して、メトロポリス法によるMCMC推定をおこなう。 library(MCMCpack) 乱数系列を初期化。 set.seed(1) 平均10、標準偏差3の乱数を1000個生成して、xに入れる。 m <- 10 s <- 3 x <- rnorm(1000, m, s) MCMC推定に使用する関数を用意する。betaは要素数2のベクトル。beta[1]が平均、beta[2]が標準偏差で、betaを推定する。関数の返り値

秋の某統計高座のテキスト作成のため、もう一度MCMCを最初からやってみる。 例題1: 正規分布することがわかっている ある母集団から、X = (5.89, 5.69, 4.71, 6.73, 4.90, 3.34, 4.60, 4.00) という標本が得られたとき、その母集団の平均と標準偏差をベイズ推定する。 コード ## ## Sample 1 ## library(R2WinBUGS) ## data X <- c(5.89, 5.69, 4.71, 6.73, 4.90, 3.34, 4.60, 4.00) ## model model <- function() { for (i in 1:N) { X[i] ~ dnorm(mu, tau); } ## priors mu ~ dnorm(0.0, 1.0E-6); tau ~ dgamma(1.0E-3, 1.0E-3); si

Hatena ID is an account used for various Hatena services.

今日はMCMC法についての解説です。 メモ程度のものですが、ご参考になれば幸いです。 日本語の良本はこれ。 マルコフ連鎖モンテカルロ法 (統計ライブラリー) 作者: 豊田秀樹出版社/メーカー: 朝倉書店発売日: 2008/05/01メディア: 単行本購入: 11人 クリック: 168回この商品を含むブログ (13件) を見る 有名な解説論文： Sampling-Based Approaches to Calculating Marginal Densities. Gelfand AE and Afrian F. M. Smith. Journal of the American Statistical Association, 85;410:398-409, 1990. 【概念】 Monte Carlo（モンテカルロ法） モンテカルロ：金持ちの町、F1もやってる 興味のある値を「頻度」を使