MCMC / Drkcore

JAGSでマルコフ連鎖モンテカルロ法

Wed, 23 Sep 2009 17:03:22 +0919

連休は旅行いったりとか子供と遊んでたんだけどそれ以外はJAGSいじってた。

マルコフ連鎖モンテカルロ法 (統計ライブラリー)
豊田秀樹
朝倉書店 / ￥ 4,410 ()
在庫あり。

5章をWikiに書いていくことにした。

マルコフ連鎖モンテカルロ法をJAGSで

分散分析のあたりを出来るようにしておくとよいと思った。

5.1　ロジスティック回帰モデル　
5.2　メタ分析　
5.3　多項ロジットモデル　
5.4　対数線形モデル　
5.5　ポアソン回帰　
5.6　2値データに対する回帰分析　
5.7　トービット回帰モデル　
5.8　変曲点のある回帰分析　
5.9　生存時間分析(ワイブル回帰)　
5.10　生存時間分析(コックス回帰)　
5.11　時系列モデル　
5.12　分散分析　
5.13　分散成分分析　
5.14　分散分析(枝分かれ配置)　
5.15　一般化可能性理論　
5.16　反復測定データの分散分析　
5.17　階層線形モデル　
5.18　項目反応理論(2母数2値モデル)　
5.19　項目反応理論(段階反応モデル)　
5.20　項目反応理論(名義反応モデル)　
5.21　項目反応理論(部分採点・評定尺度モデル)
5.22　項目反応理論(連続反応モデル)
5.23　多次元IRT
5.24　項目反応理論(混合名義反応モデル)　
5.25　項目反応理論における特異項目機能(DIF) の分析　
5.26　正規混合モデル　
5.27　潜在クラス分析　
5.28　成長曲線モデル　
5.29　非線形成長曲線モデル　
5.30　因子分析　
5.31　多母集団分析　
5.32　非線形SEM

あとこの本だけだと難しすぎるのでPRMLの下巻が必須(多分)。

パターン認識と機械学習下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 8,190 ()
在庫あり。

JAGSでメタアナリシス

Sun, 20 Sep 2009 07:22:25 +0919

マルコフ連鎖モンテカルロ法の5-2をJAGSで

マルコフ連鎖モンテカルロ法 (統計ライブラリー)
豊田秀樹
朝倉書店 / ￥ 4,410 ()
在庫あり。

モデル

model
{
  for( i in 1 : Num ) {
    rc[i] ~ dbin(pc[i], nc[i]);
    rt[i] ~ dbin(pt[i], nt[i]);
    logit(pc[i]) <- mu[i];
    logit(pt[i]) <- mu[i] + delta[i];
    mu[i] ~ dnorm(0.0,1.0E-5);
    delta[i] ~ dnorm(d, tau);
  }
  d ~ dnorm(0.0,1.0E-6);
  tau ~ dgamma(0.001,0.001);
  delta.new ~ dnorm(d, tau);
  sigma <- 1 / sqrt(tau);
}

Rのコード

library(rjags)

data <- list(rt = c(49, 44, 102, 32, 85, 246, 1570),
     nt = c(615, 758, 832, 317, 810, 2267, 8587),
     rc = c(67, 64, 126, 38, 52, 219, 1720),
     nc = c(624, 771, 850, 309, 406, 2257, 8600),
    Num = 7)

inits <- list(d = 0, delta.new = 0, tau=1,
     mu = c(0, 0, 0, 0, 0, 0, 0),
  delta = c(0, 0, 0, 0, 0, 0, 0))

m <- jags.model(
  file = "model.txt",
  data = data,  
  inits = list(inits),
  n.chains = 1
)


update(m,5000)

x <- coda.samples(
    m,
    c("d","sigma"),
    thin = 1, n.iter = 15000
)

summary(x)

plot(x)

結果

> summary(x)

Iterations = 6001:21000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 15000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean      SD  Naive SE Time-series SE
d     -0.1400 0.08028 0.0006555       0.001885
sigma  0.1285 0.08188 0.0006686       0.002196

