Bayes

使ってみよう～BUGS(乱数を使ってベイズ推定)～R+openBUGS

BUGSとは

BUGSはBayesian inference Using Gibbs Samplingの略である。

ベイズ推定をするときには次のようなものを使う。

事前分布、データ、事後分布、それらをつなぐ関数など。

つなぐ関数などをごちゃごちゃと書くのは大変なので、モデルとして、事前分布の関数、事後分布の関数、乱雑項の関数などを書いておいて、後は、乱数を発生させ、必要に応じて適切な分布に沿った乱数を使うためにGibbsサンプリングをするという手法である。

openBUGSはそのようなBUGSを支えるオープンソースであり、RはそのopenBUGSを呼び出して実行する。

準備：openBUGSを使えるようにする

openBUGSのダウンロードとインストールをする。ダウンロード先はhttp://www.openbugs.net/w/Downloads

RからopenBUGSを使えるようにする

#install.packages(c("BRugs", "coda","R2WinBUGS","Rcpp"))
library(R2WinBUGS) 

## Loading required package: coda
## Loading required package: boot

library(BRugs) 

## Welcome to BRugs connected to OpenBUGS version 3.2.3

試しに使ってみる

必要なものは、 * Rの実行コマンドと * Rが読み込むBUGS用のモデルファイル

まずは一次線形回帰に使ってみる

非BUGS推定

BUGSを実行する前に、BUGSでない手法での推定をしてみる。

# 適当にデータを作る
# 要素数
N <- 100
# 説明変数xは適当な範囲の一様乱数
x.mean <- 3
x.sd <- 1
x <- rnorm(N,x.mean,x.sd)
# 乱雑項
v <- 1
z <- rnorm(N,0,sqrt(v))
# 従属変数yはxの一次線形の値に乱雑項を加えたもの
beta1 <- 2
beta2 <- 4
y <- beta1 + beta2*x + z
# そのデータを見てみて、lm()で回帰して回帰直線を引きます
# 乱雑項zがxに依存しているので、このような回帰は「正しく」ないけれど、そんなことは知らなかったとして回帰する
lm1 <- lm(y ~ x)
# その結果
summary(lm1)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.49825 -0.68824 -0.04875  0.62549  2.81206 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.6173     0.2993   8.744  6.4e-14 ***
## x             3.7984     0.1009  37.661  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9991 on 98 degrees of freedom
## Multiple R-squared:  0.9354, Adjusted R-squared:  0.9347 
## F-statistic:  1418 on 1 and 98 DF,  p-value: < 2.2e-16

\[ y = \beta_1 \times x + \beta_2 \] のbeta1,beta2 の値がそこそこの精度で推定されていることが、summary(lm1)のCoefficientsのEstimateの出力(Intercept),xとでわかる。

# プロット
plot(x,y,pch=20) 
abline(lm1, col="blue", lwd=2) 

### BUGS推定 さて、これをBUGSでやってみる。

モデルとモデルファイル

BUGSするには、モデルが必要。

モデルはテキストファイルに書いて、Rから読み込ませる。 モデルファイルの意味や、その後に続くRの細かいことは、ひとまず置いて、実行ができるかどうかを確かめるために実施してみよう。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], S)
        mu[i] <- beta1 + beta2 * x[i]
    }
    beta1 ~ dnorm(0, 0.0001)
    beta2 ~ dnorm(0, 0.0001)
    S ~ dgamma(0.0001,0.0001)
}

この内容を“lm1.txt”という名前でRのワーキングディレクトリに保存しておけばよい。

同様のファイルをRの中から作成するには次のようにする。

lm1model <- function(){
 for(i in 1:N) { 
  y[i] ~ dnorm(mu[i], S) 
  mu[i] <- beta1 + beta2 * x[i] 
 } 
 
 beta1 ~ dnorm(0, 0.0001) 
 beta2 ~ dnorm(0, 0.0001) 
 S ~ dgamma(0.0001, 0.0001)
}
file.name <- "lm1.txt"
write.model(lm1model,file.name)

また、MCMCでの実行をするために、推定したいパラメタの初期値を与える必要がある。 beta1,beta2には事前予想値である0を与えればよいだろうし、Sについては、適当に正の値(ここでは1)を与えればよい。

実際にBUGSをRで実行するにはbugs()関数に、しかるべきルールでデータと初期値とモデルファイルのパスとopenBUGSプログラムのパスを与えるとともに、MCMCの実行条件を与える必要がある。

そのコマンドは以下のとおりである。

# データはリストで与える。y~xの回帰のためには説明変数、従属変数、データ数を与える
data1 <- list(y=y, x=x, N=N) 
# 初期値の指定
in1 <- list(beta1=0,beta2=0, S=1) 
inits <- list(in1) 
# 初期値は入れ子になったリストになる
inits

## [[1]]
## [[1]]$beta1
## [1] 0
## 
## [[1]]$beta2
## [1] 0
## 
## [[1]]$S
## [1] 1

# 推定パラメタの名前
param <- c("beta1","beta2","S") 
# テキストファイルでBUGS用のモデルを記載
model1 <- file.name
# 実行
# model1に指定したファイルのパスはもちろんRの実行ディレクトリに合わせる。Rで書きだせば、ワーキングディレクトリに作られているはず
# bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323"という指定はopenBUGSを動かすための指定なので、インストールした先のディレクトリの名前を確認しておく。以下のパスは、2014/03/26にダウンロードしたopenBUGSのデフォルトディレクトリ
# MCMCなので、初期値依存の部分burnin部分(n.burnin)とそれ以降を含めた全繰り返し処理数n.iterとを指定してある。
bugs1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=300, n.burnin=100, n.thin=2, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

実行条件と結果とを表示する。 プロットもしてみる。 興味があるのは各推定パラメタの平均値や中央値、それに区間推定値であるからそれらが示されている。

print(bugs1)

## Inference for Bugs model at "lm1.txt", fit using OpenBUGS,
##  1 chains, each with 300 iterations (first 100 discarded), n.thin = 2
##  n.sims = 100 iterations saved
##           mean  sd  2.5%   25%   50%   75% 97.5%
## beta1      2.6 0.2   2.1   2.4   2.6   2.7   3.0
## beta2      3.8 0.1   3.7   3.8   3.8   3.9   4.0
## S          1.0 0.1   0.8   0.9   1.0   1.1   1.2
## deviance 283.9 1.8 281.8 282.5 283.5 284.9 288.2
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 2.4 and DIC = 286.5
## DIC is an estimate of expected predictive error (lower deviance is better).

plot(bugs1)

bugs()関数の出力の平均・標準偏差・クオンタイル値(0.1,0.5,0.9)を取り出して、桁数を増やして表示してみる。

# 出力bugs1オブジェクトから、結果行列を取り出してpost1オブジェクトとする
post1 <- bugs1$sims.matrix
# post1の各列にパラメタとして採択された値が格納されているのでそれらの記述統計を取る
apply(post1,2,mean) 

##       beta1       beta2           S    deviance 
##   2.5606527   3.8204520   0.9946074 283.8957593

apply(post1,2,sd) 

##      beta1      beta2          S   deviance 
## 0.24469401 0.08232257 0.11738382 1.78444075

apply(post1,2,quantile,c(0,1,0.5,0.9))

##         beta1    beta2         S deviance
## 0%   2.013134 3.661867 0.7069177 281.6340
## 100% 3.035765 4.011907 1.3218921 289.1391
## 50%  2.585789 3.817109 0.9899563 283.4688
## 90%  2.868569 3.928216 1.1400827 286.7318

\[ y = beta1 + beta2 x\] の係数推定がそれなりに良いことがわかる

burnin 後の記録として残っている部分の推移を示す。 値が大きくなったり小さくなったりするので、じゃみじゃみするが、全体の傾向としては、上下に幅がある状態で安定(水平)になっていることがわかる。 水平方向に安定していれば、推定がうまく行っていると考える。 本当はn.iterを増やして、確かに水平方向に安定していることを確かめる方が良い。

par(mfrow=c(3,1)) 
plot(post1[1:bugs1$n.sims,1], xlab="iterations",type="l", col=4) 
plot(post1[1:bugs1$n.sims,2], xlab="iterations",type="l", col=4) 
plot(post1[1:bugs1$n.sims,3], xlab="iterations",type="l", col=4) 

事後分布をヒストグラムで示す。

