バイオインフォマティクス特論3

Size: px

Start display at page:

Download "バイオインフォマティクス特論3"

そうすけいなくら
3 years ago
Views:

1 藤博幸

2 平均と標準偏差を推測するの科学者の科学者のコードの改造

3 ベイズ統計で実践モデリング 4.1 平均と標準偏差を推測する第 4 章ガウス分布を使った推論

4 データは n 個の独な観測値 x 1, x 2, x n 例 : ある細胞の遺伝の発現量を温度時間などの条件を全く同じ状態で独に n 回計測ガウス分布に従うと仮定このデータが従うガウス分布の平均と標準偏差を推測しよう注意 :Jags ではガウス分布は平均と精度 ( 精度は標準偏差の逆数 ) の 2 つにパラメータで表される

5 μ σ x 1 x 2 x 3... x n

6 プレート記法によるグラフィカルモデル μ σ x i i data

7 連続値離散値観測変数 ( 確率変数 ) 観測変数依存関係

8 μ σ 尤度仮定からデータ x i は同じ正規分布に従うので n 尤度 =Πdnorm(x i, µ, 1/σ 2 ) i=1 x i i data モデルとしては x i ~ Gaussian(µ, 1/σ 2 ) と表現

9 事前分布 μ σ µ は平均 0 で精度が小さい (= 分散が大きい ) 正規分布に従うとする µ ~ Gaussian(0, 0.001) x i i data σ は上限を 10 に設定して一様分布で表現 σ ~ Uniform(0, 10) 無情報事前分布あるいはそれに近い形で表現

10 x <- seq(-10, 10, 0.01) plot(x, dnorm(x, 0, ),ty='l', ylim=c(0,10^(-10000)))

11 Jags でのモデルの記述 Gaussian.txt μ x i σ i data # Inferring the Mean and Standard Deviation # of a Gaussian model{ # Data Come From A Gaussian for (i in 1:n){ x[i] ~ dnorm(mu,lambda) } # Priors mu ~ dnorm(0,.001) sigma ~ dunif(0,10) lambda <- 1/pow(sigma,2) }

12 Gaussian_jags.R 変数のクリア # clears workspace: rm(list=ls()) Rate_1.txt と Rate_1_jags.R のあるディレクトリへの移動 # sets working directories: setwd("/users/toh/desktop/code/parameterestimation/gaussian") library(r2jags) パッケージの読み込み x <- c(1.1, 1.9, 2.3, 1.8) n <- length(x) 観測変数の設定 data <- list("x", "n") # to be passed on to JAGS

13 μ と σ の初期値の設定 myinits <- list( list(mu = 0, sigma = 1)) # parameters to be monitored: parameters <- c("mu", "sigma") # The following command calls JAGS with specific options. # For a detailed description see the R2jags documentation. samples <- jags(data, inits=myinits, parameters, model.file="gaussian.txt", n.chains=1, n.iter=1000, n.burnin=1, n.thin=1, DIC=T) # Now the values for the monitored parameters are in the "samples" object, # ready for inspection. mu <- samples$bugsoutput$sims.list$mu sigma <- samples$bugsoutput$sims.list$sigma モニターするパラメータの設定 Jags を呼び出して μ と σ をサンプリング

14 traceplot(samples) par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las=1) Nbreaks <- 80 y <- hist(mu, Nbreaks, plot=f) plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="s", lwd=2, lty=1, xlim=c(0,15), ylim=c(0,3), xlab="rate", ylab="posterior Density") par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las=1) Nbreaks <- 80 y <- hist(sigma, Nbreaks, plot=f) plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="s", lwd=2, lty=1, xlim=c(0,8), ylim=c(0,5), xlab="rate", ylab="posterior Density")

15 μ のトレース σ のトレース

16 μ の事後分布 σ の事後分布

17 μ の要約 σ の要約 summary(mu) V1 Min. : st Qu.: Median : Mean : rd Qu.: Max. : summary(sigma) V1 Min. : st Qu.: Median : Mean : rd Qu.: Max. : density(mu)$x[which(density(mu)$y==max(density(mu)$y))] [1] density(sigma)$x[which(density(sigma)$y==max(density(sigma)$y))] [1]

18 ベイズ統計で実践モデリングの科学者第 4 章ガウス分布を使った推論