2. Quantiles for each variable:

          2.5%      25%     50%     75%    97.5%
d     -0.31570 -0.18244 -0.1338 -0.0908 0.003366
sigma  0.02855  0.06916  0.1110  0.1672 0.331475

JAGSでロジスティック回帰モデル

Sat, 19 Sep 2009 17:06:08 +0919

macでwinbugsを頑張って動かすのもアレなので、潔くJAGSで。

マルコフ連鎖モンテカルロ法の5-1

マルコフ連鎖モンテカルロ法 (統計ライブラリー)
豊田秀樹
朝倉書店 / ￥ 4,410 ()
在庫あり。

生態学のデータ解析 - R で JAGSによるとjagsの場合にはBUGS言語の文（式）の終わりにはセミコロンが必要。

model
{
  for(i in 1:N){
    r[i] ~ dbin(p[i],n[i]);
    e[i] ~ dnorm(0.0,tau);
    logit(p[i]) <- beta0+beta1*x1[i]+beta21*x21[i]+beta22*x22[i]
    +beta23*x23[i]+beta24*x24[i]+beta25*x25[i]
    +beta26*x26[i]+beta27*x27[i]+e[i];
  }
  beta0 ~ dnorm(0.0,1.0E-6);
  beta1 ~ dnorm(0.0,1.0E-6);
  beta21 ~ dnorm(0.0,1.0E-6);
  beta22 ~ dnorm(0.0,1.0E-6);
  beta23 ~ dnorm(0.0,1.0E-6);
  beta24 <- 0;
  beta25 ~ dnorm(0.0,1.0E-6);
  beta26 ~ dnorm(0.0,1.0E-6);
  beta27 ~ dnorm(0.0,1.0E-6);
  tau ~ dgamma(0.001,0.001);
  sigma<-1/sqrt(tau);
}

Rのコード

library(rjags)

data <- list(
r = c(1,13,44,155,92,100,18,0,0,2,1,0,2,0),
n = c(13,70,115,301,153,141,26,6,16,25,32,8,9,4),
x1 = c(0,0,0,0,0,0,0,1,1,1,1,1,1,1),
x21 = c(1,0,0,0,0,0,0,1,0,0,0,0,0,0),
x22 = c(0,1,0,0,0,0,0,0,1,0,0,0,0,0),
x23 = c(0,0,1,0,0,0,0,0,0,1,0,0,0,0),
x24 = c(0,0,0,1,0,0,0,0,0,0,1,0,0,0),
x25 = c(0,0,0,0,1,0,0,0,0,0,0,1,0,0),
x26 = c(0,0,0,0,0,1,0,0,0,0,0,0,1,0),
x27 = c(0,0,0,0,0,0,1,0,0,0,0,0,0,1),
N = 14)

inits <- list(beta0 = 0, beta1 = 0, beta21 = 0,  beta22 = 0, beta23 = 0, beta24 = 0, beta25 = 0, beta26 = 0,  beta27 = 0, tau = 1)

m <- jags.model(
  file = "model.txt",
  data = data,  
  inits = list(inits),
  n.chains = 1
)

update(m,10000)

x <- coda.samples(
    m,
    c("beta0","beta1","beta21","beta22","beta23","beta24","beta25","beta26","beta27", "sigma"),
    thin = 1, n.iter = 40000
)

summary(x)

plot(x)

サマリーの表示

> summary(x)

Iterations = 11001:51000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 40000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean     SD Naive SE Time-series SE
beta0  -0.02576 0.4188 0.002094        0.03425
beta1  -3.11738 0.6059 0.003030        0.01815
beta21 -3.03168 1.4556 0.007278        0.03269
beta22 -1.52877 0.6757 0.003379        0.02993
beta23 -0.28549 0.6355 0.003177        0.04374
beta24  0.00000 0.0000 0.000000        0.00000
beta25  0.35897 0.5926 0.002963        0.02957
beta26  1.05996 0.6389 0.003194        0.04301
beta27  0.70287 0.7677 0.003838        0.03026
sigma   0.34110 0.3980 0.001990        0.03152

