Wolfeyes Bioinformatics beta

本書の元となった文章のことは前々から知っていたし，pdfファイルもPCにダウンロードしていたのだけれども，結局その文章を読むことはなかった．書店で本書が並んでいるのを見た時に「関心があるだけフリをしているだけではダメだ」と思って本書を購入し，まとまった時間を作って一気に読んだ．

本書はタイトルの通り，原発事故後の日本に暮らす私達が放射線とどう付き合っていけばよいかについて，科学的な知見を紹介する内容となっている．本の中で説明されるのは，放射線の基礎的な性質や放射線が生物にどのように影響を与えるのかといった仕組み，そして国際的に認められた機関が放射線と人間の関係においてどのような基準を作っているのかといった，あくまで判断材料としての基礎的な情報だ．そのため，今の生活が安心なのかそれとも危険なのかといった結論は，直接的には書かれていない．あくまで判断をするのは私達自身だという前提のもとで，できるだけ専門性は排除して初学者にも分かりやすい形で本書は構成されている．

本書を書かれた田崎晴明氏は，個人的には統計力学の先生というイメージが強い．分野が違うのに何故と最初は思っていたのだが，本書を読むにつれ，田崎氏の科学者として正しい情報を普及させたいという真摯な態度が伝わってきて，その理由が納得できた．非専門家でありながらも，放射線の説明は見事としか言いようがない．問題の本質をしっかりと捉え，噛み砕いた形でありながらも正確さを損なっていない．特に，放射線は化学反応とは別の次元の物理現象なのだという前提は，まさに放射線を取り巻く全ての問題のを理解する上で知っておくべき一番重要なことだろう．他にも，不正確な情報の中からいかに上限を見積もって判断材料にするかといったことも述べられており，科学的な考え方のエッセンスが詰まっている．そういった面でも非常に参考になる本であり，啓蒙書としての見本となるような文章だと感じた．

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

例 6.4 ランダム・ウォークのパラメータ設定によるサンプリングへの影響を，正規分布の乱数生成を例に確かめる

摂動が $\varepsilon_t \sim \mathcal{U}_{[-\delta,\delta]}$ となるようなランダム・ウォークにおいて，正規分布 $\mathcal{N}(0,1)$ をサンプリングするときのメトロポリス・ヘイスティング・アルゴリズムの受理確率は

$\rho(x^{(t)},y_t) = \min \big[1, \exp{\frac{x^{(t)2}-y_t^2}{2} } \big]$

で表される．

今回は，一様候補のパラメータを $\delta = 0.1, 1, 10$ と変更して，メトロポリス・ヘイスティング・アルゴリズムによるローカルな探索がどういう挙動を示すのかを調べる．

par(mfcol=c(3,3))
for(d in c(0.1,1,10)){
  Niter <- 2500
  x <- 0
  for(i in 2:Niter){
    y <- x[i-1] + runif(1, min=-d, max=d)
    if (runif(1) < min(1, exp(-(1/2)*(y^2 - x[i-1]^2)) )){
      x[i] <- y
    }else{
      x[i] <- x[i-1]
    }
  }
  plot(1:Niter, x, type="l", xlab="Iterations", ylab="",main=paste("d = ",d))
  hist(x, nclass=150, freq=F, xlim=c(-4,4), main="")
  curve(dnorm(x), lwd=2, add=T, col="red")
  acf(x)
}

それぞれのパラメータにおける系列の変動，分布，自己相関プロットを示したのが以下の図である．それぞれのサンプリング結果を見比べてみると，

一様候補のパラメータが小さい場合（）
- ランダム・ウォークで移動する幅が狭い
- 分布が偏っている
- 自己相関が高い（前の値からあまり変わらない）
一様候補のパラメータが大きい場合（）
- ランダム・ウォークで移動する幅が大きい
- これも同様に分布が安定しない

という結果となり，ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムにおいては摂動（ $\varepsilon_t$ ）のパラメータ設定が重要であることがわかる．

例 6.10 ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムの受理率を調べる