par(mfrow=c(1,3))
hist(post1[,1],freq=FALSE,xlab="beta[1]",main="Posterior 
distribution of beta1") 
 
hist(post1[,2],freq=FALSE,xlab="beta2",main="Posterior 
distribution of beta2") 

hist(post1[,3],freq=FALSE,xlab="S",main="Posterior 
distribution of S")

par(mfrow=c(1,1))

例で慣れるR+openBUGS

例をなぞりながら、bugs()関数の実効オプション、モデルファイルの意味合い・書き方を確認して行くことにする。

0か1かの確率の推定

今、ある商品Aの後継商品Bを作成したとして、AとBとのどちらが好まれるかを知りたいとする。Bの方がより好まれる確率(pb)は、非常に高くて1かもしれないけれど、逆に全く評価されず、Aの方が好まれる確率が1(Bが好まれる確率が0)かもしれない。 まったく予想がつかないとする。

モニターを募って、何人(n人）かの人にA/Bどちらを好むかを答えてもらった上で、Bをより好む人(nb人）の割合を推定したいとする。

ここでは、問題を単純にするために、必ずどちらか片方を好むものとする。

まず、モニター開始前のpbは0-1区間で一様分布であると想定する。

n人中nb人がBを好むと答える確率は、pbを用いて \[ \frac{n!}{nb!(n-nb)}pb^{nb}(1-pb)^{n-nb} \] と表せるだろう。 二項分布である。

この式のように表されると考えるのは、「各モニターが前後のモニターに影響されていない」などの条件をつけたときに成り立つ式であることなどを考えれば、「モデルとしてこのような式を採用する」と仰々しく述べてもよい。 これがモデルである。

model
{
    nb ~ dbin(p, n)
    p ~ dbeta(1, 1)
}

この内容を“model1.txt”という名前でRのワーキングディレクトリに保存しておけばよい。

同様のファイルをRの中から作成するには次のようにする。

model1 <- function(){
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dbeta(1,1)
}
file.name <- "model1.txt"
write.model(model1,file.name)

data1 <- list(n=10,nb=4) 
in1 <- list(p=0.5) 
inits <- list(in1) 
param <- c("p") 
file.name <- "model1.txt"
model1 <- file.name
n.iter <- 300
bugs1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=n.iter, n.burnin=0, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

print(bugs1)

## Inference for Bugs model at "model1.txt", fit using OpenBUGS,
##  1 chains, each with 300 iterations (first 0 discarded)
##  n.sims = 300 iterations saved
##          mean  sd 2.5% 25% 50% 75% 97.5%
## p         0.4 0.1  0.2 0.3 0.4 0.5   0.7
## deviance  3.8 1.3  2.8 2.9 3.3 4.1   7.4
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = 4.8
## DIC is an estimate of expected predictive error (lower deviance is better).

plot(bugs1)

post1 <- bugs1$sims.matrix
apply(post1,2,mean) 

##         p  deviance 
## 0.4246177 3.7769153

apply(post1,2,sd) 

##         p  deviance 
## 0.1476779 1.2539364

hist(post1[,1],freq=FALSE,xlab="p",main="Posterior 
distribution of p") 

これはベータ分布\[ \beta(nb+1,n-nb+1)\]になるはずである。 ベータ分布乱数と比較してみる。

r.beta <- rbeta(n.iter,data1$nb+1,data1$n-data1$nb+1)
plot(sort(post1[,1]),sort(r.beta))
abline(0,1,col=2)

ではこの単純な例でやっている内容を確認してみる。

モデルファイル

まずはモデルファイルについて確認する。 モデルファイルでは、分布関数を用いて表現する。openBUGSの分布関数は(基本的には)winBUGSのそれと同じである。 winBUGSの分布関数については、http://d.hatena.ne.jp/ryamada22/20140329 を参照。

さて、以下のモデルファイルでは、dbin(),dbeta()という2つの分布関数が用いられている。二項分布関数とベータ関数である。

model
{
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dbeta(1,1)
}

1行目は、Bを好む確率がpのときに、n人に訊いたら、nb人が「Bを好む」ことが二項分布に従って定まることを示している。

今、BUGSを回すときにはnとnbとを与えて、適切なpの分布をしミューレーションで示そうとしていることを示している。

他方、2行目はpとしてどんな分布を事前分布を想定しているかを表している。 2つのパラメタが1,1であるようなベータ分布であることを示している。

実行に際しては、このベータ分布からpをランダムに発生し、それを二項分布に与え、nのもとでnbになるかどうかでそのpの乱数を採択するかどうかを決める。

dbeta(1,1)とは

p <- seq(from=0,to=1,length=1000)
db <- dbeta(p,1,1)
plot(p,db)

でわかるように、0から1までの一様分布であるから、dunif(0,1)という分布を指定しても同じである。

実際にやってみる。

model
{
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dunif(0,1)
}

model1.1 <- function(){
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dunif(0,1)
}
file.name <- "model1.1.txt"
write.model(model1.1,file.name)

data1 <- list(n=10,nb=4) 
in1 <- list(p=0.5) 
inits <- list(in1) 
param <- c("p") 
file.name <- "model1.1.txt"
model1 <- file.name
# 実行
n.iter <- 300
bugs1.1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=n.iter, n.burnin=0, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

一様分布を事前分布に指定しても得られる値の出現順は異なるが、分布としては同じ結果が得られている。

post1.1 <- bugs1.1$sims.matrix
(head(post1[,1]))

## [1] 0.4422682 0.5461071 0.2374524 0.3250296 0.2934158 0.5204837

(head(post1.1[,1]))

## [1] 0.5250736 0.7339078 0.5559581 0.4560897 0.4569803 0.4559864

plot(sort(post1[,1]),sort(post1.1[,1]))
abline(0,1,col=2)

このように同じ事前分布なのに、あえてベータ分布を用いたのは、dbeta(1,1)という一様分布でないような分布を使う場合を考えると、二項分布の共役事前分布として、異なる形の事前分布に変えやすいからである。

実行条件

bugs()関数にはデータ、初期値、パラメタ名、モデルファイル、debug、program、bugs.directoryのほかにいくつかの情報を与える。 それぞれを説明する

• n.chains
• MCMCでは初期値と乱数列を使って処理をするが、初期値を変えて何度か実行する方がよいこともあり、その回数を指定する
• n.iter
• 推定対象は多数の値が作る分布として判断するわけだが、その作るべき値の数を示す。ただし、そのうち、初めの方のn.burnin個は捨てるので、実際に推定される分布として使用可能なのはn.iter-n.burnin個になる
• n.burnin
• 上述したように、初めのn.burnin個は捨てる
• 今回のように一番初め値から、推定に用いることが適切な場合にはn.burnin=0でも構わない
• n.thin
• 非常に多くの値を発生させるのは、その方が、初期値の影響が混入することを避けることができるからであるが、あまりにたくさんだと大変なので、n.iterが大きいときには、何個かおきに記録してそれ以外を捨てることも有効である。n.iter=1ならば全部を記録し、n.iter=2ならば一つおきに捨てる、とそういう具合である。

n.chains,n.iter,n.burnin,n.thinの値を変えながら実行して、出力がどのように変わるかを確認してみる。

n.chains=1でn.iter=100,n.burnin=10,n.thin=3としてみる。

bugs1.2 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=100, n.burnin=10, n.thin=3, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs1.2)

## Inference for Bugs model at "model1.1.txt", fit using OpenBUGS,
##  1 chains, each with 100 iterations (first 10 discarded), n.thin = 3
##  n.sims = 30 iterations saved
##          mean  sd 2.5% 25% 50% 75% 97.5%
## p         0.4 0.1  0.2 0.4 0.4 0.5   0.7
## deviance  3.7 1.5  2.8 2.8 3.2 3.8   7.6
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = 4.8
## DIC is an estimate of expected predictive error (lower deviance is better).

n.iter-n.burnin = n.thin*bugs1$n.simsとなっていることがわかる。

次にn.chains=3、n.thin=2でやってみる。初期値をn.chains通り、指定する必要がある。

in1 <- list(p=0.5)
in2 <- list(p=0.3)
in3 <- list(p=0.7)
inits <- list(in1,in2,in3)
inits

## [[1]]
## [[1]]$p
## [1] 0.5
## 
## 
## [[2]]
## [[2]]$p
## [1] 0.3
## 
## 
## [[3]]
## [[3]]$p
## [1] 0.7

このように初期値を3通りリストとして準備している。

さて、実行してみる。

bugs1.3 <- bugs(data1, inits, param, model.file=model1, 
n.chains=3, n.iter=100, n.burnin=10, n.thin=2, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs1.3)

## Inference for Bugs model at "model1.1.txt", fit using OpenBUGS,
##  3 chains, each with 100 iterations (first 10 discarded), n.thin = 2
##  n.sims = 135 iterations saved
##          mean  sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## p         0.4 0.1  0.2 0.3 0.4 0.5   0.7    1    59
## deviance  3.7 1.2  2.8 2.9 3.3 4.1   7.0    1   140
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = 4.8
## DIC is an estimate of expected predictive error (lower deviance is better).

各chainについてn.iter回をn.burnin込みで実施していることが記載されている。

正規分布の平均と分散とを推定する

比較的ふんだんに発現していると思われる遺伝子mRNAをある検出手法で測定したときの値の分布について推定することを考える。 検出手法についてもあまりよくわかっていないので、平均がいくつで分散がいくつになるか、さっぱり予想がつかないが、正規分布に従うだろうことは予想してもよいとする。

このような場合のモデルは以下のようになる。 “model2.txt”というファイルで保存しておくこととする。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu, S)
    }
    mu ~ dnorm(0, 0.0001)
    S ~ dgamma(0.0001,0.0001)
}

model2 <- function(){
   # Model
   for (i in 1:N) {
        y[i] ~ dnorm(mu, tau)
    }
   
   # Prior
    mu ~ dnorm(0, 0.0001)
    tau ~ dgamma(0.0001,0.0001)
   
}
file.name <- "model2.txt"
write.model(model2,file.name)

このモデル設定では、N個の観察値を一つのmu,tauのセットに基づいて正規分布に従うことを示している。

muには、第1引数(v1 平均)0、第2引数(v2 tau) 0.0001の正規分布を事前分布として想定している。 このv2が小さいということは、分散が大きいということであるので、この事前分布は、0を中心として、非常に裾が長い分布である。 要するに、いろいろな値をほぼ同確率で想定していることになる。

v2の値を変化させて、どのような事前分布になるのかを表示してみる。 v2の値が小さくなるほど平坦な分布となることを、実数と対数とで示す。

横軸の中央付近にピークがあることに着目するのではなく、そこにピークはあるものの、それはごく限られた幅に納まり、それ以外の広い範囲で平坦になっていることが、この分布を採用した理由である。

x <- seq(from=-100,to=100,length=1000)
v2s <- 10^(0:(-5))
d <- matrix(0,length(x),length(v2s))
for(i in 1:length(v2s)){
  v2 <- v2s[i]
  d[,i] <- dnorm(x,0,1/sqrt(v2))
}
par(mfcol=c(1,2))
matplot(x,d,type="l")
matplot(x,log(d),type="l")

par(mfcol=c(1,1))

また、標準偏差の逆数としのtauの事前分布はガンマ分布を用いている。 ガンマ分布が選ばれているのは、tauには正の値をとりたいからである。 そのような理由でガンマ分布を採用したうえで、次にその形を決めたい。

dgamma(a,b)で指定されるガンマ分布は、平均が\[\frac{a}{b}\]分散が\[\frac{a}{b^2}\]であるが、事前分布としてはさまざまな値をとりたい、言い換えると分散を大きくしたい。 ガンマ分布の分散は\[\frac{a}{b^2}\]であるので、a=bとして小さい値を与えることで、平均は1で分散の大きい分布が得られる。 これが、tauの事前分布の選び方である。

事前分布としての指定に際してはa,bに等しく、小さい値を与えているが、これは、平均\(\frac{a}{b}=\frac{a}{a}=1\)が1で分散が大きい分布である。

a,b=aの値を色々に変えてガンマ分布の様子をプロットしてみる。

x <- seq(from=0,to=20,length=1000)
as <- 10^(0:(-7))
d <- matrix(0,length(x),length(as))
for(i in 1:length(as)){
  d[,i] <- dgamma(x,as[i],as[i])
}
par(mfcol=c(1,2))
matplot(x,d,type="l")
matplot(x,log(d),type="l")

par(mfcol=c(1,1))

0付近にピークが残るも、aが小さくなるにつれ、平坦な分布となることがわかる。 ではデータを作成して推定してみる。

N <- 100
mu <- 5
sd <- 2
y <- rnorm(N,mu,sd)
data1 <- list(y=y,N=N)
in1 <- list(mu=0,tau=1) 
inits <- list(in1) 
param <- c("mu","tau") 
file.name <- "model2.txt"
model2 <- file.name
n.iter <- 1000
bugs2 <- bugs(data1, inits, param, model.file=model2, 
n.chains=1, n.iter=n.iter, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

print(bugs2)

## Inference for Bugs model at "model2.txt", fit using OpenBUGS,
##  1 chains, each with 1000 iterations (first 100 discarded)
##  n.sims = 900 iterations saved
##           mean  sd  2.5%   25%   50%   75% 97.5%
## mu         5.2 0.2   4.7   5.0   5.1   5.3   5.6
## tau        0.2 0.0   0.2   0.2   0.2   0.2   0.3
## deviance 439.3 1.9 437.4 437.9 438.8 440.1 444.4
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 2.0 and DIC = 441.3
## DIC is an estimate of expected predictive error (lower deviance is better).

ここで推定されたのは\[\tau=\frac{1}{sd^2}\]であるので、データ生成時に与えたsd=2に近い値が推定されたかを確認するために計算しておく。 また、サンプル平均、サンプルに基づくsdも計算しておく。

post2 <- bugs2$sims.matrix
mean(post2[,1])

## [1] 5.150721

mean(1/sqrt(post2[,2]))

## [1] 2.173052

# サンプル平均とサンプルからの標準偏差
mean(y)

## [1] 5.150687

sd(y)

## [1] 2.165273

その上で分布を示す。

赤は真値、緑はサンプルから算出した平均とsdである。

par(mfcol=c(1,2))
hist(post2[,1],freq=FALSE,xlab="mu",main="Posterior 
distribution of mu")
abline(v=mu,col=2)
abline(v=mean(y),col=3)
hist(1/sqrt(post2[,2]),freq=FALSE,xlab="sd",main="Posterior 
distribution of sd")
abline(v=sd,col=2)
abline(v=sd(y),col=3)

par(mfcol=c(1,1))

線形回帰をやってみる

\[ \begin{equation} y = beta0 + beta1 \times x1 + beta2 \times x2 + \epsilon\times z\\ z ~ N(0,1) \end{equation} \] という線形回帰をしてみよう。冒頭の「ひとまず回す」の例では、説明変数が1つの線形回帰だったが、今回は2つにしてみる。

具体性を持たせるためにx1を年齢、x2を体重、yを最高血圧として考えよう。 値とその分布の妥当性はここでは問わ無いことにする。

まずはデータを作る。

N <- 500
# 年齢
x1 <- sample(20:100,N,replace=TRUE)
# 体重
x2 <- rnorm(N,65,10)

beta0 <- 30
beta1 <- 0.5
beta2 <- 1
epsilon <- 10

y <- beta0 + beta1*x1+beta2*x2+epsilon * rnorm(N)
hist(y)

# yの値で色を付けて、散布図にする
st.y <- (y-min(y))/(max(y)-min(y))
col <- rgb(st.y,1-st.y,1)
plot(x1,x2,pch=20,col=col,xlab="age",ylab="weight",main="systolic pressure")

モデルファイルを作り、“model3.txt”として保存する。

項の数が増えたが、N人それぞれについて、個人の年齢x1[i]と体重x2[i]に応じて、最高血圧の予想値mu[i]を求めその値を平均とした正規分布に従うことを示してる。 パラメタbeta0,beta1,beta2,epsilonに事前予想が立たないと仮定すれば、 beta0,beta1,beta2は0を中心とした分散の大きな正規分布を事前分布とし、 乱雑項も分散が不明な正規分布とするのがよいので、そのような分布が指定してある。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], tau)
        mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
    }
    beta0 ~ dnorm(0,0.0001)
    beta1 ~ dnorm(0,0.0001)
    beta2 ~ dnorm(0,0.0001)
    tau ~ dgamma(0.0001,0.0001)
}

model3 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], tau)
        mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
    }
    beta0 ~ dnorm(0,0.0001)
    beta1 ~ dnorm(0,0.0001)
    beta2 ~ dnorm(0,0.0001)
    tau ~ dgamma(0.0001,0.0001)
}
file.name <- "model3.txt"
write.model(model3,file.name)

