R: knitrをためす

RStudioからknitrでレポートを自動作成してみた - Take a Risk: 林岳彦の研究メモをよんで、さっそく自分でもためしてみる。

RStudio 0.96.330を使用。こういうRmdファイルをかいてみた。

rjagsを使用する
========================================================
## 準備

パッケージの呼び出し。R2WinBUGSはデータとモデルファイルのために呼び出す。
```{r packages}
library(rjags)
library(R2WinBUGS)
```

乱数のタネを設定する。
```{r seed}
set.seed(119)
```

モデルファイル
```{r modelFile}
model.file <- system.file(package = "R2WinBUGS", "model", "schools.txt")
```

データの設定
```{r data}
data(schools)
J <- nrow(schools) 
y <- schools$estimate 
sigma.y <- schools$sd 
data <- list(J = J, y = y, sigma.y = sigma.y)
```

初期値の設定
```{r inits}
inits <- vector("list", 3)
inits[[1]] <- list(.RNG.name = "base::Mersenne-Twister",
                   .RNG.seed = 123,
                   theta = rnorm(J, 0, 100),
                   mu.theta = rnorm(1, 0, 100), 
                   sigma.theta = runif(1, 0, 100))
inits[[2]] <- list(.RNG.name = "base::Mersenne-Twister",
                   .RNG.seed = 1234,
                   theta = rnorm(J, 0, 100),
                   mu.theta = rnorm(1, 0, 100), 
                   sigma.theta = runif(1, 0, 100))
inits[[3]] <- list(.RNG.name = "base::Mersenne-Twister",
                   .RNG.seed = 12345,
                   theta = rnorm(J, 0, 100),
                   mu.theta = rnorm(1, 0, 100), 
                   sigma.theta = runif(1, 0, 100))
```

保存するパラメータの設定
```{r parameters}
parameters <- c("theta", "mu.theta", "sigma.theta")
```

## 実行

MCMCの計算を実行する。
```{r MCMC}
model <- jags.model(model.file, data, inits,
                    n.chains = 3, n.adapt = 1000)
post <- coda.samples(model, parameters,
                     n.iter = 5000, thin = 5)
```

## 結果の表示
Gelman-Rubinの収束判定
```{r Rhat}
gelman.diag(post)
```

サマリーの表示
```{r summary}
summary(post)
```

軌跡と確率密度を表示する。
```{r plot, fig.width=7, fig.height=6}
plot(post)
```

このRmdファイルの出力はこちら→rjags_knitr.html

タグ：R