メトロポリス・ヘイスティング・アルゴリズムにおけるアルゴリズム効率化の調整は，本書の2章で扱った受理・棄却アルゴリズムの時の最大受理率に基づいた調整の場合とは異なる．まず，上の例 6.4で示した正規分布からのサンプリングにおける受理率を見てみる．パラメータごとの受理率は以下のようになった（Rのコードは省略）．

$\delta = 0.1$ ： $r = 0.9792$
$\delta = 1$ ： $r = 0.8132$
$\delta = 10$ ： $r = 0.1716$

この結果から，受理率が高ければサンプリングの結果も良くなるわけではないことがわかる．このことから，ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムにおいては，受理率を高めるようにアルゴリズムを調節してはいけないと判断できる．このアルゴリズムにおいて，どのような受理率を取ればよいかといった疑問に関しては，以下の練習問題6.14のc.で詳しく見ていく．

練習問題 6.14 ランダム・ウォーク候補がサンプリング結果に与える影響について，複数の候補分布を用いてより詳細に調べる

a. P.205の図6.7を再現する

図6.7の再現と $\delta$ の解釈は，上の例6.4の解答を参考．

b. コーシー候補とラプラス候補を検討する

今度は，一様分布の代わりにコーシー分布やラプラス分布を候補分布として使った場合の，ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムの挙動を調べる．

コーシー候補

ランダム・ウォーク候補は $\varepsilon_t \sim \mathcal{C}(x_0,\gamma)$ とし，コーシー分布のパラメータは $x_0 = 0$ ， $\gamma = 0.01,\ 1,\ 10$ を考える．

# Cauchy candidate
par(mfcol=c(3,3))
for(g in c(0.01,1,10)){
  Niter <- 2500
  x <- 0
  for(i in 2:Niter){
    y <- x[i-1] + rcauchy(1, location=0, scale=g)
    if (runif(1) < min(1, exp(-(1/2)*(y^2 - x[i-1]^2)) )){
      x[i] <- y
    }else{
      x[i] <- x[i-1]
    }
  }
  plot(1:Niter, x, type="l", xlab="Iterations", ylab="",main=paste("Cauchy candidate,","gamma=",g))
  hist(x, nclass=150, freq=F, xlim=c(-4,4), main="")
  curve(dnorm(x), lwd=2, add=T, col="red")
  acf(x)
}

ラプラス候補

ランダム・ウォーク候補は $\varepsilon_t \sim \mathcal{L}(\mu,b)$ とし，ラプラス分布のパラメータは $\mu=0$ ， $b = 0.1,\ 1,\ 10$ を考える．

# Laplace candidate
par(mfcol=c(3,3))
for(a in c(0.1,1,10)){
  Niter <- 2500
  x <- 0
  for(i in 2:Niter){
    y <- x[i-1] + ifelse(runif(1)>0.5, 1, -1)*rexp(1)/a
    if (runif(1) < min(1, exp(-(1/2)*(y^2 - x[i-1]^2)) )){
      x[i] <- y
    }else{
      x[i] <- x[i-1]
    }
  }
  plot(1:Niter, x, type="l", xlab="Iterations", ylab="",main=paste("Laplace candidate,","b=",a))
  hist(x, nclass=150, freq=F, xlim=c(-4,4), main="")
  curve(dnorm(x), lwd=2, add=T, col="red")
  acf(x)
}

以上の結果から，コーシー候補やラプラス候補においても，ランダム・ウォークの値の取り方によって適切なサンプリングができるかどうかが決まってくることがわかる．また，どちらも自己相関が低い値の時にサンプリング結果が一番良いことが見て取れる．これは例6.4の一様候補の場合とは違う結果となったが，恐らくコーシー分布やラプラス分布はともに一峰型の正規分布のような形であることが効いていて，ランダム・ウォークの移動が適度に調節されているからだろう．

c. パラメータの調節により受理確率を0.25に近づける

