Site admin: motohiro PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 7.4.33
HTML convert time: 0.064 sec.
[[Rの備忘録]] 授業で使います． #contents * ポアソンの共役事前分布はガンマ分布である [#l4b46a58] (Gelman et al.: Bayesian Data Analysis, p.52, 2nd Edition)． * ポアソンの共役事前分布はガンマ分布とする [#l4b46a58] Gelman et al.: Bayesian Data Analysis, p.52, 2nd Edition． ポアソン分布とガンマ分布の関係については繁桝算男『ベイズ統計入門』64ページを参照． -ポアソン分布 #mimetex(p(y|\theta) = \frac{\theta^y \, e^{-\theta}}{y!}) - 尤度 #mimetex(\bf{y} = y_i, \dots , y_n) #mimetex(p(\bf{y}|\theta) = l(\theta|\bf{y}) = \Pi \frac{1}{y_i!} \theta^{y_i} \, e^{-\theta}) #mimetex(p(\bf{y}|\theta) = l(\theta|\bf{y}) = \Pi \frac{1}{y_i!} \theta^{y_i} \, e^{-\theta}) #mimetex( \propto \theta^{t(y)} \, e^{-n \theta}) 指数属で表現すると #mimetex( \propto e^{-n\theta}\, e^{t(y) log \theta} ) 文系学生相手にするので，指数法則を指摘する． #mimetex(\theta^y = e^{log(\theta^y)} = e^{y log \theta}) 指数属とは関数が次の形で表されること(Lee: Bayesian Statistics, pp.60)． #mimetex(p(\bf{x}| \theta) = g(\bf{x}) \, h(\theta) \, e^{t(\bf{x}) \, \phi(\theta)}) 上の式では&mimetex(\frac{1}{y!});が消えているので注意 -自然共役事前分布 #mimetex(p(\theta) \propto (e^{-\theta})^n e^{v log \theta} ) - 事後分布 #mimetex(\theta|\bf{y} \sim Gamma(\alpha + n \bar{y}, \beta + n)) - 医学関係ではポアソン分布を次のように書いている． #mimetex(y_i \sim Poisson(x_i \theta)) &mimetex(y_i); と &mimetex(x_i); が観測数で，&mimetex(\theta); が未知のパラメータだが，&mimetex(x); を exposure，&mimetex(\theta); をrate と表現しているので，分野違いの人間は，何度みても慣れないというか，馴染めない． &mimetex(y_i); と &mimetex(x_i); が観測数で，&mimetex(\theta); が未知のパラメータだが，&mimetex(x); を exposure，&mimetex(\theta); をrate と表現しており，分野違いの人間はいつまでたっても馴染みにくい． * R での実行例 [#pd06f7ab] - Albert Bayesian: Computation with R より -- 米国のある病院での心臓手術が失敗するケース数&mimetex(y); を考える．ここで失敗とは手術後30以内に死亡することである．これはポアソン分布に従うとする．手術(exposure)の総数を&mimetex(e);，exposureあたりの死亡率を&mimetex(\lambda); とすると， -- 米国のある病院での心臓手術が失敗するケース数&mimetex(y); を考える．ここで失敗とは手術後30以内に死亡することである．これはポアソン分布に従うとする．手術(exposure)の総数を&mimetex(e);，exposureあたりの死亡率を&mimetex(\lambda); とすると， #mimetex(y \sim Poisson(e \lambda)) -- &mimetex(\lambda);の最尤推定量は &mimetex( \bar{\lambda} = y/e);である． -- 上で述べたように事前共役分布はガンマ分布だが，二つのパラメータを設定する必要がある．ここでは，いま対象としている病院と同じ手術レベルの他の病院のデータを参考にして，&mimetex(\alpha = 16); と&mimetex(\beta = 15174); とする．すなわち事前分布は #mimetex(\theta^{16-1} \, e^{-15174\, \theta} = Gamma(16,15174 )) -- すると事後分布は &mimetex(16 + y, 15174 + e); のガンマ分布となる．（事前）予測密度関数&mimetex(f(x)); は，&mimetex(f(y|\lambda)); をサンプルのポアソン分布，&mimetex(g(\lambda)); を事前分布，&mimetex(g(y|\lambda)); を事後分布とすると， -- すると事後分布は &mimetex(16 + y, 15174 + e); のガンマ分布となる． -- さて，ここで新たに手術例が66件で，うち失敗が1件の病院と，手術例が1767件で，うち失敗が4件の病院があるとする． my.alpha <- 16 my.beta <- 15174 lam <- my.alpha/my.beta lambdaA <- rgamma(1000, shape = my.alpha + 1, rate = my.beta + 66) lambdaB <- rgamma(1000, shape = my.alpha + 1767, rate = my.beta + 4) lambda <- seq(0, max(c(lambdaA, lambdaB)), length = 500) par(mfrow = c(2,1))# , mar = rep(1, 4)) hist(lambdaA, freq = FALSE, main = "", ylim =c(0, 1600)) lines(lambda, dgamma(lambda, shape = my.alpha, my.beta), col = "blue", lwd = 3) lines(lambda, dgamma(lambda, shape = my.alpha+ 1, my.beta + 66), col = "red", lwd = 3) legend(0.0015, 1500, legend= c("prior", "posterior"), col = c("blue","red"), lwd = 3) hist(lambdaB, freq = FALSE, main = "", ylim =c(0, 1600)) lines(lambda, dgamma(lambda, shape = my.alpha, my.beta), col = "blue", lwd = 3) lines(lambda, dgamma(lambda, shape = my.alpha+ 4, my.beta + 1767), col = "red", lwd = 3) legend(0.0015, 1500, legend= c("prior", "posterior"), col = c("blue","red"), lwd = 3) #ref(bayes.poisson.png,nowrap,nolink) * 予測分布と負の二項分布 [#e0786ce2] -（事前）予測密度関数&mimetex(f(x)); は，&mimetex(f(y|\lambda)); をサンプルのポアソン分布，&mimetex(g(\lambda)); を事前分布，&mimetex(g(y|\lambda)); を事後分布とすると， #mimetex(f(y) = \frac{f(y|\lambda) \, g(\lambda)}{g(y|\lambda)}) -- さて，ここで新たに手術例が66件で，うち失敗が1件の病院について，その失敗率を知りたいとする．単純には 1/66 であるが，失敗例が1件だけなので，この値は信頼に足るだろうか． #mimetex(p(y) = \frac{Poisson(y|\theta) Gamma(\theta|\alpha, \beta)}{Gamma(\theta|\alpha+y, 1 + \beta)}) #mimetex(=\frac{\Gamma (\alpha + y) \beta^\alpha}{\Gamma(\alpha)y! \, (1+\beta)^{\alpha+y}}) #mimetex(y \sim Neg-bin(\alpha, \beta)) これはポアソン分布の&mimetex(\theta); がガンマ分布に従う混合分布である． * ポアソン回帰 [#z331fdcf] 豊田秀樹『マルコフ連鎖モンテカルロ法』は良い本なんだけど，掲載サンプルではBRugsパッケージの利用を前提としている．ところがBRugsはCRANのパッケージから外されている．ポアソン回帰ならば，''MCMCpack''でも間に合う． 以下，『マルコフ連鎖モンテカルロ法』p.107掲載データ（生化学専攻大学院生の論文執筆数を性別，既婚・独身，子供の数で説明する）を glm と MCMCpoisson で実行する． データ&ref(MCMCpoisson.R);を読み込んで library(MCMCpack) source("MCMCpoisson.R") posterior <- MCMCpoisson(y ~ x1 + x2 + x3) plot(posterior) summary(posterior) Iterations = 1001:11000 Thinning interval = 1 Number of chains = 1 Sample size per chain = 10000 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE (Intercept) 0.6389 0.05432 0.0005432 0.002321 x1 -0.2852 0.05355 0.0005355 0.002204 x2 0.1300 0.06095 0.0006095 0.002382 x3 -0.1623 0.03891 0.0003891 0.001549 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% (Intercept) 0.5323 0.60324 0.6397 0.6746 0.74692 x1 -0.3908 -0.32266 -0.2837 -0.2498 -0.18105 x2 0.0082 0.09005 0.1311 0.1703 0.25447 x3 -0.2396 -0.18926 -0.1622 -0.1360 -0.08486 一般線形化モデル pois.glm <- glm(y ~ x1 + x2 + x3,family=poisson()) summary(pois.glm) Call: glm(formula = y ~ x1 + x2 + x3, family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -2.0776 -1.6856 -0.3729 0.4936 6.9981 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.63656 0.05457 11.665 < 2e-16 *** x1 -0.28549 0.05433 -5.255 1.48e-07 *** x2 0.13271 0.06092 2.179 0.0294 * x3 -0.16118 0.03934 -4.097 4.19e-05 *** (Dispersion parameter for poisson family taken to be 1) Null deviance: 1817.4 on 914 degrees of freedom Residual deviance: 1776.7 on 911 degrees of freedom AIC: 3452.5 #ref(MCMCpoisson.png,nowrap,nolink) Number of Fisher Scoring iterations: 5 * 参考文献 [#y1bca176] Gelman et al.: Bayesian Data Analysis, 2nd Edition,Chapman Lee: Bayesian Statistics, Hodder Arnold Albert Bayesian: Computation with R, Springer 繁桝算男『ベイズ統計入門』東京大学 豊田秀樹『マルコフ連鎖モンテカルロ法』朝倉書店 こういう[[ページ:http://www001.upp.so-net.ne.jp/ito-hi/stat/R2.html]]も参考に．