トップ
新規
一覧
単語検索
最終更新
ヘルプ
ログイン
アールメカブ
Bayes_Poisson
をテンプレートにして作成
開始行:
[[Rの備忘録]]
授業で使います．
#contents
* ポアソンの共役事前分布はガンマ分布とする [#l4b46a58]
Gelman et al.: Bayesian Data Analysis, p.52, 2nd Edition．
ポアソン分布とガンマ分布の関係については繁桝算男『ベイズ...
-ポアソン分布
#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{...
#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 Statis...
#mimetex(p(\bf{x}| \theta) = g(\bf{x}) \, h(\theta) \, e^...
上の式では&mimetex(\frac{1}{y!});が消えているので注意
-自然共役事前分布
#mimetex(p(\theta) \propto (e^{-\theta})^n e^{v log \thet...
- 事後分布
#mimetex(\theta|\bf{y} \sim Gamma(\alpha + n \bar{y}, \be...
- 医学関係ではポアソン分布を次のように書いている．
#mimetex(y_i \sim Poisson(x_i \theta))
&mimetex(y_i); と &mimetex(x_i); が観測数で，&mimetex(\th...
* R での実行例 [#pd06f7ab]
- Albert Bayesian: Computation with R より
-- 米国のある病院での心臓手術が失敗するケース数&mimetex(y...
-- 米国のある病院での心臓手術が失敗するケース数&mimetex(y...
#mimetex(y \sim Poisson(e \lambda))
-- &mimetex(\lambda);の最尤推定量は &mimetex( \bar{\lambd...
-- 上で述べたように事前共役分布はガンマ分布だが，二つのパ...
#mimetex(\theta^{16-1} \, e^{-15174\, \theta} = Gamma(16,...
-- すると事後分布は &mimetex(16 + y, 15174 + e); のガンマ...
-- さて，ここで新たに手術例が66件で，うち失敗が1件の病院...
my.alpha <- 16
my.beta <- 15174
lam <- my.alpha/my.beta
lambdaA <- rgamma(1000, shape = my.alpha + 1, rate = my....
lambdaB <- rgamma(1000, shape = my.alpha + 1767, rate = ...
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|\la...
#mimetex(f(y) = \frac{f(y|\lambda) \, g(\lambda)}{g(y|\la...
#mimetex(p(y) = \frac{Poisson(y|\theta) Gamma(\theta|\alp...
#mimetex(=\frac{\Gamma (\alpha + y) \beta^\alpha}{\Gamma(...
#mimetex(y \sim Neg-bin(\alpha, \beta))
これはポアソン分布の&mimetex(\theta); がガンマ分布に従う...
* ポアソン回帰 [#z331fdcf]
豊田秀樹『マルコフ連鎖モンテカルロ法』は良い本なんだけど...
以下，『マルコフ連鎖モンテカルロ法』p.107掲載データ（生化...
データ&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 variab...
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/st...
終了行:
[[Rの備忘録]]
授業で使います．
#contents
* ポアソンの共役事前分布はガンマ分布とする [#l4b46a58]
Gelman et al.: Bayesian Data Analysis, p.52, 2nd Edition．
ポアソン分布とガンマ分布の関係については繁桝算男『ベイズ...
-ポアソン分布
#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{...
#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 Statis...
#mimetex(p(\bf{x}| \theta) = g(\bf{x}) \, h(\theta) \, e^...
上の式では&mimetex(\frac{1}{y!});が消えているので注意
-自然共役事前分布
#mimetex(p(\theta) \propto (e^{-\theta})^n e^{v log \thet...
- 事後分布
#mimetex(\theta|\bf{y} \sim Gamma(\alpha + n \bar{y}, \be...
- 医学関係ではポアソン分布を次のように書いている．
#mimetex(y_i \sim Poisson(x_i \theta))
&mimetex(y_i); と &mimetex(x_i); が観測数で，&mimetex(\th...
* R での実行例 [#pd06f7ab]
- Albert Bayesian: Computation with R より
-- 米国のある病院での心臓手術が失敗するケース数&mimetex(y...
-- 米国のある病院での心臓手術が失敗するケース数&mimetex(y...
#mimetex(y \sim Poisson(e \lambda))
-- &mimetex(\lambda);の最尤推定量は &mimetex( \bar{\lambd...
-- 上で述べたように事前共役分布はガンマ分布だが，二つのパ...
#mimetex(\theta^{16-1} \, e^{-15174\, \theta} = Gamma(16,...
-- すると事後分布は &mimetex(16 + y, 15174 + e); のガンマ...
-- さて，ここで新たに手術例が66件で，うち失敗が1件の病院...
my.alpha <- 16
my.beta <- 15174
lam <- my.alpha/my.beta
lambdaA <- rgamma(1000, shape = my.alpha + 1, rate = my....
lambdaB <- rgamma(1000, shape = my.alpha + 1767, rate = ...
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|\la...
#mimetex(f(y) = \frac{f(y|\lambda) \, g(\lambda)}{g(y|\la...
#mimetex(p(y) = \frac{Poisson(y|\theta) Gamma(\theta|\alp...
#mimetex(=\frac{\Gamma (\alpha + y) \beta^\alpha}{\Gamma(...
#mimetex(y \sim Neg-bin(\alpha, \beta))
これはポアソン分布の&mimetex(\theta); がガンマ分布に従う...
* ポアソン回帰 [#z331fdcf]
豊田秀樹『マルコフ連鎖モンテカルロ法』は良い本なんだけど...
以下，『マルコフ連鎖モンテカルロ法』p.107掲載データ（生化...
データ&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 variab...
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/st...
ページ名:
