问题描述
我正在尝试进行模型拟合,并从Stan开始,但是已经迁移到JAGS尝试进行一些预测,并且遇到了一些奇怪的边界影响,这些影响是我没能进入Stan的参数估计值。
我在Stan中所做的工作与在JAGS中所做的工作之间的主要区别在于:(1)在Stan中,我使用的是RK ODE求解器,而在JAGS中,我必须使用欧拉逼近法; (2)在Stan中,我正在进行全局拟合,而在JAGS中,我正在一步一步拟合。
这是我的JAGS模型脚本:
# Functional response project
## One-step-ahead model with JAGS instead of Stan because I'm stupid
# data {
# for (i in 2:ts) {
# New[i] <- y[i] - y[i-1]
# }
# }
model { # Keep priors simple for Now
r~dunif(0,5) # Intrinsic growth rate of parasite
O~dunif(0,1) # Rate at which immune system correctly recognizes parasites
h~dunif(0,1) # Handelling time
b~dunif(0,100) # Immigration rate of immune cells to infection site
c~dunif(0,1) # Growth rate of immune system as a response to infection
u~dunif(0,1) # Natural mortality rate of immune cells
#r~dunif(2,3) # Intrinsic growth rate of parasite
#O~dunif(0,0.5) # Rate at which immune system correctly recognizes parasites
#h~dunif(0,0.5) # Handelling time
#b~dunif(30,40) # Immigration rate of immune cells to infection site
#c~dunif(0,1) # Growth rate of immune system as a response to infection
#u~dunif(0,1) # Natural mortality rate of immune cells
# State vars
P[1] <- 80
H[1] <- 200 # Dat
y_hat_P[1] <- 80
y_hat_H[1] <- 200 # Est
for (i in 2:L_T){
P[i] <- P[i-1] + (P[i-1]*r - H[i-1]*(O*P[i-1]/(1 + O*P[i-1]*h)))*dt
H[i] <- H[i-1] + (b + H[i-1]*(c*(O*P[i-1]/(1 + O*P[i-1]*h)) - u))*dt
y_hat_P[i] <- P[i];
y_hat_H[i] <- H[i]; # Expectation
C_P[i] ~ dpois(y_hat_P[i]);
C_H[i] ~ dpois(y_hat_H[i]);
}
}
这是我的实现:
# Data org.
IC=list(r=2.5,O=0.012,h=0.075,b=35,c=0.3,u=0.41); inits=list(IC,IC,IC)
data = SmallTS_Data
L_T = length(data$t)
C_P = as.integer(data$P); C_H = as.integer(data$H)
dt <- 0.001 # Same for all subsequnt models; small step size for "better" Euler approx.
# First trial
summary(Fit7)
plot(Fit7)
# Trials w/out setting seed
# Also,increased iter by 10x
Fit7.2 = run.jags(#model = here('./Current_Code/OSA_JAGS_Para.R'),model = "/Users/mjarviscross/Desktop/GitHub/FcnalRespProj_2/Current_Code/OSA_JAGS_SmallTS.R",data = list('L_T' = L_T,'C_P' = C_P,'C_H' = C_H,'dt' = dt),monitor = c('r','O','h','b','c','u'),n.chains = 4,method = 'rjags',burnin = 100000,adapt = 100000,sample = 10000,inits = inits,summarise = TRUE,plots = TRUE)
summary(Fit7.2)
plot(Fit7.2)
我尝试混合使用“ inits”,因此每条链都从一个不同的地方开始,没有运气。在上面的块中,init值是参数的真实值。
我还尝试记录/求幂一些参数,以使参数接近0的采样分布更容易,但这也无济于事。
奇怪的是,其中4/6个参数遇到边界效应(r,h,b和u),而O和c没有。
这是输出:
Lower95 Median Upper95 Mean SD Mode
r 4.977625e+00 4.994850857 4.99999888 4.992436257 0.0077408904 4.997662853
O 2.951738e-02 0.029777153 0.03030740 0.029825597 0.0002001345 0.029726595
h 1.283357e-07 0.002377487 0.01023549 0.003333945 0.0030660994 0.001091999
b 9.604828e+01 98.843977266 99.99991291 98.522219100 1.2420867280 99.477325634
c 7.300956e-02 0.082557460 0.09259498 0.082891674 0.0048552133 0.082253983
u 9.845811e-01 0.996027921 0.99999981 0.994573133 0.0049174393 0.998227142
MCerr MC%ofSD SSeff AC.10 psrf
r 3.120977e-04 4.0 615 0.2528493 1.016287
O 2.513034e-05 12.6 63 0.8550039 1.074630
h 4.088013e-04 13.3 56 0.8525586 1.021393
b 8.206921e-02 6.6 229 0.5700555 1.037147
c 3.548487e-04 7.3 187 0.6381600 1.051069
u 2.474836e-04 5.0 395 0.3892634 1.015181
任何帮助将不胜感激! :)
P.S。如果有人对将这种一步一步的模型转换为Stan有任何建议,我很想听听!我已经在Stan论坛上花了很短的时间运气了!
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)