2. Quantiles for each variable:

           2.5%      25%      50%      75%   97.5%
beta0  -1.11529 -0.15412  0.02137  0.16865  0.7406
beta1  -4.48731 -3.46881 -3.06188 -2.70570 -2.0882
beta21 -6.39669 -3.83976 -2.86884 -2.05688 -0.6382
beta22 -2.90198 -1.87085 -1.54216 -1.21238 -0.0738
beta23 -1.27301 -0.63488 -0.39022 -0.08097  1.4261
beta24  0.00000  0.00000  0.00000  0.00000  0.0000
beta25 -0.87361  0.10421  0.34604  0.59078  1.7201
beta26  0.08425  0.72545  0.96155  1.25666  2.7809
beta27 -0.80720  0.33750  0.71085  1.10396  2.1284
sigma   0.02816  0.08266  0.19387  0.45606  1.3892

5章は事例がたくさん載っているので、この勢いで一通りやっておこう。

5.1　ロジスティック回帰モデル　
5.2　メタ分析　
5.3　多項ロジットモデル　
5.4　対数線形モデル　
5.5　ポアソン回帰　
5.6　2値データに対する回帰分析　
5.7　トービット回帰モデル　
5.8　変曲点のある回帰分析　
5.9　生存時間分析(ワイブル回帰)　
5.10　生存時間分析(コックス回帰)　
5.11　時系列モデル　
5.12　分散分析　
5.13　分散成分分析　
5.14　分散分析(枝分かれ配置)　
5.15　一般化可能性理論　
5.16　反復測定データの分散分析　
5.17　階層線形モデル　
5.18　項目反応理論(2母数2値モデル)　
5.19　項目反応理論(段階反応モデル)　
5.20　項目反応理論(名義反応モデル)　
5.21　項目反応理論(部分採点・評定尺度モデル)
5.22　項目反応理論(連続反応モデル)
5.23　多次元IRT
5.24　項目反応理論(混合名義反応モデル)　
5.25　項目反応理論における特異項目機能(DIF) の分析　
5.26　正規混合モデル　
5.27　潜在クラス分析　
5.28　成長曲線モデル　
5.29　非線形成長曲線モデル　
5.30　因子分析　
5.31　多母集団分析　
5.32　非線形SEM

Machine Learning: An Algorithmic Perspective

Tue, 08 Sep 2009 20:11:26 +0919

一ヶ月近くかかってやっと到着

Machine Learning: An Algorithmic Perspective (Chapman & Hall/Crc Machine Learning & Patrtern Recognition)
Stephen Marsland
Chapman & Hall / ￥ 6,707 ()
通常1～3か月以内に発送

機械学習をpythonで的な本だと思っていたが、パラパラめくってみると、プログラムを丁寧に解説している本ではなくて、理論の説明を助けるためにコードを載せている感じ。

Introduction.
Linear Discriminants.
The Multi-Layer Perceptron.
Radial Basis Functions and Splines.
Support Vector Machines.
Learning with Trees.
Decision by Committee: Ensemble Learning.
Probability and Learning. Unsupervised Learning.
Dimensionality Reduction.
Optimisation and Search.
Evolutionary Learning.
Reinforcement Learning.
Markov Chain Monte Carlo (MCMC) Methods.
Graphical Models.

そのうち読む。

SIR(sampling-importance-resampling)

Tue, 18 Aug 2009 22:13:35 +0919

水産資源学におけるベイズ統計の応用ワークショップのスライドの22,23枚目

g(x)からM個のサンプルx1*, …, xM*を発生
w*= f(x*)/g(x*)を計算
x1*, …, xM*からw1*, …, wM*に比例する確率でサンプルをm個とってくる
得られたm個のサンプルはf(x)からの事後サンプル

となっているのだが、23枚目のRのコードは

n2 <-sample(length(x1), 1000, replace=T, prob=dlnorm(x1,1.1,0.6))

