model {
for (i in 1:N) {
Y[i] ~ dnorm(mu, tau)
Y.rep[i] ~ dnorm(mu, tau)
T.data.obs[i] <- pow((Y[i] - mean(Y[]))/sd(Y[]),3)
T.data.rep[i] <- pow((Y.rep[i] - mean(Y.rep[]))/sd(Y.rep[]),3)
T.para.obs[i] <- pow((Y[i] - mu)/sigma,3)
T.para.rep[i] <- pow((Y.rep[i] - mu)/sigma,3)
}
mu ~ dunif(-100, 100)
tau ~ dgamma(0.001, 0.001)
sigma <- 1/sqrt(tau)
T.data.obs.tot <- sum(T.data.obs[])
T.data.rep.tot <- sum(T.data.rep[])
T.para.obs.tot <- sum(T.para.obs[])
T.para.rep.tot <- sum(T.para.rep[])
P.data <- step(T.data.obs.tot - T.data.rep.tot)
P.para <- step(T.para.obs.tot - T.para.rep.tot)
}
Inits:
list(mu = 0, tau = 1)
Data:
list(N=66, Y=c(
28, 26, 33, 24, 34, -44, 27, 16, 40, -2,
29, 22, 24, 21, 25, 30, 23, 29, 31, 19,
24, 20, 36, 32, 36, 28, 25, 21, 28, 29,
37, 25, 28, 26, 30, 32, 36, 26, 30, 22,
36, 23, 27, 27, 28, 27, 31, 27, 26, 33,
26, 32, 32, 24, 39, 28, 24, 25, 32, 25,
29, 27, 28, 29, 16, 23))
node mean sd MC error 2.5% median 97.5% start sample
P.data 0.0 0.0 1.0E-12 0.0 0.0 0.0 1 10000
P.para 0.0 0.0 1.0E-12 0.0 0.0 0.0 1 10000
T.data.obs.tot -289.8 0.0 1.0E-12 -289.8 -289.8 -289.8 1 10000
T.data.rep.tot 0.02773 18.65 0.1727 -37.1 -0.01383 37.07 1 10000
T.para.obs.tot -293.3 81.25 0.8554 -475.1 -285.5 -157.4 1 10000
T.para.rep.tot -0.444 31.57 0.2758 -62.93 -0.4164 62.88 1 10000
mu 26.23 1.34 0.01353 23.6 26.23 28.88 1 10000
sigma 10.87 0.9811 0.01047 9.169 10.79 13.0 1 10000
tau 0.00866 0.001532 1.653E-5 0.005921 0.008591 0.0119 1 10000