さて、実行する。

data1 <- list(y=y,x1=x1,x2=x2,N=N) 
in1 <- list(beta0=0,beta1=0,beta2=0,tau=1) 
inits <- list(in1) 
param <- c("beta0","beta1","beta2","tau") 
file.name <- "model3.txt"
model1 <- file.name
# 実行
bugs3 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs3)

## Inference for Bugs model at "model3.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##            mean  sd   2.5%    25%    50%    75%  97.5%
## beta0      30.7 3.4   23.6   28.2   31.2   33.2   36.2
## beta1       0.5 0.0    0.5    0.5    0.5    0.5    0.5
## beta2       1.0 0.0    0.9    1.0    1.0    1.0    1.1
## tau         0.0 0.0    0.0    0.0    0.0    0.0    0.0
## deviance 3717.4 2.8 3713.7 3715.2 3716.6 3718.8 3724.8
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 4.2 and DIC = 3722.0
## DIC is an estimate of expected predictive error (lower deviance is better).

post3 <- bugs3$sims.matrix
apply(post3,2,quantile,c(0.1,0.5,0.9))

##        beta0     beta1     beta2        tau deviance
## 10% 25.86782 0.4815349 0.9378688 0.00930471 3714.432
## 50% 31.21748 0.5048564 0.9866686 0.01004688 3716.610
## 90% 34.67386 0.5293288 1.0588734 0.01092468 3721.427

beta0,beta1,beta2が良い推定値になっていることがわかる。

招待選手と一般参加選手の2群～混合分布～

今、全部でN人がエントリーしている100m走大会があるとする。 このうち、pi[1]の割合で招待選手(群 g=1)、pi[2]の割合で一般参加選手(群 g=2)であるとする。

招待選手のタイムは平均mm[1]標準偏差ss[1]、一般参加選手はmm[2],ss[2]の正規分布に従うとする。 ただし、この招待選手というのは、100m走の速さを基準に招待されたわけではなく、マラソンが得意であることで招待されたという。 招待選手はアスリートなので、一般参加選手よりも100m走のパフォーマンスがよいようにも思えるし、大差ないかもしれないし、かえって遅いかもしれない。 pi[1],pi[2]の事前確率は0-1の範囲で全く不明だという。

データを作成してみる。

N <- 1000
pi <- c(0.2,0.8)
n1 <- N*pi[1]
n2 <- N-n1

mm <- c(12,16)
ss <- c(1,2)
# 100m走のタイム
y1 <- rnorm(n1,mm[1],ss[1])
y2 <- rnorm(n2,mm[2],ss[2])
y12 <- c(y1,y2)
# グループ
g12 <- c(rep(1,n1),rep(2,n2))
# 選手番号をランダムにつける
ord <- sample(1:N)
# 選手番号順に並べ替える。これが観察データ
y <- y12[ord]
g <- g12[ord]

h <- hist(y,plot=FALSE)
hist(y1,breaks=h$breaks,density=20,col=2,ylim=c(0,max(h$counts)))
hist(y2,breaks=h$breaks,density=17,col=3,add=TRUE)

ではモデルを立ててみる。