さて，例6.10で述べたように，ランダム・ウォークを使ったメトロポリス・ヘイスティング・アルゴリズムにおいては，受理率の最大化がアルゴリズムの効率化に繋がらないことがわかった．受理率はどれくらいを目安にすればよいかという疑問に関しては，本書では「高次元では0.25，低次元では0.5程度」というRoberts et al.(2007)の提案を引用している．そこで，先ほどの3つのランダム・ウォーク候補について受理率を0.25に近づけるために，パラメータの値を0から10まで0.1刻みで設定したときの，一様候補・コーシー候補・ラプラス候補を用いたアルゴリズムにおける受理率の変化を以下の図でプロットした．

赤の点線は受理率0.25を表している．この図を見ると，一様候補とコーシー候補はパラメータの値が大きくなるにつれて受理率も小さくなり，逆にラプラス候補の場合ではパラメータの値が大きくなるにつれて受理率は大きくなることがわかる．いずれの場合においても，受理率を0.25に近づけることは可能である．

MacPortsでJAGSを入れる場合はちょっと注意

MacでJAGSをインストールする場合，

SourceForgeからパッケージをダウンロードしてインストールする方法
MacPortsを使ってインストールする方法

の2通りがある．どちらで入れてもJAGS自体に違いはないのだが，Rのrjagsパッケージをインストールする際に少し注意しなければいけないポイントがある．1つ目の方法でパッケージからインストールした場合だと，いつも通りに

# パッケージからインストールした場合
> install.packages("rjags")
> library(rjags)

とすれば良いのだが，MacPortsでインストールした場合では，

# MacPortsからインストールした場合
> install.packages("rjags", type="source", configure.args = '--with-jags-lib=/opt/local/lib --with-jags-include=/opt/local/include/JAGS --with-jags-modules=/opt/local/lib/JAGS/modules-3')
> library(rjags)

のように，パッケージのインストールの際に色々とオプションを付ける必要がある．というのも，パッケージからインストールした場合では/usr/local/以下にインストールされるのに対し，MacPortsからインストールした場合では/opt/local/以下にインストールされるからだ．Rでrjagsパッケージをインストールする際には，どうやら/usr/local/以下のモジュールや共有ライブラリを参照するようなので，MacPortsで入れた場合にはそもそもディレクトリが違うためファイルが見つからないといったエラーが出る．

他の解決方法

ちなみに，/usr/local/lib以下にシンボリックリンクを貼って強引に解決する方法もある．自分は最初こちらで対処していたのだが，長期的に考えた場合にJAGSのバージョンアップによってdylibの名前が変わる可能性などもあるので，上のようにパッケージのインストールの段階で対処する方が良いだろう．

$ cd /usr/local/lib
$ sudo ln -s /opt/local/lib/libjags.3.dylib .
$ sudo ln -s /opt/local/lib/JAGS .

参考：MacPortsの場合でinstall.package(“rjags”)で入れてしまったときのrjags読み込み時のエラー

> library(rjags)
 要求されたパッケージ coda をロード中です
 要求されたパッケージ lattice をロード中です
Error :  .onLoadはloadNamespace()（rjagsに対する）の中で失敗しました、詳細は:
 call: dyn.load(file, DLLpath = DLLpath, ...)
 error:  共有ライブラリ '/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rjags/libs/x86_64/rjags.so' を読み込めません：
  dlopen(/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rjags/libs/x86_64/rjags.so, 10): Library not loaded: /usr/local/lib/libjags.3.dylib
  Referenced from: /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rjags/libs/x86_64/rjags.so
  Reason: image not found
 エラー：  '‘rjags’' に対するパッケージもしくは名前空間のロードが失敗しました

実行環境

Mac OS X 10.8.2
R version 2.15.1
JAGS 3.3.0
rjags version 3-7

参考

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

例 6.3 二次モデルの線形回帰におけるパラメータの事後分布から，メトロポリス・ヘイスティング・アルゴリズムを用いてサンプリングする

この例題では，carsのデータセットにおいて線形回帰による二次曲線のあてはめをおこない，最尤法で推定したパラメータによる事後分布から複数のサンプリングを行うことを目標にする．今回使用するデータセットは車の速度とその時の制動距離に関するデータで，Rの入門書でよく例に使われるものである．carsのデータセットに関するHelpのExampleには，本問題と似たようなプロットを描写するコード例が載っている．