19 実験スキルのきく異なる 7 の科学者が全員同じ量について測定をう直感的には最初の 2 は適性をいた研究者この量の真の値は 10 をわずかに下回る程度と考えられるこの測定量についての事後分布を求めて測定される量について推論する副次的には 7 の科学者の測定スキルをついて推論する

20 σ i 仮定 (1) 全ての科学者の測定はガウス分布に従う (2) 全員同じ量を測定しているのでそれぞれのガウス分布は同じ平均値を持つ λ i (3) 標準偏差の違いで実験スキルの違いを表現 μ x i i 番のデータプレート表記でのグラフィカルモデル

21 σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 σ 7 λ 1 λ 2 λ 3 λ 4 λ 5 λ 6 λ 7 x 1 x 2 x 3 x 4 x 5 x 6 x 7 μ

22 # The Seven Scientists model{ # Data Come From Gaussians With Common Mean But Different Precisions for (i in 1:n){ x[i] ~ dnorm(mu,lambda[i]) 尤度の要素の設定 } # Priors mu ~ dnorm(0,.001) for (i in 1:n){ 事前分布の設定 lambda[i] ~ dgamma(.001,.001) sigma[i] <- 1/sqrt(lambda[i]) } }

23 事前分布の形 x <- seq(0, 1.0, 0.001) plot(x, dgamma(x, 0.001, 0.001), ty='l') 先の例では標準偏差の事前分布として様分布を使ったが今回は精度の事前分布としてガンマ分布を使

24 SevenScientists_jags.R # clears workspace: rm(list=ls()) 変数のクリア Rate_1.txt と Rate_1_jags.R のあるディレクトリへの移動 # sets working directories: setwd("/users/toh/desttop/code/parameterestimation/gaussian") library(r2jags) パッケージの読み込み x <- c( ,3.570,8.191,9.898,9.603,9.945,10.056) n <- length(x) data <- list("x", "n") # to be passed on to JAGS 観測変数の設定

25 myinits <- list( list(mu = 0, lambda = rep(1,n))) # parameters to be monitored: parameters <- c("mu", "sigma") μ と λ の初期値の設定モニターするパラメータの設定 Jags を呼び出して μ と σ をサンプリング # The following command calls JAGS with specific options. # For a detailed description see the R2jags documentation. samples <- jags(data, inits=myinits, parameters, model.file ="SevenScientists.txt", n.chains=1, n.iter=1000, n.burnin=1, n.thin=1, DIC=T) # Now the values for the monitored parameters are in the "samples" object, # ready for inspection.

26 MCMC のサンプリングのトレース traceplot(samples) σ は 7 分出てくる

27 μ の事後分布のヒストグラム par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las=1) Nbreaks <- 80 y <- hist(mu, Nbreaks, plot=f) plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="s", lwd=2, lty=1, xlim=c(0,15), ylim=c(0,3), xlab="rate", ylab="posterior Density")

28 μ の要約 summary(mu) V1 Min. : st Qu.: Median : Mean : rd Qu.: Max. :10.658

29 σ の事後分布のヒストグラム Nbreaks <- 80 hist(sigma[,1], Nbreaks) hist(sigma[,7], Nbreaks)

30 σ の要約 summary(sigma) V1 V2 V3 V4 V5 V6 V7 Min. : Min. : Min. : Min. : Min. : Min. : Min. : st Qu.: st Qu.: st Qu.: st Qu.: st Qu.: st Qu.: st Qu.: Median : Median : Median : Median : Median : Median : Median : Mean : Mean : Mean : Mean : Mean : Mean : Mean : rd Qu.: rd Qu.: rd Qu.: rd Qu.: rd Qu.: rd Qu.: rd Qu.: Max. : Max. : Max. : Max. : Max. : Max. : Max. : とだけでなく 3 もスキルに問題がありそう

31 モード ( 最頻値 ) を出 for (i in 1:n) { print(density(sigma[,i])$x[which(density(sigma[,i])$y==max(density(sigma[,i])$y))]) } [1] [1] [1] [1] [1] [1] [1] 結果が平均やメジアン ( 中央値 ) ほど安定していない MCMC のステップ数 ( サンプル数 ) が少ないせいか?

32 ステップ数 ( サンプル数 ) が少ないバーンインやシンニングが考慮されていない正確な不変分布 ( 事後分布 ) が得られていないのでは?

