未加预测:
model {
for (i in 1:n){
Y[i]~dnorm(s[i],tauC)
s[i]<-theta[1]/(1+theta[2]*exp(-theta[3]*t[i]))}
theta[1]~dnorm(mu[1],tau[1])
theta[2]~dnorm(mu[2],tau[2])
theta[3]~dnorm(mu[3],tau[3])
tauC~dgamma(0.01,0.01)
sigmaC<-1/sqrt(tauC)
for (k in 1:3) {
mu[k]~dnorm(0,1.0E-6)
tau[k]~dgamma(0.01,0.01)}
}
#Data
list(n=20,t=c(5,10,20,21,25,30,34,35,40,45,48,50,55,60,62,65,70,75,76,80),
Y=c(0.51,1.44,2.71,3.21,3.57,5.35,5.88,6.56,7.48,7.82,8.05,8.38,8.59,8.84,9.01,9.18,9.81,10.04,10.06,10.11)
)
#Init
list(theta=c(1,1,1),mu=c(1,1,1),tau=c(1,1,1),tauC=1)
不加预测
node mean sd MC error 2.5% median 97.5% start sample
sigmaC 0.3456 0.06482 0.001814 0.2473 0.3368 0.4981 4001 7000
theta[1] 9.818 0.1768 0.009688 9.469 9.815 10.18 4001 7000
theta[2] 16.01 3.108 0.2881 11.58 15.34 24.48 4001 7000
theta[3] 0.09341 0.00665 6.033E-4 0.08179 0.09268 0.1097 4001 7000
加预测后:
model {
for (i in 1:n){
Y[i]~dnorm(s[i],tauC)
s[i]<-theta[1]/(1+theta[2]*exp(-theta[3]*t[i]))}
theta[1]~dnorm(mu[1],tau[1])
theta[2]~dnorm(mu[2],tau[2])
theta[3]~dnorm(mu[3],tau[3])
tauC~dgamma(0.01,0.01)
sigmaC<-1/sqrt(tauC)
for (k in 1:3) {
mu[k]~dnorm(0,1.0E-6)
tau[k]~dgamma(0.01,0.01)}
for(m in 1:20){
pred.Y[m]~dnorm(pred.s[m],tauC)
pred.s[m]<-theta[1]/(1+theta[2]*exp(-theta[3]*t.new[m]))
}
}
#Data
list(n=20,t=c(5,10,20,21,25,30,34,35,40,45,48,50,55,60,62,65,70,75,76,80),
Y=c(0.51,1.44,2.71,3.21,3.57,5.35,5.88,6.56,7.48,7.82,8.05,8.38,8.59,8.84,9.01,9.18,9.81,10.04,10.06,10.11),
t.new=c(5,10,20,21,25,30,34,35,40,45,48,50,55,60,62,65,70,75,76,80),
)
#Init
list(theta=c(1,1,1),mu=c(1,1,1),tau=c(1,1,1),tauC=1,
pred.Y=c(0.51,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA))
MCMC运行的结果:(预测的结果是错的啊)
node mean sd MC error 2.5% median 97.5% start sample
pred.Y[1] 6.59 3.43 0.1568 -0.5105 6.707 13.07 4001 6000
pred.Y[2] 6.608 3.462 0.156 -0.4805 6.652 13.38 4001 6000
pred.Y[3] 6.588 3.365 0.1369 -0.04587 6.667 13.15 4001 6000
pred.Y[4] 6.65 3.334 0.1293 0.1088 6.672 13.21 4001 6000
pred.Y[5] 6.632 3.24 0.07678 0.3863 6.647 13.12 4001 6000
pred.Y[6] 6.784 3.235 0.03608 0.3023 6.807 13.24 4001 6000
pred.Y[7] 6.896 3.222 0.05258 0.4493 6.982 13.14 4001 6000
pred.Y[8] 6.872 3.235 0.04595 0.3636 6.893 13.25 4001 6000
pred.Y[9] 6.9 3.231 0.05063 0.4337 6.967 13.15 4001 6000
pred.Y[10] 6.881 3.234 0.05925 0.5227 6.979 13.2 4001 6000
pred.Y[11] 6.876 3.22 0.05991 0.2422 6.951 13.24 4001 6000
pred.Y[12] 6.959 3.223 0.05387 0.5343 7.057 13.12 4001 6000
pred.Y[13] 6.952 3.202 0.0574 0.5909 6.993 13.19 4001 6000
pred.Y[14] 6.976 3.238 0.06121 0.3893 7.093 13.28 4001 6000
pred.Y[15] 6.853 3.222 0.06024 0.4255 6.906 12.99 4001 6000
pred.Y[16] 6.897 3.23 0.06107 0.3054 6.967 13.09 4001 6000
pred.Y[17] 6.95 3.255 0.06438 0.4199 7.008 13.58 4001 6000
pred.Y[18] 6.898 3.29 0.05809 0.3202 6.983 13.25 4001 6000
pred.Y[19] 6.842 3.267 0.05323 0.401 6.875 13.25 4001 6000
pred.Y[20] 6.965 3.256 0.05402 0.4672 7.07 13.41 4001 6000
sigmaC 3.084 0.6551 0.0428 1.287 3.062 4.406 4001 6000
theta[1] 6.921 0.7761 0.04164 5.473 6.876 8.649 4001 6000
theta[2] 112.5 1021.0 111.6 -1967.0 88.81 2557.0 4001 6000
theta[3] 805.8 674.8 72.01 0.3036 703.8 2365.0 4001 6000