この問題では，車のスピード $x$ と制動距離 $y$ の関係について，以下式で表される二次モデルで線形回帰をおこなう．

$y_{ij}=a+bx_i+cx_i^2+\varepsilon_{ij}$

data(cars)

y <- cars$dist
x <- cars$speed
x2 <- x^2
fit <- lm(y ~ x + x2)
s <- seq(4,25,by=0.1)
yhat <- fit$coef[1] + fit$coef[2]*s + fit$coef[3]*s^2

plot(x,y, ylim=c(0,120), xlab="speed", ylab="dist", main="data(cars)")
lines(s, yhat, lwd=4, col="red")

以下の図は，車のスピードをx軸に，制動距離をy軸に取って，carsに含まれる50個のデータをプロットし，オレンジの線で回帰直線を表示している．また，この時の線形回帰で求められた各パラメータの推定値をsummary()関数を用いて表示している．この情報から，それぞれのパラメータの推定値は

$a = 2.47014$
$b = 0.91329$
$c = 0.09996$

ということがわかる．次の練習問題 6.11では，この推定値を元に，事後分布のサンプリングをメトロポリス・ヘイスティング・アルゴリズムで行う．

> summary(fit)

Call:
lm(formula = y ~ x + x2)

Residuals:
    Min      1Q  Median      3Q     Max
-28.720  -9.184  -3.188   4.628  45.152

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.47014   14.81716   0.167    0.868
x            0.91329    2.03422   0.449    0.656
x2           0.09996    0.06597   1.515    0.136

Residual standard error: 15.18 on 47 degrees of freedom
Multiple R-squared: 0.6673,  Adjusted R-squared: 0.6532
F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12

練習問題 6.11 メトロポリス・ヘイスティング・アルゴリズムを用いて事後分布からのサンプリングをする（例6.3の続き）

a. P.202の図6.6を描写する

ここで $\varepsilon_{ij} \sim \mathcal{N}(0,\sigma^2)$ とすると，この時の尤度関数は以下式のように表される．

$p(y_{ij}|x_i,a,b,c,\sigma^2) = \prod_{ij} \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(y_{ij}-a-bx_i-cx_i^2)^2}{2\sigma^2})$

これを整理すると．

$p(y_{ij}|x_i,a,b,c,\sigma^2) \sim \bigg( \frac{1}{\sqrt{\sigma^2}} \bigg)^N \exp{(-\frac{-1}{2\sigma^2} \sum_{ij}(y_{ij}-a-bx_i-cx_{ij})^2 )}$

となる．この尤度関数を $x_i,a,b,c,\sigma^2$ による事後分布とみなす．それぞれのパラメータの分布は，先ほどの線形回帰で求めた推定値と標準誤差を用いる．

$a \sim \mathcal{N}(2.47, (14.8)^2)$
$b \sim \mathcal{N}(0.913, (2.03)^2)$
$c \sim \mathcal{N}(0.100, (1.52)^2)$
$\sigma^{-2} \sim \mathcal{G}(n/2, (n-3)(15.18)^2)$

なお，上の推定値や標準誤差は本書に例示されている値と若干異なっているが，今回は自分の計算環境で計算した推定値を用いることにした．

# 例6.3での線形回帰
data(cars)
y <- cars$dist
x <- cars$speed
x2 <- x^2
fit <- lm(y ~ x + x2)

loglik <- function(a,b,c,sigma){-(n/2)*log(sigma)-sum((y-a-b*x-c*x2)^2)/(2*sigma)}
logq <- function(a,b,c,sigma){
  dnorm(a,mean=a.est[1],sd=a.std,log=T)
  + dnorm(b,mean=b.est[1],sd=b.std,log=T)
  + dnorm(c,mean=c.est[1],sd=c.std,log=T)
  - (n/2)*log(sigma) - (n-3)*(15.18)^2/(2*sigma)
}

