人生で初めて著書を贈呈していただきました。
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2015/10/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (10件) を見る
レビュー書く。
COI：著者の一人より贈呈されました。
謹呈という帯が入っていました。ありがとうございました。というわけで優し目のコメントになっております。
良かった点
ベイズ推論とMCMCのフリーソフト、という特集を掲げているだけあって、MCMCの具体的な例をあげてRコードも適宜載せつつ解説されていた。
特に例については実際にありえそうなテーマを出して、
というわけでMCMCやR, Python などその他ソフトウェアに馴染みのある人ならば理解しやすい。
インストール、みたいな初級者ページを完全すっ飛ばしているので無駄は少ない(その分初級者お断り感ある)。
increment_log_prob の使い方に馴染みがなかったため、ようやく使い方がわかった。
執筆陣はわりと豪華であり、しかも最近よく目にする()若手()の人も執筆陣として多数名を連ねている様子。
悪かった点
MCMC は実践的なやり方についてなので、理論を大事に派の人はたぶん物足りない。
初級者は完全に置いてけぼりになる。
サポートページは、まあ分かってたけど分かっている人向けな感じ。
脳とディープラーニングについては本当に概略だけというかそれだけ。
たぶんこの執筆陣のブログなどをよく見ている人にとってはほとんど既出の内容じゃないのか?
途中でしれっとサポートページに触れられてここからスクリプト拾える仕様だということに気づく。
雑誌的な機動力、というのでお手軽感はあるが、
総評
MCMC をやりはじめなぐらいで、行き詰まっているレベルの人は必須。
上級者は小一時間あれば読める。
初級者は日本語でおk. な感じになる。
基本的にJAGS で書かれているが、Rstan しか使ったことのない人としてはRstan で書きなおしてみる。中の人の例。
給食と身長の話題
data { int<lower=0> N_pref; int<lower=0> N_age; real<lower=0> Y[N_age, N_pref]; int<lower=0, upper=1> type[N_age, N_pref]; int<lower=0, upper=1> Age[N_age, N_pref]; real<lower=0> Tau_se[N_age, N_pref]; } parameters { real beta[3]; real r[N_age, N_pref]; } transformed parameters { real<lower=0> tau[N_age]; real<lower=0> SD[N_age]; real mu[N_age, N_pref]; for(i in 1:N_age){ tau[i] <- 1 / (SD[i] * SD[i]); for(j in 1:N_pref){ mu[i, j] <- beta[1] + r[1, j] + (beta[2] + beta[3] * type[i, j] + r[2, j]) * Age[i, j]; } } } model { for(b in 1:3) beta[b] ~ uniform(-1, 1); for (i in 1:N_age){ increment_log_prob(normal_log(SD, 0, 2)); SD ~ uniform(0, 2); for (j in 1:N_pref){ r[i, j] ~ normal(0, 1 / (SD[i] * SD[i])); Y[i, j] ~ normal(mu[i, j], Tau_se[i, j]); } } }
このモデルはtau がうまくいきませんって言われてうまくいかなかった。
サイトからスクリプトをパクってJAGSでやった。
dat <- read.delim("height.txt", stringsAsFactors=FALSE) summary(glm(height ~ age + age:type, data=dat))
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 152.9050 0.3978 384.341 < 2e-16 *** age 5.6490 0.6891 8.198 2.61e-07 *** age:type -1.7220 0.7957 -2.164 0.045 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 1.582744)
線形回帰では推定結果は一致する。
library(rstan) stanmodel <- stan_model("height_model02.stan") standata <- list(N_pref=length(unique(dat$site)), N_age=length(unique(dat$age)), Y=matrix(dat$height, nr=2, byrow=TRUE), type=matrix(dat$type, nr=2, byrow=TRUE), Age=matrix(dat$age, nr=2, byrow=TRUE), Tau_se=matrix(dat$size/(dat$sd^2), nr=2, byrow=TRUE) ) fit <- sampling(stanmodel, data=standata, chains=3, warmup=500, iter=1500, seed=1234) ex <- extract(fit) cols <- rainbow(N.pref) xl <- c(-1, 1)*3.5 yl <- c(0, 1) par(mfrow=c(2, 1), mar=c(2, 4, 2, 2), las=1) for(j in 1:2){ idx <- grep(paste("r\\[", j, ",", sep=""), colnames(a)) den <- apply(a[,idx], 2, density) plot(0, type="n", xlim=xl, ylim=yl, xlab="", ylab="") for(i in seq(den)) lines(den[[i]]$x, den[[i]]$y, col=cols[i], lwd=2) legend("topleft", legend=LETTERS[seq(N.pref)], col=cols, lty=1, lwd=2) }
1次元のCARモデル
increment_log_prob の使い方がわかったようなわからんかったような…
data { int<lower=0> N; real<lower=0> y[N]; } parameters { real alpha[N]; real sigma[2]; real tau[2]; } model { for(i in 1:2){ alpha[i] ~ normal(0, 10); sigma[i] ~ uniform(0, 100); } for(i in 1:N) #y[i] ~ normal(alpha[i], 1/sigma[1]^2); # 正規分布 increment_log_prob(normal_log(y[i], alpha[i], 1/sigma[1]^2)); # 対数正規分布だけどなんか違う for(i in 3:N) alpha[i] ~ normal(2*alpha[i-1] - alpha[i-2], 1/sigma[2]^2); }
X <- 1961:1990 Y <- c(4.71, 7.70, 7.97, 8.35, 5.70, 7.33, 3.10, 4.98, 3.75, 3.35, 1.84, 3.28, 2.77, 2.72, 2.54, 3.23, 2.45, 1.90, 2.56, 2.12, 1.78, 3.18, 2.64, 1.86, 1.69, 0.81, 1.02, 1.40, 1.31, 1.57) stanmodel <- stan_model("trend.stan") standata <- list(N=length(X), y=Y) fit <- sampling(stanmodel, data=standata, chains=3, warmup=500, iter=1500, seed=1234) ex <- extract(fit) ci <- apply(ex$alpha, 2, quantile, c(0.025, 0.975)) plot(X, Y, xlab="year", ylab="width", ylim=c(0, 10), type="n") polygon(c(X, rev(X)), c(ci[1,], rev(ci[2,])), col=grey(0.8)) points(X, Y, type="o", pch=16)
2次元のCARモデル
面倒だったのでパッケージ使った。
library(CARBayes) library(MASS) dat <- read.csv("Qglauca.csv") ## 隣接行列 n.x <- length(levels(as.factor(dat$X))) n.y <- length(levels(as.factor(dat$Y))) n.sites <- n.x * n.y W <- matrix(0, nrow = n.sites, ncol = n.sites) for (x in 0:(n.x - 1)) { for (y in 1:n.y) { if (x > 0) W[x * n.y + y, (x - 1) * n.y + y] <- 1 if (x < n.x - 1) W[x * n.y + y, (x + 1) * n.y + y] <- 1 if (y > 1) W[x * n.y + y, x * n.y + y - 1] <- 1 if (y < n.y) W[x * n.y + y, x * n.y + y + 1] <- 1 } } fit.iar <- S.CARiar(N ~ 1, family = "poisson", data = dat, W = W, burnin = 2000, n.sample = 32000, thin = 10) geweke.diag(fit.iar$samples$beta) plot(fit.iar$samples$beta) A <- list(matrix(dat$N, max(dat$X), max(dat$Y), byrow=TRUE), matrix(fit.iar$fitted.values, max(dat$X), max(dat$Y), byrow=TRUE)) par(mfrow=c(1, 2)) for(i in seq(A)) image(A[[i]], col=grey(seq(1, 0, length=12)))