各選手のグループは確率pにより1か2に決まる。 そのpは一様分布を事前確率とする。 各選手のタイムはグループ1か2かによって平均、標準偏差を変える必要があるが、それらに基づく正規分布に従う。 2つのグループのタイムを決める平均と標準偏差はm1,m2,s1,s2とするが、m1,m2は9秒から30秒の範囲で一様分布を想定しよう。s1,s2は全く不明として、正の値を分散の大きいガンマ分布で与えることとする。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau[i])
        m[i] <- mm[g[i]]
        tau[i] <- tt[g[i]]
        g[i] ~ dcat(pi[])
    }
    pi[1:2] ~ ddirch(alpha[])
    mm[1] ~ dunif(9, 30)
    mm[2] ~ dunif(9, 30)
    tt[1] ~ dgamma(1.00000E-04, 1.00000E-04)
    tt[2] ~ dgamma(1.00000E-04, 1.00000E-04)
}

model4 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau[i])
        m[i] <- mm[g[i]]
        tau[i] <- tt[g[i]]
        g[i] ~ dcat(pi[])
    }
    pi[1:2] ~ ddirch(alpha[])
    mm[1] ~ dunif(9,30)
    mm[2] ~ dunif(9,30)
    tt[1] ~ dgamma(0.0001,0.0001)
    tt[2] ~ dgamma(0.0001,0.0001)
}
file.name <- "model4.txt"
write.model(model4,file.name)

さて実行しよう。

alpha <- c(1,1)
data1 <- list(y=sort(y),g=c(1,rep(NA,N-2),2),alpha=alpha,N=N) 
in1 <- list(pi=c(0.5,0.5),mm=c(15,15),tt=c(1,1)) 
inits <- list(in1) 
param <- c("pi","mm","tt") 
file.name <- "model4.txt"
model1 <- file.name
# 実行
bugs4 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs4)

## Inference for Bugs model at "model4.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##            mean   sd   2.5%    25%    50%    75%  97.5%
## pi[1]       0.2  0.0    0.1    0.2    0.2    0.2    0.3
## pi[2]       0.8  0.0    0.7    0.8    0.8    0.8    0.9
## mm[1]      12.0  0.2   11.7   11.9   12.0   12.1   12.5
## mm[2]      15.8  0.2   15.5   15.7   15.8   15.9   16.2
## tt[1]       1.1  0.3    0.6    0.9    1.1    1.3    1.8
## tt[2]       0.2  0.0    0.2    0.2    0.2    0.2    0.3
## deviance 4032.0 93.4 3855.6 3970.2 4033.5 4093.2 4227.0
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 4364.2 and DIC = 8396.2
## DIC is an estimate of expected predictive error (lower deviance is better).

2群の割合と2群の平均がよい感じに推定されていることがわかる。 推定されたtauから標準偏差を計算し、データ生成に用いたそれがよく推定されていることも確認できる。

post4 <- bugs4$sims.matrix
apply(post4,2,quantile,c(0.1,0.5,0.9))

##         pi[1]     pi[2]    mm[1]    mm[2]     tt[1]     tt[2] deviance
## 10% 0.1287006 0.7586582 11.78398 15.60043 0.7605146 0.2069366 3906.375
## 50% 0.1729164 0.8270836 11.97887 15.78823 1.1006213 0.2319064 4033.527
## 90% 0.2413418 0.8712994 12.25860 16.05125 1.5520641 0.2665965 4150.270

mean(1/sqrt(post4[,5]))

## [1] 0.9675272

mean(1/sqrt(post4[,6]))

## [1] 2.07019

n次線形回帰をしてみる

\[ y = \sum_{i=0}^n b_i x^i + \epsilon z \] (ただしzは標準正規分布に従う乱雑項) のようなデータから、n+1個のbの値を推定してみる。

まずはデータ作成。ただし、三角関数を使って作成する。

N <- 100
x <- sort(rnorm(N))
y <- 3*sin(x)
epsilon <- 0.5
z <- rnorm(N)
y <- y + epsilon * z
plot(x,y,pch=20)

モデルファイルを作る。

model5 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau)
        m[i] <- b[1]+b[2]*x[i] + b[3]*x[i]*x[i] + b[4] * pow(x[i],3)
    }
    for(i in 1:(n+1)){
      b[i] ~ dnorm(0,0.0001)
    }
    tau ~ dgamma(0.0001,0.0001)
}
file.name <- "model5.txt"
write.model(model5,file.name)

実行する。

n <- 3
data1 <- list(y=y,x=x,N=N,n=n) 
in1 <- list(b=rep(0,n+1),tau=1) 
inits <- list(in1) 
param <- c("b","tau") 
file.name <- "model5.txt"
model1 <- file.name
# 実行
bugs5 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Warning in readChar(modelFile, 10^3): can only read in bytes in a non-UTF-8
## MBCS locale

## model is syntactically correct
## data loaded
## model compiled
## Initializing chain 1: 
## model is initialized
## model is already initialized

## Sampling has been started ...

## 100 updates took 0 s
## deviance set
## monitor set for variable 'b'
## monitor set for variable 'tau'
## monitor set for variable 'deviance'
## 1000 updates took 0 s

print(bugs5)

## Inference for Bugs model at "model5.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##           mean  sd  2.5%   25%   50%   75% 97.5%
## b[1]       0.1 0.1   0.0   0.1   0.1   0.2   0.3
## b[2]       2.8 0.1   2.6   2.8   2.8   2.9   3.0
## b[3]      -0.1 0.0  -0.2  -0.2  -0.1  -0.1   0.0
## b[4]      -0.4 0.0  -0.4  -0.4  -0.4  -0.3  -0.3
## tau        4.0 0.6   2.9   3.6   4.0   4.4   5.2
## deviance 145.5 3.1 141.3 143.2 144.9 147.1 153.1
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 4.9 and DIC = 150.4
## DIC is an estimate of expected predictive error (lower deviance is better).

推定された係数を用いて、推定線を描いてみる。

post5 <- bugs5$sims.matrix
out5.coef <- apply(post5,2,mean)
plot(x,y,pch=20)
x. <- seq(from=min(x)-0.1*(max(x)-min(x)),to=max(x)+0.1*(max(x)-min(x)),length=100)
pred.y <- rep(out5.coef[1],length(x.))
for(i in 2:(n+1)){
  pred.y <- pred.y + out5.coef[i] * x.^(i-1)
}
points(x.,pred.y,col=2,type="l")

カーネル分布推定

N個の値を観測してそこに密度分布を推定するとする。 正規分布をカーネル関数としてカーネル法でやるなら、個々の観測に、ある分散の正規分布を与え、その混合分布として求めることとなる。

混合分布なら、すでに、「100m走大会～招待選手と一般参加選手」の例でやったので、それを少し改変すればよい。

100m走大会では2群だったが、今度は、データ数Nのすべてを個別の群として扱うので群数はNである。

それぞれの群の寄与度は1/Nで固定である。

推定したいのは、個々の値を中心とした正規分布のその分散をいくつにするかである。

すべての値について分散の値を共通にするのであれば、推定パラメタは一つだけである。

それを書けば以下の通り。

yは観測値のベクトル。 mmもyと同じベクトルであるが、これはBUGS中に採取されるのではなくて、決まった値。 tauは推定パラメタである。 gはN個の群のラベルである。 モデルファイルの中では、個々の値が、1,…,Nのいずれかの群に属するかを等確率で採取し、取られた群g[i]に応じて、正規分布の平均値m[i]が決まり、tauに応じて値が採取している。

piは1,…,Nを等確率で選ぶために与えるベクトルである。

model6 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau)
        m[i] <- mm[g[i]]
        g[i] ~ dcat(pi[])
    }

    tau ~ dgamma(0.0001,0.0001)
}
file.name <- "model6.txt"
write.model(model6,file.name)

データを作る。

y <- c(rnorm(100),runif(100,3,8),rexp(200,4)+10)
y <- sort(y)
plot(density(y))

N <- length(y)
alpha <- rep(100,N)
data1 <- list(y=y,mm=y,N=N,g=1:N,pi=rep(1/N,N)) 
in1 <- list(tau=1) 
inits <- list(in1) 
param <- c("tau") 
file.name <- "model6.txt"
model1 <- file.name
# 実行
bugs6 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Warning in readChar(modelFile, 10^3): can only read in bytes in a non-UTF-8
## MBCS locale

## model is syntactically correct
## data loaded
## model compiled
## Initializing chain 1: 
## model is initialized
## model is already initialized

## Sampling has been started ...

## 100 updates took 0 s
## deviance set
## monitor set for variable 'tau'
## monitor set for variable 'deviance'
## 1000 updates took 0 s

print(bugs6)

## Inference for Bugs model at "model6.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##               mean       sd      2.5%       25%       50%       75%
## tau      2000861.2 139442.3 1737639.4 1902228.9 2000254.2 2094249.0
## deviance    -274.3     28.0    -325.8    -293.6    -275.2    -255.1
##              97.5%
## tau      2270276.1
## deviance    -218.9
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = -273.4
## DIC is an estimate of expected predictive error (lower deviance is better).

使ってみよう～BUGS(乱数を使ってベイズ推定)～R+openBUGS

BUGSとは

BUGSはBayesian inference Using Gibbs Samplingの略である。

ベイズ推定をするときには次のようなものを使う。

事前分布、データ、事後分布、それらをつなぐ関数など。

つなぐ関数などをごちゃごちゃと書くのは大変なので、モデルとして、事前分布の関数、事後分布の関数、乱雑項の関数などを書いておいて、後は、乱数を発生させ、必要に応じて適切な分布に沿った乱数を使うためにGibbsサンプリングをするという手法である。