Niter <- 4000
n <- length(x)
fit.param <- summary(fit)$coef
a.est <- fit.param[1,1]
a.std <- fit.param[1,2]
b.est <- fit.param[2,1]
b.std <- fit.param[2,2]
c.est <- fit.param[3,1]
c.std <- fit.param[3,2]
sig <- 1/rgamma(1, n/2, rate=(n-3)*(15.18)^2/2)

for(i in 2:Niter){
  candidate <- c(rnorm(1, mean=a.est[1], sd=a.std),
                 rnorm(1, mean=b.est[1], sd=b.std),
                 rnorm(1, mean=c.est[1], sd=c.std),
                 1/rgamma(1, n/2, rate=(n-3)*(15.18)^2/2)
                 )
  fy <- loglik(candidate[1], candidate[2], candidate[3], candidate[4])
  fx <- loglik(a.est[i-1], b.est[i-1], c.est[i-1], sig[i-1])
  qxy <- logq(a.est[i-1], b.est[i-1], c.est[i-1], sig[i-1])
  qyx <- logq(candidate[1], candidate[2], candidate[3], candidate[4])
  rho <- min(1, exp(fy-fx+qxy-qyx))

  if (runif(1)<rho){
    a.est[i] <- candidate[1]
    b.est[i] <- candidate[2]
    c.est[i] <- candidate[3]
    sig[i] <- candidate[4]
  }else{
    a.est[i] <- a.est[i-1]
    b.est[i] <- b.est[i-1]
    c.est[i] <- c.est[i-1]
    sig[i] <- sig[i-1]
  }
}

s <- seq(4,25,by=0.1)
yhat <- fit$coef[1] + fit$coef[2]*s + fit$coef[3]*s^2

plot(x,y, ylim=c(0,120), xlab="speed", ylab="dist", main="data(cars)", type="n")
for(i in (Niter-500):Niter){
  lines(x, a.est[i]+b.est[i]*x+c.est[i]*x2, lwd=2, col="grey")
}
lines(s, yhat, lwd=4, col="orange")
points(x, y, pch=19)

以下の図は，先ほどの線形回帰のプロットに加えて，サンプリングから得られた値を用いて回帰曲線を灰色の曲線として描写したものである．

b. 各パラメータの収束をモニタリングする

先ほどのサンプリング結果を， $a,b,c,\sigma^2$ について，4000回のシミュレーションで生成された系列と，そのときの自己相関プロット，そして系列の分布をプロットした．

par(mfrow=c(4,3))
plot(1:Niter, a.est, type="l", xlab="Iterations", main="a")
acf(a.est)
plot(density(a.est), main="")
plot(1:Niter, b.est, type="l", xlab="Iterations", main="b")
acf(b.est)
plot(density(b.est), main="")
plot(1:Niter, c.est, type="l", xlab="Iterations", main="c")
acf(c.est)
plot(density(c.est), main="")
plot(1:Niter, sig, type="l", xlab="Iterations", main=expression(sigma^2))
acf(sig)
plot(density(sig), main="")

c. 95%信用区間を求める

> quantile(a.est, c(0.025,0.975))
     2.5%     97.5%
-14.90933  18.98784

先日の伊庭さんのニコ生MCMC講義の中で出てきた「N次元の単位球の体積は5次元で最大になる」というのをRで試してみた．

体積

N次元の単位球の体積は $V_n = \frac{\pi^{n/2}}{\Gamma(1+n/2)}$ で表される．

plot(apply(matrix(seq(1,15)),1,function(n){pi^(n/2)/gamma(n/2+1)}),type="b",xlab="N dimension",ylab="V",main="Volume of the unit ball in N dimensions")

たしかに5次元で最大になっている．もはや想像が出来無い…．

表面積

ちなみに，単位球の表面積は7次元で最大になる．（ $A_n = nV_n = \frac{2\pi^{n/2}}{\Gamma(n/2)}$ ）

plot(apply(matrix(seq(1,15)),1,function(n){2*pi^(n/2)/gamma(n/2)}),type="b",xlab="N dimension",ylab="A",main="Surface area of the unit ball in N dimensions")