とwを求めずにいきなりサンプリングしているので？となったが、一様分布だから重みを計算する必要ないのね。

わざわざ書くならこうか。

x1 <-runif(100000,0,60)
w <- dlnorm(x1,1.1,0.6) / (dunif(x1,0,60))
n2 <-sample(length(x1), 1000, replace=T, prob=w)
x2 <-x1[n2]
plot(density(x2))

パターン認識と機械学習下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 8,190 ()
在庫あり。

マルコフ連鎖モンテカルロ法

Mon, 17 Aug 2009 21:29:46 +0919

夏休み前にBUGSの具体的な事例が載っていたので買って積んでおいた本を読みはじめた。さっきまでWinBUGSをmacから使えるように頑張っていたけど、JAGSもpymcもあるし、無理して使わんでもまぁええかと諦めた。

マルコフ連鎖モンテカルロ法 (統計ライブラリー)
豊田秀樹
朝倉書店 / ￥ 4,410 ()
在庫あり。

統計科学のフロンティアシリーズのやつも薦められたので読んでみたい。

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア)
伊庭幸人,種村正美
岩波書店 / ￥ 4,620 ()
在庫あり。

伊庭先生と言えば「ベイズ統計と統計物理」にはお世話になった。

岩波講座物理の世界物理と情報（3）ベイズ統計と統計物理
伊庭幸人
岩波書店 / ￥ 1,470 ()
在庫あり。

で、この本買ったのが、GIWとなんかの共同開催の時で、基調講演が赤池先生だったのは覚えている。

pythonでMCMC

Wed, 15 Jul 2009 22:01:52 +0919

最近、Journal of Pharmaceutical SciencesみたいなDMPKな雑誌を読むようになってきた。その中でScipyとpymcっていうpythonでMCMCができるモジュールを使って解析してた論文を見たので試したくなった。

sudo easy_install-2.6 pymc

でさくっと入って、テストしてみたらなんか色々こけてるようなのだけど、とりあえずサンプル写経

import pymc
import numpy as np
n = 5*np.ones(4,dtype=int)
x = np.array([-.86,-.3,-.05,.73])
alpha = pymc.Normal('alpha',mu=0,tau=.01)
beta = pymc.Normal('beta',mu=0,tau=.01)

@pymc.deterministic
def theta(a=alpha, b=beta):
    """theta = logit^{-1}(a+b)"""
    return pymc.invlogit(a+b*x)

d = pymc.Binomial('d', n=n, p=theta,value=np.array([0.,1.,3.,5.]),observed=True)

これをmymodel.pyって名前で保存しておいて

対話環境から

>>> import mymodel
>>> import pymc
>>> S = pymc.MCMC(mymodel, db='pickle')
>>> S.sample(iter=10000, burn=5000, thin=2)
>>> pymc.Matplot.plot(S)
>>> pymc.Matplot.savefig("test.png")

うーんMCMCわからん。精進せねば。

パターン認識と機械学習上 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 6,825 ()
在庫あり。

パターン認識と機械学習下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 8,190 ()
在庫あり。

マルコフ連鎖モンテカルロ

Sun, 19 Nov 2006 22:09:17 +0919

SVMの流れで、MCMCとかベイズあたりを読んでる。

で、昔買ったけど最初の方の章だけ読んで放り投げてあったベイズ統計と統計物理を読んでたら、MCMCの話とかちゃんと書いてあって、面白かった。

岩波講座物理の世界物理と情報〈3〉ベイズ統計と統計物理
伊庭幸人
岩波書店 / 1470円 ( 2003-08-27 )

で、尤度とか、事前分布、事後分布がまだイメージできないナァと、ちょっとググッたら、Bayesianってどういう考え方なんだろうというページがみつかった。

事後確率 = (尤度 * 事前確率)/エビデンス

エビデンスの部分の積分が困難なので、事後分布のシミュレーションをMCMCでやると。

ベイズ統計学のMCMCとの出会い

こんな感じの理解でいいのか?