openBUGSはそのようなBUGSを支えるオープンソースであり、RはそのopenBUGSを呼び出して実行する。

準備：openBUGSを使えるようにする

openBUGSのダウンロードとインストールをする。ダウンロード先はhttp://www.openbugs.net/w/Downloads ## RからopenBUGSを使えるようにする

#install.packages(c("BRugs", "coda", "R2WinBUGS","Rcpp"))
library(R2WinBUGS) 
library(BRugs) 

試しに使ってみる

必要なものは、 * Rの実行コマンドと * Rが読み込むBUGS用のモデルファイル

まずは一次線形回帰に使ってみる

非BUGS推定

BUGSを実行する前に、BUGSでない手法での推定をしてみる。

# 適当にデータを作る
# 要素数
N <- 100
# 説明変数xは適当な範囲の一様乱数
x.mean <- 3
x.sd <- 1
x <- rnorm(N,x.mean,x.sd)
# 乱雑項
v <- 1
z <- rnorm(N,0,sqrt(v))
# 従属変数yはxの一次線形の値に乱雑項を加えたもの
beta1 <- 2
beta2 <- 4
y <- beta1 + beta2*x + z
# そのデータを見てみて、lm()で回帰して回帰直線を引きます
# 乱雑項zがxに依存しているので、このような回帰は「正しく」ないけれど、そんなことは知らなかったとして回帰する
lm1 <- lm(y ~ x)
# その結果
summary(lm1)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92730 -0.63625 -0.04369  0.67342  2.33084 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.77788    0.28092   6.329 7.44e-09 ***
## x            4.06809    0.08977  45.315  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9591 on 98 degrees of freedom
## Multiple R-squared:  0.9545, Adjusted R-squared:  0.954 
## F-statistic:  2053 on 1 and 98 DF,  p-value: < 2.2e-16

\[ y = \beta_1 \times x + \beta_2 \] のbeta1,beta2 の値がそこそこの精度で推定されていることが、summary(lm1)のCoefficientsのEstimateの出力(Intercept),xとでわかる。

# プロット
plot(x,y,pch=20) 
abline(lm1, col="blue", lwd=2) 

### BUGS推定 さて、これをBUGSでやってみる。

モデルとモデルファイル

BUGSするには、モデルが必要。

モデルはテキストファイルに書いて、Rから読み込ませる。 モデルファイルの意味や、その後に続くRの細かいことは、ひとまず置いて、実行ができるかどうかを確かめるために実施してみよう。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], S)
        mu[i] <- beta1 + beta2 * x[i]
    }
    beta1 ~ dnorm(0, 0.0001)
    beta2 ~ dnorm(0, 0.0001)
    S ~ dgamma(0.0001,0.0001)
}

この内容を“lm1.txt”という名前でRのワーキングディレクトリに保存しておけばよい。

同様のファイルをRの中から作成するには次のようにする。

lm1model <- function(){
 for(i in 1:N) { 
  y[i] ~ dnorm(mu[i], S) 
  mu[i] <- beta1 + beta2 * x[i] 
 } 
 
 beta1 ~ dnorm(0, 0.0001) 
 beta2 ~ dnorm(0, 0.0001) 
 S ~ dgamma(0.0001, 0.0001)
}
file.name <- "lm1.txt"
write.model(lm1model,file.name)

また、MCMCでの実行をするために、推定したいパラメタの初期値を与える必要がある。 beta1,beta2には事前予想値である0を与えればよいだろうし、Sについては、適当に正の値(ここでは1)を与えればよい。

実際にBUGSをRで実行するにはbugs()関数に、しかるべきルールでデータと初期値とモデルファイルのパスとopenBUGSプログラムのパスを与えるとともに、MCMCの実行条件を与える必要がある。

そのコマンドは以下のとおりである。

# データはリストで与える。y~xの回帰のためには説明変数、従属変数、データ数を与える
data1 <- list(y=y, x=x, N=N) 
# 初期値の指定
in1 <- list(beta1=0,beta2=0, S=1) 
inits <- list(in1) 
# 初期値は入れ子になったリストになる
inits

## [[1]]
## [[1]]$beta1
## [1] 0
## 
## [[1]]$beta2
## [1] 0
## 
## [[1]]$S
## [1] 1

# 推定パラメタの名前
param <- c("beta1","beta2","S") 
# テキストファイルでBUGS用のモデルを記載
model1 <- file.name
# 実行
# model1に指定したファイルのパスはもちろんRの実行ディレクトリに合わせる。Rで書きだせば、ワーキングディレクトリに作られているはず
# bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323"という指定はopenBUGSを動かすための指定なので、インストールした先のディレクトリの名前を確認しておく。以下のパスは、2014/03/26にダウンロードしたopenBUGSのデフォルトディレクトリ
# MCMCなので、初期値依存の部分burnin部分(n.burnin)とそれ以降を含めた全繰り返し処理数n.iterとを指定してある。
bugs1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=300, n.burnin=100, n.thin=2, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

実行条件と結果とを表示する。 プロットもしてみる。 興味があるのは各推定パラメタの平均値や中央値、それに区間推定値であるからそれらが示されている。

print(bugs1)

## Inference for Bugs model at "lm1.txt", fit using OpenBUGS,
##  1 chains, each with 300 iterations (first 100 discarded), n.thin = 2
##  n.sims = 100 iterations saved
##           mean  sd  2.5%   25%   50%   75% 97.5%
## beta1      1.7 0.2   1.3   1.5   1.8   1.9   2.1
## beta2      4.1 0.1   4.0   4.0   4.1   4.1   4.2
## S          1.1 0.1   0.8   1.0   1.1   1.2   1.3
## deviance 275.7 1.8 273.6 274.3 275.3 276.7 279.9
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 2.4 and DIC = 278.4
## DIC is an estimate of expected predictive error (lower deviance is better).

plot(bugs1)

bugs()関数の出力の平均・標準偏差・クオンタイル値(0.1,0.5,0.9)を取り出して、桁数を増やして表示してみる。

# 出力bugs1オブジェクトから、結果行列を取り出してpost1オブジェクトとする
post1 <- bugs1$sims.matrix
# post1の各列にパラメタとして採択された値が格納されているのでそれらの記述統計を取る
apply(post1,2,mean) 

##      beta1      beta2          S   deviance 
##   1.726603   4.087248   1.079478 275.707267

apply(post1,2,sd) 

##     beta1     beta2         S  deviance 
## 0.2300033 0.0734137 0.1273984 1.7827548

apply(post1,2,quantile,c(0,1,0.5,0.9))

##         beta1    beta2         S deviance
## 0%   1.216159 3.943534 0.7673317 273.4430
## 100% 2.181787 4.259980 1.4338640 280.9098
## 50%  1.758131 4.083649 1.0741993 275.2811
## 90%  2.024725 4.184156 1.2374160 278.5672

\[ y = beta1 + beta2 x\] の係数推定がそれなりに良いことがわかる

burnin 後の記録として残っている部分の推移を示す。 値が大きくなったり小さくなったりするので、じゃみじゃみするが、全体の傾向としては、上下に幅がある状態で安定(水平)になっていることがわかる。 水平方向に安定していれば、推定がうまく行っていると考える。 本当はn.iterを増やして、確かに水平方向に安定していることを確かめる方が良い。

par(mfrow=c(3,1)) 
plot(post1[1:bugs1$n.sims,1], xlab="iterations",type="l", col=4) 
plot(post1[1:bugs1$n.sims,2], xlab="iterations",type="l", col=4) 
plot(post1[1:bugs1$n.sims,3], xlab="iterations",type="l", col=4) 

事後分布をヒストグラムで示す。

