########################################################
# Section 7.5 Modeling a Prior Belief of
Exchangeability
########################################################
model
{
#
likelihood
for(i
in 1:n)
{
y[i]
~ dpois(e.lambda[i])
e.lambda[i]
<- e[i]*lambda[i]
lambda[i] ~ dgamma(alpha, beta)
loge[i]
<- log(e[i])
}
beta <-
alpha*mu.inv
# prior for mu
mu.inv ~
dgamma(0.001, 0.001)
mu <-
1/mu.inv
# prior for
alpha
alpha <- x-z0
z0 <- 0.53 # median of alpha
x ~ dpar(1,
z0) # f(x) = z0/x^2, x > z0
(x=alpha+z0)
lalpha <-
log(alpha)
lmu <- log(mu)
# Compute
shrinkage size toward the pooled estimate
for(i in 1:n)
{
B[i]
<- alpha/(alpha+e[i]*mu)
}
# Posterior
predictive analysis
# Is one of
E(pright[i]) and E(pleft[i]) less than 0.05?
# If yes, the
model may not be valid to describe i-th observation.
for(i in 1:n)
{
ypred[i]
~ dpois(e.lambda[i])
pright[i]
<- step(ypred[i] - y[i])
pleft[i] <- step(y[i] - ypred[i])
}
}
# Initial values. Initialize other values using
"gen inits".
list( mu.inv=1, x=1 )
# Data
list(n=94)
# Data
e[] y[]
532 0
584 0
672 2
722 1
904 1
1236
0
950 0
1405
1
776 3
1013
0
739 0
1770
1
821 0
1115
2
1164
3
1164
0
1303
0
1774
3
3585
1
1193
1
1213
1
1232
1
1517
4
1520
3
1862
3
1888
1
1247
0
1381
2
1643
2
1660
4
1827
4
1486
3
1593 2
2265
4
1524
1
1759
3
1309
0
1529
4
1677
1
1654
2
1785
3
1979
4
1767
4
2465
2
1750
2
2458
4
2383
2
2717
3
2282
0
2115
0
2852
2
2856
5
3174
5
2369
1
2557
1
3859
3
2641
1
2741
1
3055
3
3513
1
2728
2
3354
6
3814
0
4014
2
2612
2
2815
1
4294
2
3450
8
3628
6
4219
1
3932
6
4082
4
4203
1
4022
3
4636
5
5571
2
6436
4
5344
3
4445
4
4705
5
5039
2
6043
6
5121
8
11260
5
5789
0
6044
6
5569
8
6130
7
6249
3
7002
3
7851
9
9573
7
12050 18
12131 17
END