33 burn in を100 ß--- n.burnin thinningを50 ß--- n.thin MCMCの全ステップ数を ß--- n.iter chainを2にして初期値を与える ß-- n.chains SevenScientistsv2_jags.R (=SevenScientists_jags.R をコピーしたもの ) を書き換える 2 回目の Rate_1_jags.R を参考にできるモデルは同じなので SevenScientists.txt は変更しなくて良い

34 (1) 実部分を書き換える samples <- jags(data, inits=myinits, parameters, model.file ="SevenScientists.txt", n.chains=1, n.iter=1000, n.burnin=1, n.thin=1, DIC=T) samples <- jags(data, inits=myinits, parameters, model.file ="SevenScientists.txt", n.chains=2, n.iter=100000, n.burnin=100, n.thin=50, DIC=T)

35 (2) 2 つの連鎖の初期値を与える myinits <- list( list(mu = 0, lambda = rep(1,n))) myinits <- list( list(mu = 0, lambda = rep(1,n)), list(mu = 10.0, lamda = rep(10, n)))

36 結果の要約 print(samples) Inference for Bugs model at "SevenScientists.txt", fit using jags, 2 chains, each with 1e+05 iterations (first 100 discarded), n.thin = 50 n.sims = 3996 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff mu sigma[1] sigma[2] sigma[3] sigma[4] sigma[5] sigma[6] sigma[7] deviance 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 = var(deviance)/2) pd = 14.1 and DIC = 37.3 DIC is an estimate of expected predictive error (lower deviance is better).

37 mixing の確認 traceplot traceplot(samples)

38 MCMC の 2 つの連鎖をそれぞれ取り出す mu <- samples$bugsoutput$sims.array[,,2] max(mu) # [1] min(mu) # [1] plot(mu[,1],ylim=c(7,14),ty='l',col=2) par(new=t) plot(mu[,2],ylim=c(7,14),ty='l',col=4)

40 summary(mu[,1]) Min. 1st Qu. Median Mean 3rd Qu. Max summary(mu[,2]) Min. 1st Qu. Median Mean 3rd Qu. Max density(mu[,1])$x[which(density(mu[,1])$y==max(density(mu[,1])$y))] [1] density(mu[,2])$x[which(density(mu[,2])$y==max(density(mu[,2])$y))] [1]

41 Nbreaks <- 80 x <- hist(mu[,1],nbreaks,plot=f) y <- hist(mu[,2],nbreaks,plot=f) plot(c(x$breaks, max(x$breaks)), c(0,x$density,0), type="s", lwd=2, lty=1, xlim=c(7,12), ylim=c(0,7), xlab="mu", ylab="posterior Density", col=2) par(new=t) plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="s", lwd=2, lty=1, xlim=c(7,12), ylim=c(0,7), xlab="mu", ylab="posterior Density", col=4) 2 つの MCMC で得られた μ の事後分布がほぼ同じ形に収束していることがわかる

42 sigma <- samples$bugsoutput$sims.array[,,3:9] max(sigma[,,7]) # [1] min(sigma[,,7]) # [1] plot(sigma[,1,7],ylim=c(0,790),ty='l',col=2) par(new=t) plot(sigma[,2,7],ylim=c(0,790),ty='l',col=4)

44 for (i in 1:7) { print(i) for (j in 1:2) { print(summary(sigma[,j,i])) print(density(sigma[,j,i])$x[which(density(sigma[,j,i])$y==max(density(sigma[,j,i])$y))]) } print("//") } [1] 1 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//" [1] 2 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//"

45 [1] 3 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//" [1] 4 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//" [1] 5 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//"

46 [1] 6 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//" [1] 7 Min. 1st Qu. Median Mean 3rd Qu. Max [1] Min. 1st Qu. Median Mean 3rd Qu. Max [1] [1] "//"

バイオインフォマティクス特論4

バイオインフォマティクス特論4 藤博幸 1-3-1. ピアソン相関係数 1-3-2. 致性のカッパ係数 1-3-3. 時系列データにおける変化検出ベイズ統計で実践モデリング 5.1 ピアソン係数第 5 章データ解析の例データは n ペアの独な観測値の対例 : 特定の薬剤の投与量と投与から t 時間後の注する遺伝の発現量 2 つの変数間の線形の関係性はピアソンの積率相関係数 r で表現される t 時間後の注する遺伝