par(mfrow=c(1,3))
hist(post1[,1],freq=FALSE,xlab="beta[1]",main="Posterior 
distribution of beta1") 
 
hist(post1[,2],freq=FALSE,xlab="beta2",main="Posterior 
distribution of beta2") 

hist(post1[,3],freq=FALSE,xlab="S",main="Posterior 
distribution of S")

par(mfrow=c(1,1))

例で慣れるR+openBUGS

例をなぞりながら、bugs()関数の実効オプション、モデルファイルの意味合い・書き方を確認して行くことにする。

0か1かの確率の推定

今、ある商品Aの後継商品Bを作成したとして、AとBとのどちらが好まれるかを知りたいとする。Bの方がより好まれる確率(pb)は、非常に高くて1かもしれないけれど、逆に全く評価されず、Aの方が好まれる確率が1(Bが好まれる確率が0)かもしれない。 まったく予想がつかないとする。

モニターを募って、何人(n人）かの人にA/Bどちらを好むかを答えてもらった上で、Bをより好む人(nb人）の割合を推定したいとする。

ここでは、問題を単純にするために、必ずどちらか片方を好むものとする。

まず、モニター開始前のpbは0-1区間で一様分布であると想定する。

n人中nb人がBを好むと答える確率は、pbを用いて \[ \frac{n!}{nb!(n-nb)}pb^{nb}(1-pb)^{n-nb} \] と表せるだろう。 二項分布である。

この式のように表されると考えるのは、「各モニターが前後のモニターに影響されていない」などの条件をつけたときに成り立つ式であることなどを考えれば、「モデルとしてこのような式を採用する」と仰々しく述べてもよい。 これがモデルである。

model
{
    nb ~ dbin(p, n)
    p ~ dbeta(1, 1)
}

この内容を“model1.txt”という名前でRのワーキングディレクトリに保存しておけばよい。

同様のファイルをRの中から作成するには次のようにする。

model1 <- function(){
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dbeta(1,1)
}
file.name <- "model1.txt"
write.model(model1,file.name)

data1 <- list(n=10,nb=4) 
in1 <- list(p=0.5) 
inits <- list(in1) 
param <- c("p") 
file.name <- "model1.txt"
model1 <- file.name
n.iter <- 300
bugs1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=n.iter, n.burnin=0, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

print(bugs1)

## Inference for Bugs model at "model1.txt", fit using OpenBUGS,
##  1 chains, each with 300 iterations (first 0 discarded)
##  n.sims = 300 iterations saved
##          mean  sd 2.5% 25% 50% 75% 97.5%
## p         0.4 0.1  0.2 0.3 0.4 0.5   0.7
## deviance  3.8 1.3  2.8 2.9 3.3 4.1   7.4
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = 4.8
## DIC is an estimate of expected predictive error (lower deviance is better).

plot(bugs1)

post1 <- bugs1$sims.matrix
apply(post1,2,mean) 

##         p  deviance 
## 0.4246177 3.7769153

apply(post1,2,sd) 

##         p  deviance 
## 0.1476779 1.2539364

hist(post1[,1],freq=FALSE,xlab="p",main="Posterior 
distribution of p") 

これはベータ分布\[ \beta(nb+1,n-nb+1)\]になるはずである。 ベータ分布乱数と比較してみる。

r.beta <- rbeta(n.iter,data1$nb+1,data1$n-data1$nb+1)
plot(sort(post1[,1]),sort(r.beta))
abline(0,1,col=2)

ではこの単純な例でやっている内容を確認してみる。

モデルファイル

まずはモデルファイルについて確認する。 モデルファイルでは、分布関数を用いて表現する。openBUGSの分布関数は(基本的には)winBUGSのそれと同じである。 winBUGSの分布関数については、http://d.hatena.ne.jp/ryamada22/20140329 を参照。

さて、以下のモデルファイルでは、dbin(),dbeta()という2つの分布関数が用いられている。二項分布関数とベータ関数である。

model
{
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dbeta(1,1)
}

1行目は、Bを好む確率がpのときに、n人に訊いたら、nb人が「Bを好む」ことが二項分布に従って定まることを示している。

今、BUGSを回すときにはnとnbとを与えて、適切なpの分布をしミューレーションで示そうとしていることを示している。

他方、2行目はpとしてどんな分布を事前分布を想定しているかを表している。 2つのパラメタが1,1であるようなベータ分布であることを示している。

実行に際しては、このベータ分布からpをランダムに発生し、それを二項分布に与え、nのもとでnbになるかどうかでそのpの乱数を採択するかどうかを決める。

dbeta(1,1)とは

p <- seq(from=0,to=1,length=1000)
db <- dbeta(p,1,1)
plot(p,db)

でわかるように、0から1までの一様分布であるから、dunif(0,1)という分布を指定しても同じである。

実際にやってみる。

model
{
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dunif(0,1)
}

model1.1 <- function(){
   # Model
   nb ~ dbin(p,n)
   
   # Prior
   p ~ dunif(0,1)
}
file.name <- "model1.1.txt"
write.model(model1.1,file.name)

data1 <- list(n=10,nb=4) 
in1 <- list(p=0.5) 
inits <- list(in1) 
param <- c("p") 
file.name <- "model1.1.txt"
model1 <- file.name
# 実行
n.iter <- 300
bugs1.1 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=n.iter, n.burnin=0, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

一様分布を事前分布に指定しても得られる値の出現順は異なるが、分布としては同じ結果が得られている。

post1.1 <- bugs1.1$sims.matrix
(head(post1[,1]))

## [1] 0.6107764 0.3699931 0.4849625 0.4464058 0.4110734 0.3640060

(head(post1.1[,1]))

## [1] 0.2741216 0.2076584 0.2349845 0.2095576 0.6004191 0.2650980

plot(sort(post1[,1]),sort(post1.1[,1]))
abline(0,1,col=2)

このように同じ事前分布なのに、あえてベータ分布を用いたのは、dbeta(1,1)という一様分布でないような分布を使う場合を考えると、二項分布の共役事前分布として、異なる形の事前分布に変えやすいからである。

実行条件

bugs()関数にはデータ、初期値、パラメタ名、モデルファイル、debug、program、bugs.directoryのほかにいくつかの情報を与える。 それぞれを説明する

• n.chains
• MCMCでは初期値と乱数列を使って処理をするが、初期値を変えて何度か実行する方がよいこともあり、その回数を指定する
• n.iter
• 推定対象は多数の値が作る分布として判断するわけだが、その作るべき値の数を示す。ただし、そのうち、初めの方のn.burnin個は捨てるので、実際に推定される分布として使用可能なのはn.iter-n.burnin個になる
• n.burnin
• 上述したように、初めのn.burnin個は捨てる
• 今回のように一番初め値から、推定に用いることが適切な場合にはn.burnin=0でも構わない
• n.thin
• 非常に多くの値を発生させるのは、その方が、初期値の影響が混入することを避けることができるからであるが、あまりにたくさんだと大変なので、n.iterが大きいときには、何個かおきに記録してそれ以外を捨てることも有効である。n.iter=1ならば全部を記録し、n.iter=2ならば一つおきに捨てる、とそういう具合である。

n.chains,n.iter,n.burnin,n.thinの値を変えながら実行して、出力がどのように変わるかを確認してみる。

n.chains=1でn.iter=100,n.burnin=10,n.thin=3としてみる。

bugs1.2 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=100, n.burnin=10, n.thin=3, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs1.2)

## Inference for Bugs model at "model1.1.txt", fit using OpenBUGS,
##  1 chains, each with 100 iterations (first 10 discarded), n.thin = 3
##  n.sims = 30 iterations saved
##          mean  sd 2.5% 25% 50% 75% 97.5%
## p         0.4 0.1  0.2 0.4 0.4 0.5   0.7
## deviance  3.7 1.5  2.8 2.8 3.2 3.8   7.6
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = 4.8
## DIC is an estimate of expected predictive error (lower deviance is better).

n.iter-n.burnin = n.thin*bugs1$n.simsとなっていることがわかる。

次にn.chains=3、n.thin=2でやってみる。初期値をn.chains通り、指定する必要がある。

in1 <- list(p=0.5)
in2 <- list(p=0.3)
in3 <- list(p=0.7)
inits <- list(in1,in2,in3)
inits

## [[1]]
## [[1]]$p
## [1] 0.5
## 
## 
## [[2]]
## [[2]]$p
## [1] 0.3
## 
## 
## [[3]]
## [[3]]$p
## [1] 0.7

このように初期値を3通りリストとして準備している。

さて、実行してみる。

bugs1.3 <- bugs(data1, inits, param, model.file=model1, 
n.chains=3, n.iter=100, n.burnin=10, n.thin=2, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs1.3)

## Inference for Bugs model at "model1.1.txt", fit using OpenBUGS,
##  3 chains, each with 100 iterations (first 10 discarded), n.thin = 2
##  n.sims = 135 iterations saved
##          mean  sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## p         0.4 0.1  0.2 0.3 0.4 0.5   0.7    1    59
## deviance  3.7 1.2  2.8 2.9 3.3 4.1   7.0    1   140
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = 4.8
## DIC is an estimate of expected predictive error (lower deviance is better).

各chainについてn.iter回をn.burnin込みで実施していることが記載されている。

正規分布の平均と分散とを推定する

比較的ふんだんに発現していると思われる遺伝子mRNAをある検出手法で測定したときの値の分布について推定することを考える。 検出手法についてもあまりよくわかっていないので、平均がいくつで分散がいくつになるか、さっぱり予想がつかないが、正規分布に従うだろうことは予想してもよいとする。

このような場合のモデルは以下のようになる。 “model2.txt”というファイルで保存しておくこととする。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu, S)
    }
    mu ~ dnorm(0, 0.0001)
    S ~ dgamma(0.0001,0.0001)
}

model2 <- function(){
   # Model
   for (i in 1:N) {
        y[i] ~ dnorm(mu, tau)
    }
   
   # Prior
    mu ~ dnorm(0, 0.0001)
    tau ~ dgamma(0.0001,0.0001)
   
}
file.name <- "model2.txt"
write.model(model2,file.name)

このモデル設定では、N個の観察値を一つのmu,tauのセットに基づいて正規分布に従うことを示している。

muには、第1引数(v1 平均)0、第2引数(v2 tau) 0.0001の正規分布を事前分布として想定している。 このv2が小さいということは、分散が大きいということであるので、この事前分布は、0を中心として、非常に裾が長い分布である。 要するに、いろいろな値をほぼ同確率で想定していることになる。

v2の値を変化させて、どのような事前分布になるのかを表示してみる。 v2の値が小さくなるほど平坦な分布となることを、実数と対数とで示す。

横軸の中央付近にピークがあることに着目するのではなく、そこにピークはあるものの、それはごく限られた幅に納まり、それ以外の広い範囲で平坦になっていることが、この分布を採用した理由である。

x <- seq(from=-100,to=100,length=1000)
v2s <- 10^(0:(-5))
d <- matrix(0,length(x),length(v2s))
for(i in 1:length(v2s)){
  v2 <- v2s[i]
  d[,i] <- dnorm(x,0,1/sqrt(v2))
}
par(mfcol=c(1,2))
matplot(x,d,type="l")
matplot(x,log(d),type="l")

par(mfcol=c(1,1))

また、標準偏差の逆数としのtauの事前分布はガンマ分布を用いている。 ガンマ分布が選ばれているのは、tauには正の値をとりたいからである。 そのような理由でガンマ分布を採用したうえで、次にその形を決めたい。

dgamma(a,b)で指定されるガンマ分布は、平均が\[\frac{a}{b}\]分散が\[\frac{a}{b^2}\]であるが、事前分布としてはさまざまな値をとりたい、言い換えると分散を大きくしたい。 ガンマ分布の分散は\[\frac{a}{b^2}\]であるので、a=bとして小さい値を与えることで、平均は1で分散の大きい分布が得られる。 これが、tauの事前分布の選び方である。

事前分布としての指定に際してはa,bに等しく、小さい値を与えているが、これは、平均\(\frac{a}{b}=\frac{a}{a}=1\)が1で分散が大きい分布である。

a,b=aの値を色々に変えてガンマ分布の様子をプロットしてみる。

x <- seq(from=0,to=20,length=1000)
as <- 10^(0:(-7))
d <- matrix(0,length(x),length(as))
for(i in 1:length(as)){
  d[,i] <- dgamma(x,as[i],as[i])
}
par(mfcol=c(1,2))
matplot(x,d,type="l")
matplot(x,log(d),type="l")

par(mfcol=c(1,1))

0付近にピークが残るも、aが小さくなるにつれ、平坦な分布となることがわかる。 ではデータを作成して推定してみる。

N <- 100
mu <- 5
sd <- 2
y <- rnorm(N,mu,sd)
data1 <- list(y=y,N=N)
in1 <- list(mu=0,tau=1) 
inits <- list(in1) 
param <- c("mu","tau") 
file.name <- "model2.txt"
model2 <- file.name
n.iter <- 1000
bugs2 <- bugs(data1, inits, param, model.file=model2, 
n.chains=1, n.iter=n.iter, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323") 

## Sampling has been started ...

print(bugs2)

## Inference for Bugs model at "model2.txt", fit using OpenBUGS,
##  1 chains, each with 1000 iterations (first 100 discarded)
##  n.sims = 900 iterations saved
##           mean  sd  2.5%   25%   50%   75% 97.5%
## mu         4.7 0.2   4.3   4.6   4.7   4.8   5.1
## tau        0.3 0.0   0.2   0.3   0.3   0.3   0.4
## deviance 411.9 1.9 410.0 410.5 411.4 412.7 417.0
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 2.0 and DIC = 413.9
## DIC is an estimate of expected predictive error (lower deviance is better).

ここで推定されたのは\[\tau=\frac{1}{sd^2}\]であるので、データ生成時に与えたsd=2に近い値が推定されたかを確認するために計算しておく。 また、サンプル平均、サンプルに基づくsdも計算しておく。

post2 <- bugs2$sims.matrix
mean(post2[,1])

## [1] 4.712235

mean(1/sqrt(post2[,2]))

## [1] 1.894946

# サンプル平均とサンプルからの標準偏差
mean(y)

## [1] 4.712201

sd(y)

## [1] 1.888164

その上で分布を示す。

赤は真値、緑はサンプルから算出した平均とsdである。

par(mfcol=c(1,2))
hist(post2[,1],freq=FALSE,xlab="mu",main="Posterior 
distribution of mu")
abline(v=mu,col=2)
abline(v=mean(y),col=3)
hist(1/sqrt(post2[,2]),freq=FALSE,xlab="sd",main="Posterior 
distribution of sd")
abline(v=sd,col=2)
abline(v=sd(y),col=3)

par(mfcol=c(1,1))

線形回帰をやってみる

\[ \begin{equation} y = beta0 + beta1 \times x1 + beta2 \times x2 + \epsilon\times z\\ z ~ N(0,1) \end{equation} \] という線形回帰をしてみよう。冒頭の「ひとまず回す」の例では、説明変数が1つの線形回帰だったが、今回は2つにしてみる。

具体性を持たせるためにx1を年齢、x2を体重、yを最高血圧として考えよう。 値とその分布の妥当性はここでは問わ無いことにする。

まずはデータを作る。

N <- 500
# 年齢
x1 <- sample(20:100,N,replace=TRUE)
# 体重
x2 <- rnorm(N,65,10)

beta0 <- 30
beta1 <- 0.5
beta2 <- 1
epsilon <- 10

y <- beta0 + beta1*x1+beta2*x2+epsilon * rnorm(N)
hist(y)

# yの値で色を付けて、散布図にする
st.y <- (y-min(y))/(max(y)-min(y))
col <- rgb(st.y,1-st.y,1)
plot(x1,x2,pch=20,col=col,xlab="age",ylab="weight",main="systolic pressure")

モデルファイルを作り、“model3.txt”として保存する。

項の数が増えたが、N人それぞれについて、個人の年齢x1[i]と体重x2[i]に応じて、最高血圧の予想値mu[i]を求めその値を平均とした正規分布に従うことを示してる。 パラメタbeta0,beta1,beta2,epsilonに事前予想が立たないと仮定すれば、 beta0,beta1,beta2は0を中心とした分散の大きな正規分布を事前分布とし、 乱雑項も分散が不明な正規分布とするのがよいので、そのような分布が指定してある。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], tau)
        mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
    }
    beta0 ~ dnorm(0,0.0001)
    beta1 ~ dnorm(0,0.0001)
    beta2 ~ dnorm(0,0.0001)
    tau ~ dgamma(0.0001,0.0001)
}

model3 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], tau)
        mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
    }
    beta0 ~ dnorm(0,0.0001)
    beta1 ~ dnorm(0,0.0001)
    beta2 ~ dnorm(0,0.0001)
    tau ~ dgamma(0.0001,0.0001)
}
file.name <- "model3.txt"
write.model(model3,file.name)

さて、実行する。

data1 <- list(y=y,x1=x1,x2=x2,N=N) 
in1 <- list(beta0=0,beta1=0,beta2=0,tau=1) 
inits <- list(in1) 
param <- c("beta0","beta1","beta2","tau") 
file.name <- "model3.txt"
model1 <- file.name
# 実行
bugs3 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs3)

## Inference for Bugs model at "model3.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##            mean  sd   2.5%    25%    50%    75%  97.5%
## beta0      31.7 3.4   24.7   29.2   32.3   34.2   37.4
## beta1       0.5 0.0    0.4    0.5    0.5    0.5    0.5
## beta2       1.0 0.0    0.9    0.9    1.0    1.0    1.1
## tau         0.0 0.0    0.0    0.0    0.0    0.0    0.0
## deviance 3767.8 2.8 3764.2 3765.7 3767.1 3769.3 3775.2
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 4.1 and DIC = 3772.0
## DIC is an estimate of expected predictive error (lower deviance is better).

post3 <- bugs3$sims.matrix
apply(post3,2,quantile,c(0.1,0.5,0.9))

##        beta0     beta1     beta2         tau deviance
## 10% 26.79475 0.4641732 0.9222129 0.008414926 3764.901
## 50% 32.32948 0.4895096 0.9736940 0.009081334 3767.094
## 90% 35.75325 0.5160894 1.0429780 0.009878016 3771.762

beta0,beta1,beta2が良い推定値になっていることがわかる。

招待選手と一般参加選手の2群～混合分布～

今、全部でN人がエントリーしている100m走大会があるとする。 このうち、pi[1]の割合で招待選手(群 g=1)、pi[2]の割合で一般参加選手(群 g=2)であるとする。

招待選手のタイムは平均mm[1]標準偏差ss[1]、一般参加選手はmm[2],ss[2]の正規分布に従うとする。 ただし、この招待選手というのは、100m走の速さを基準に招待されたわけではなく、マラソンが得意であることで招待されたという。 招待選手はアスリートなので、一般参加選手よりも100m走のパフォーマンスがよいようにも思えるし、大差ないかもしれないし、かえって遅いかもしれない。 pi[1],pi[2]の事前確率は0-1の範囲で全く不明だという。

データを作成してみる。

N <- 1000
pi <- c(0.2,0.8)
n1 <- N*pi[1]
n2 <- N-n1

mm <- c(12,16)
ss <- c(1,2)
# 100m走のタイム
y1 <- rnorm(n1,mm[1],ss[1])
y2 <- rnorm(n2,mm[2],ss[2])
y12 <- c(y1,y2)
# グループ
g12 <- c(rep(1,n1),rep(2,n2))
# 選手番号をランダムにつける
ord <- sample(1:N)
# 選手番号順に並べ替える。これが観察データ
y <- y12[ord]
g <- g12[ord]

h <- hist(y,plot=FALSE)
hist(y1,breaks=h$breaks,density=20,col=2,ylim=c(0,max(h$counts)))
hist(y2,breaks=h$breaks,density=17,col=3,add=TRUE)

ではモデルを立ててみる。

各選手のグループは確率pにより1か2に決まる。 そのpは一様分布を事前確率とする。 各選手のタイムはグループ1か2かによって平均、標準偏差を変える必要があるが、それらに基づく正規分布に従う。 2つのグループのタイムを決める平均と標準偏差はm1,m2,s1,s2とするが、m1,m2は9秒から30秒の範囲で一様分布を想定しよう。s1,s2は全く不明として、正の値を分散の大きいガンマ分布で与えることとする。

model
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau[i])
        m[i] <- mm[g[i]]
        tau[i] <- tt[g[i]]
        g[i] ~ dcat(pi[])
    }
    pi[1:2] ~ ddirch(alpha[])
    mm[1] ~ dunif(9, 30)
    mm[2] ~ dunif(9, 30)
    tt[1] ~ dgamma(1.00000E-04, 1.00000E-04)
    tt[2] ~ dgamma(1.00000E-04, 1.00000E-04)
}

