負の二項分布のベイズ推定(2) [統計]
負の二項分布のベイズ推定の続き。今度はもうちょっと複雑なモデルを。
用意したデータ
BUGSコード
結果
glm.nbの結果と比較してみる。
用意したデータ
set.seed(13) n <- 100 size <- 7 x <- runif(n, 0, 4) mu <- exp(0.8 + x * 0.4) num <- rnbinom(n, size=size, mu=mu)
BUGSコード
var N, size, mu[N], p[N], num[N], x[N] model { for (i in 1:N) { num[i] ~ dnegbin(p[i], size); p[i] <- size / (mu[i] + size); log(mu[i]) <- b + b.x * (x[i] - mx); } mx <- mean(x[]); beta <- b - b.x * mx; # priors b ~ dnorm(0.0, 1.0E-6); b.x ~ dnorm(0.0, 1.0E-6); size ~ dgamma(1.0E-3, 1.0E-3); }
結果
Iterations = 0:4995 Thinning interval = 5 Number of chains = 3 Sample size per chain = 1000 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE size 7.3777 3.37353 0.061592 0.069462 beta 0.6688 0.13873 0.002533 0.002526 b.x 0.4656 0.05378 0.000982 0.000962 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% size 3.4809 5.149 6.5979 8.6796 16.1034 beta 0.3840 0.576 0.6691 0.7643 0.9406 b.x 0.3583 0.429 0.4649 0.5032 0.5737
glm.nbの結果と比較してみる。
> summary(glm.nb(num ~ x)) Call: glm.nb(formula = num ~ x, init.theta = 6.59845061324799, link = log) Deviance Residuals: Min 1Q Median 3Q Max -2.8104 -0.8820 -0.1725 0.5181 2.6780 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.67062 0.13768 4.871 1.11e-06 *** x 0.46487 0.05393 8.619 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for Negative Binomial(6.5985) family taken to be 1) Null deviance: 200.04 on 99 degrees of freedom Residual deviance: 121.05 on 98 degrees of freedom AIC: 499.59 Number of Fisher Scoring iterations: 1 Theta: 6.60 Std. Err.: 2.39 2 x log-likelihood: -493.592
タグ:負の二項分布
2008-05-01 18:37
nice!(0)
コメント(0)
トラックバック(0)
コメント 0