SSブログ

負の二項分布のベイズ推定(2) [統計]

負の二項分布のベイズ推定の続き。今度はもうちょっと複雑なモデルを。

用意したデータ
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 

nice!(0)  コメント(0)  トラックバック(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0