model4 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau[i])
        m[i] <- mm[g[i]]
        tau[i] <- tt[g[i]]
        g[i] ~ dcat(pi[])
    }
    pi[1:2] ~ ddirch(alpha[])
    mm[1] ~ dunif(9,30)
    mm[2] ~ dunif(9,30)
    tt[1] ~ dgamma(0.0001,0.0001)
    tt[2] ~ dgamma(0.0001,0.0001)
}
file.name <- "model4.txt"
write.model(model4,file.name)

さて実行しよう。

alpha <- c(1,1)
data1 <- list(y=sort(y),g=c(1,rep(NA,N-2),2),alpha=alpha,N=N) 
in1 <- list(pi=c(0.5,0.5),mm=c(15,15),tt=c(1,1)) 
inits <- list(in1) 
param <- c("pi","mm","tt") 
file.name <- "model4.txt"
model1 <- file.name
# 実行
bugs4 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Sampling has been started ...

print(bugs4)

## Inference for Bugs model at "model4.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##            mean   sd   2.5%    25%    50%    75%  97.5%
## pi[1]       0.6  0.1    0.5    0.6    0.6    0.7    0.7
## pi[2]       0.4  0.1    0.3    0.3    0.4    0.4    0.5
## mm[1]      16.6  0.2   16.3   16.5   16.6   16.8   17.0
## mm[2]      13.1  0.3   12.6   12.8   13.0   13.2   13.6
## tt[1]       0.3  0.0    0.3    0.3    0.3    0.4    0.4
## tt[2]       0.4  0.1    0.3    0.3    0.4    0.4    0.5
## deviance 3862.0 65.2 3755.1 3815.7 3855.3 3902.3 4014.2
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 2128.6 and DIC = 5990.6
## DIC is an estimate of expected predictive error (lower deviance is better).

2群の割合と2群の平均がよい感じに推定されていることがわかる。 推定されたtauから標準偏差を計算し、データ生成に用いたそれがよく推定されていることも確認できる。

post4 <- bugs4$sims.matrix
apply(post4,2,quantile,c(0.1,0.5,0.9))

##         pi[1]     pi[2]    mm[1]    mm[2]     tt[1]     tt[2] deviance
## 10% 0.5489903 0.3092796 16.42602 12.71375 0.3042551 0.3142937 3784.277
## 50% 0.6310976 0.3689024 16.63914 13.02444 0.3459111 0.3938328 3855.280
## 90% 0.6907205 0.4510097 16.86821 13.43095 0.3949724 0.4934491 3949.119

mean(1/sqrt(post4[,5]))

## [1] 1.701282

mean(1/sqrt(post4[,6]))

## [1] 1.598109

n次線形回帰をしてみる

\[ y = \sum_{i=0}^n b_i x^i + \epsilon z \] (ただしzは標準正規分布に従う乱雑項) のようなデータから、n+1個のbの値を推定してみる。

まずはデータ作成。ただし、三角関数を使って作成する。

N <- 100
x <- sort(rnorm(N))
y <- 3*sin(x)
epsilon <- 0.5
z <- rnorm(N)
y <- y + epsilon * z
plot(x,y,pch=20)

モデルファイルを作る。

model5 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau)
        m[i] <- b[1]+b[2]*x[i] + b[3]*x[i]*x[i] + b[4] * pow(x[i],3)
    }
    for(i in 1:(n+1)){
      b[i] ~ dnorm(0,0.0001)
    }
    tau ~ dgamma(0.0001,0.0001)
}
file.name <- "model5.txt"
write.model(model5,file.name)

実行する。

n <- 3
data1 <- list(y=y,x=x,N=N,n=n) 
in1 <- list(b=rep(0,n+1),tau=1) 
inits <- list(in1) 
param <- c("b","tau") 
file.name <- "model5.txt"
model1 <- file.name
# 実行
bugs5 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Warning in readChar(modelFile, 10^3): can only read in bytes in a non-UTF-8
## MBCS locale

## model is syntactically correct
## data loaded
## model compiled
## Initializing chain 1: 
## model is initialized
## model is already initialized

## Sampling has been started ...

## 100 updates took 0 s
## deviance set
## monitor set for variable 'b'
## monitor set for variable 'tau'
## monitor set for variable 'deviance'
## 1000 updates took 0 s

print(bugs5)

## Inference for Bugs model at "model5.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##           mean  sd  2.5%   25%   50%   75% 97.5%
## b[1]       0.0 0.1  -0.1   0.0   0.0   0.0   0.1
## b[2]       3.0 0.1   2.8   2.9   3.0   3.0   3.1
## b[3]       0.0 0.0  -0.1   0.0   0.0   0.0   0.1
## b[4]      -0.4 0.0  -0.5  -0.4  -0.4  -0.4  -0.3
## tau        4.9 0.7   3.6   4.4   4.8   5.3   6.3
## deviance 126.3 3.1 122.1 124.0 125.7 127.9 133.8
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 4.9 and DIC = 131.1
## DIC is an estimate of expected predictive error (lower deviance is better).

推定された係数を用いて、推定線を描いてみる。

post5 <- bugs5$sims.matrix
out5.coef <- apply(post5,2,mean)
plot(x,y,pch=20)
x. <- seq(from=min(x)-0.1*(max(x)-min(x)),to=max(x)+0.1*(max(x)-min(x)),length=100)
pred.y <- rep(out5.coef[1],length(x.))
for(i in 2:(n+1)){
  pred.y <- pred.y + out5.coef[i] * x.^(i-1)
}
points(x.,pred.y,col=2,type="l")

カーネル分布推定

N個の値を観測してそこに密度分布を推定するとする。 正規分布をカーネル関数としてカーネル法でやるなら、個々の観測に、ある分散の正規分布を与え、その混合分布として求めることとなる。

混合分布なら、すでに、「100m走大会～招待選手と一般参加選手」の例でやったので、それを少し改変すればよい。

100m走大会では2群だったが、今度は、データ数Nのすべてを個別の群として扱うので群数はNである。

それぞれの群の寄与度は1/Nで固定である。

推定したいのは、個々の値を中心とした正規分布のその分散をいくつにするかである。

すべての値について分散の値を共通にするのであれば、推定パラメタは一つだけである。

それを書けば以下の通り。

yは観測値のベクトル。 mmもyと同じベクトルであるが、これはBUGS中に採取されるのではなくて、決まった値。 tauは推定パラメタである。 gはN個の群のラベルである。 モデルファイルの中では、個々の値が、1,…,Nのいずれかの群に属するかを等確率で採取し、取られた群g[i]に応じて、正規分布の平均値m[i]が決まり、tauに応じて値が採取している。

piは1,…,Nを等確率で選ぶために与えるベクトルである。

model6 <- function()
{
    for (i in 1:N) {
        y[i] ~ dnorm(m[i], tau)
        m[i] <- mm[g[i]]
        g[i] ~ dcat(pi[])
    }

    tau ~ dgamma(0.0001,0.0001)
}
file.name <- "model6.txt"
write.model(model6,file.name)

データを作る。

y <- c(rnorm(100),runif(100,3,8),rexp(200,4)+10)
y <- sort(y)
plot(density(y))

N <- length(y)
alpha <- rep(100,N)
data1 <- list(y=y,mm=y,N=N,g=1:N,pi=rep(1/N,N)) 
in1 <- list(tau=1) 
inits <- list(in1) 
param <- c("tau") 
file.name <- "model6.txt"
model1 <- file.name
# 実行
bugs6 <- bugs(data1, inits, param, model.file=model1, 
n.chains=1, n.iter=1100, n.burnin=100, n.thin=1, debug=TRUE, program="OpenBUGS", bugs.directory="C:/Program Files(x86)/OpenBUGS/OpenBUGS323")

## Warning in readChar(modelFile, 10^3): can only read in bytes in a non-UTF-8
## MBCS locale

## model is syntactically correct
## data loaded
## model compiled
## Initializing chain 1: 
## model is initialized
## model is already initialized

## Sampling has been started ...

## 100 updates took 0 s
## deviance set
## monitor set for variable 'tau'
## monitor set for variable 'deviance'
## 1000 updates took 0 s

print(bugs6)

## Inference for Bugs model at "model6.txt", fit using OpenBUGS,
##  1 chains, each with 1100 iterations (first 100 discarded)
##  n.sims = 1000 iterations saved
##               mean       sd      2.5%       25%       50%       75%
## tau      2000861.2 139442.3 1737639.4 1902228.9 2000254.2 2094249.0
## deviance    -274.3     28.0    -325.8    -293.6    -275.2    -255.1
##              97.5%
## tau      2270276.1
## deviance    -218.9
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.0 and DIC = -273.4
## DIC is an estimate of expected predictive error (lower deviance is better).

推定された係数を用いて、推定線を描いてみる。density()の結果と似ている。

post6 <- bugs6$sims.matrix
out6.coef <- apply(post6,2,mean)
x <- seq(from=min(y)-(max(y)-min(y))*0.3,to=max(y)+(max(y)-min(y))*0.3,length=100)
pred.d <- rep(0,length(x))
for(i in 1:length(y)){
  pred.d <- pred.d + dnorm(x-y[i],1/sqrt(out6.coef[1]))/length(y)
}
par(mfcol=c(1,2))
plot(x,pred.d,type="l")
plot(density(y))

par(mfcol=c(1,1))

Permalink | コメント