如何使函数中的变量等于之前 jags 试验中该变量的总和
How to make a variable in the function equal to the sum of that variable in previous trials in jags
在我的模型中，每个主题和每个试验都有 10 个选项（从 1-10）可供选择 (expectancy)。我根据下图中的规则计算了每个选项的值，因此每个选项的值根据每次试验中 shock 和 v 之间的差异更新（乘以 alpha） .然后，我使用 softmax 规则将每个选项的 v 转换为具有相同函数的特定概率，在此威胁中：JAGS errors: "Resolving undeclared variables" and "Invalid vector argument to exp".
我想这里的问题是我无法让 jags 更新相同选择的值。
数据：expectancy = 每次试验中 1-10 的数字。 shock=每次试验中的数字为 1 或 0。 （我在下面提供了示例数据）
第二个情节是如何在 stan 中完成 2 choices/1 主题情况。
RW_model <- function(){
# data
for(i in 1:nsubjects) # for each person
{
# initial value for v
v [i,1,expectancy[i,1]] <- 0
for (j in 2:ntrials) # for each trial
{
# expectancy chosen
expectancy[i,j] ~ dcat(mu[i,j,1:10])
predk[i,j] ~ dcat(mu[i,j,1:10])
# softmax rule to calculate values of each expectancy for each subject
# tau is the value sensitivity parameter
mu[i,j,1:10] <- exp_v[i,j,1:10] / sum(exp_v[i,j,1:10])
exp_v[i,j,expectancy[i,j-1]] <- exp(v[i,j,expectancy[i,j-1]]/tau[i])
# prediction error: difference between feedback and learned values of the chosen expectancy
pe [i,j-1] <- shock [i,j-1] - v [i,j-1,expectancy[i,j-1]]
# value updating process for expectancy
v [i,j,expectancy[i,j-1]] <- v [i,j-1,expectancy[i,j-1]] + alpha [i] * pe [i,j-1]
}
}
# priors
for (i in 1:nsubjects){
tau [i] ~ dunif (0,3)
alpha [i] ~ dunif (0,1)
}
}
# example data/ initial value/ parameters
nsubjects <- 42
ntrials <- 14
shock <- matrix(c(0,0,1,1,0,0,1,1,0,0,1,0,1,0),nrow=42,ncol = 14,byrow = T)
expectancy <- matrix(c(1,2,3,4,5,6,7,7,8,8,7,10,10,00),nrow=42,ncol = 14,byrow = T)
data <- list('shock','nsubjects','ntrials','expectancy')
myinits <- list(list(tau = runif (42,0,3),
alpha = runif (42,0,1)))
parameters <- c("tau",'alpha','v','predk')
# jags sampling
samples <- jags(data, inits=myinits, parameters,
model.file = RW_model,
n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
由于设置中没有随机节点，我不确定 MCMC 模拟会给你带来什么，但这里有一些有效的代码。据我所知，主要问题是当 table 说 c(Vmax-Vn) 时它没有使用 c() 作为连接运算符。 c 是定义为 0.3 的常量。
dat <- list(
v = c(0, rep(NA, 8)),
vn = rep(NA, 8),
vmax=1,
c=0.3,
Ntrials = 9
)
jmod <- "model{
for(i in 2:Ntrials){
v[i] <- c*(vmax - v[(i-1)]) + v[(i-1)]
}
}"
library(runjags)
out <- run.jags(jmod, data=dat, monitor="v")
编辑：已编辑问题的答案。
你必须看看这是否在做你想要的，但它产生了一个结果。据我所知，这是将 Stan 代码直接扩展到多个主题和每个试验的多个选择移植到 JAGS。一、型号：
RW_model <- function(){
for(i in 1:Nsubjects){
for(j in 1:Ntrials){
expectancy[i,j] ~ dcat(p[i, 1:Nposs, j])
mu[i,1:Nposs, j] <- tau[i] * v[i,1:Nposs, j]
for(k in 1:Nposs){
q[i,k,j] <- exp(mu[i,k,j])
p[i,k,j] <- q[i,k,j]/sum(q[i,,j])
}
pe[i,j] <- shock[i,j] - v[i,expectancy[i,j], j]
for(k in 1:Nposs){
v[i,k,j+1] <- v[i,k,j] + ifelse(expectancy[i,j] == k,
alpha[i] * pe[i,j],
0)
}
}
}
for (i in 1:Nsubjects){
tau [i] ~ dunif (0,3)
alpha [i] ~ dunif (0,1)
}
}
接下来，我们可以补一些数据。这里的数据是针对 20 个受试者，expectancy 中的 5 个可能的选择和 14 个试验。
set.seed(40120)
# v has to be a Nsubjects x Nchoices x Ntrials+1 matrix
v <- array(NA, dim=c(20, 5, 15))
# The first trial of v is initialized to 0 fo all subjects and choices
v[,,1] <- matrix(0, nrow=20, ncol=5)
# expectancy and shock are both Nsubjects x Ntrials matrices
expectancy <- matrix(sample(1:5, 20*14, replace=TRUE), ncol=14)
shock <- matrix(sample(c(0,1), 20*14, replace=TRUE), ncol=14)
dlist <- list(
Nsubjects = 20,
Ntrials = 14,
Nposs = 5,
expectancy = expectancy,
shock = shock,
v=v
)
最后，我们可以指定初始值和 运行 模型。
myinits <- list(list(tau = runif (nrow(expectancy),0,3),
alpha = runif (nrow(expectancy)0,1)))
parameters <- c("tau",'alpha')
library(R2jags)
# jags sampling
samples <- jags(dlist, inits=myinits, parameters,
model.file = RW_model,
n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
您可以从输出中看到它生成了样本，并根据受试者估计了不同的 tau 和 alpha 参数。
samples$BUGSoutput
Inference for Bugs model at "/var/folders/55/q9y1hbcx13b5g50kks_p0mb00000gn/T//Rtmp32zeG4/model1726bdee0c99.txt", fit using jags,
1 chains, each with 1000 iterations (first 500 discarded)
n.sims = 500 iterations saved
mean sd 2.5% 25% 50% 75% 97.5%
alpha[1] 0.6 0.3 0.1 0.4 0.7 0.9 1.0
alpha[2] 0.4 0.3 0.0 0.1 0.3 0.7 0.9
alpha[3] 0.3 0.2 0.0 0.1 0.2 0.4 0.9
alpha[4] 0.4 0.3 0.0 0.1 0.3 0.6 0.9
alpha[5] 0.4 0.3 0.0 0.1 0.3 0.6 1.0
alpha[6] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[7] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[8] 0.3 0.2 0.0 0.1 0.2 0.5 0.9
alpha[9] 0.3 0.3 0.0 0.1 0.3 0.5 0.9
alpha[10] 0.4 0.2 0.0 0.2 0.4 0.5 0.9
alpha[11] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[12] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[13] 0.3 0.3 0.0 0.1 0.3 0.5 0.9
alpha[14] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[15] 0.5 0.3 0.0 0.2 0.5 0.7 1.0
alpha[16] 0.4 0.3 0.0 0.1 0.3 0.7 0.9
alpha[17] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[18] 0.4 0.3 0.0 0.2 0.4 0.7 1.0
alpha[19] 0.4 0.3 0.0 0.2 0.4 0.7 1.0
alpha[20] 0.3 0.3 0.0 0.0 0.2 0.5 0.9
deviance 909.3 5.1 901.1 905.8 908.8 912.5 919.9
tau[1] 1.1 0.6 0.2 0.7 1.1 1.5 2.3
tau[2] 0.8 0.7 0.0 0.3 0.6 1.2 2.7
tau[3] 1.0 0.8 0.0 0.3 0.8 1.6 2.8
tau[4] 1.1 0.7 0.0 0.4 0.9 1.6 2.7
tau[5] 0.9 0.7 0.0 0.3 0.7 1.3 2.7
tau[6] 1.0 0.8 0.0 0.3 0.9 1.6 2.8
tau[7] 1.0 0.7 0.1 0.5 0.9 1.5 2.7
tau[8] 1.1 0.8 0.0 0.4 0.9 1.7 2.9
tau[9] 1.0 0.8 0.0 0.3 0.8 1.6 2.7
tau[10] 1.6 0.8 0.1 0.9 1.7 2.3 2.9
tau[11] 1.1 0.9 0.0 0.3 0.9 1.8 2.8
tau[12] 1.0 0.8 0.0 0.4 0.8 1.6 2.9
tau[13] 0.9 0.8 0.0 0.3 0.6 1.4 2.7
tau[14] 1.1 0.8 0.0 0.5 0.9 1.6 2.7
tau[15] 1.2 0.7 0.1 0.7 1.1 1.6 2.7
tau[16] 1.0 0.8 0.1 0.4 0.9 1.5 2.8
tau[17] 1.1 0.8 0.1 0.4 0.8 1.6 2.9
tau[18] 1.2 0.8 0.1 0.5 1.1 1.7 2.7
tau[19] 1.2 0.8 0.1 0.5 1.0 1.8 2.8
tau[20] 1.0 0.9 0.0 0.3 0.8 1.5 2.9
DIC info (using the rule, pD = var(deviance)/2)
pD = 12.9 and DIC = 922.1
DIC is an estimate of expected predictive error (lower deviance is better).
由于我不知道该模型会带来什么，您必须验证它是否确实按照您的预期进行，但看起来这是一个很好的起点。
在我的模型中，每个主题和每个试验都有 10 个选项（从 1-10）可供选择 (expectancy)。我根据下图中的规则计算了每个选项的值，因此每个选项的值根据每次试验中 shock 和 v 之间的差异更新（乘以 alpha） .然后，我使用 softmax 规则将每个选项的 v 转换为具有相同函数的特定概率，在此威胁中：JAGS errors: "Resolving undeclared variables" and "Invalid vector argument to exp".
我想这里的问题是我无法让 jags 更新相同选择的值。
数据：expectancy = 每次试验中 1-10 的数字。 shock=每次试验中的数字为 1 或 0。 （我在下面提供了示例数据）
第二个情节是如何在 stan 中完成 2 choices/1 主题情况。
RW_model <- function(){
# data
for(i in 1:nsubjects) # for each person
{
# initial value for v
v [i,1,expectancy[i,1]] <- 0
for (j in 2:ntrials) # for each trial
{
# expectancy chosen
expectancy[i,j] ~ dcat(mu[i,j,1:10])
predk[i,j] ~ dcat(mu[i,j,1:10])
# softmax rule to calculate values of each expectancy for each subject
# tau is the value sensitivity parameter
mu[i,j,1:10] <- exp_v[i,j,1:10] / sum(exp_v[i,j,1:10])
exp_v[i,j,expectancy[i,j-1]] <- exp(v[i,j,expectancy[i,j-1]]/tau[i])
# prediction error: difference between feedback and learned values of the chosen expectancy
pe [i,j-1] <- shock [i,j-1] - v [i,j-1,expectancy[i,j-1]]
# value updating process for expectancy
v [i,j,expectancy[i,j-1]] <- v [i,j-1,expectancy[i,j-1]] + alpha [i] * pe [i,j-1]
}
}
# priors
for (i in 1:nsubjects){
tau [i] ~ dunif (0,3)
alpha [i] ~ dunif (0,1)
}
}
# example data/ initial value/ parameters
nsubjects <- 42
ntrials <- 14
shock <- matrix(c(0,0,1,1,0,0,1,1,0,0,1,0,1,0),nrow=42,ncol = 14,byrow = T)
expectancy <- matrix(c(1,2,3,4,5,6,7,7,8,8,7,10,10,00),nrow=42,ncol = 14,byrow = T)
data <- list('shock','nsubjects','ntrials','expectancy')
myinits <- list(list(tau = runif (42,0,3),
alpha = runif (42,0,1)))
parameters <- c("tau",'alpha','v','predk')
# jags sampling
samples <- jags(data, inits=myinits, parameters,
model.file = RW_model,
n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
由于设置中没有随机节点，我不确定 MCMC 模拟会给你带来什么，但这里有一些有效的代码。据我所知，主要问题是当 table 说 c(Vmax-Vn) 时它没有使用 c() 作为连接运算符。 c 是定义为 0.3 的常量。
dat <- list(
v = c(0, rep(NA, 8)),
vn = rep(NA, 8),
vmax=1,
c=0.3,
Ntrials = 9
)
jmod <- "model{
for(i in 2:Ntrials){
v[i] <- c*(vmax - v[(i-1)]) + v[(i-1)]
}
}"
library(runjags)
out <- run.jags(jmod, data=dat, monitor="v")
编辑：已编辑问题的答案。
你必须看看这是否在做你想要的，但它产生了一个结果。据我所知，这是将 Stan 代码直接扩展到多个主题和每个试验的多个选择移植到 JAGS。一、型号：
RW_model <- function(){
for(i in 1:Nsubjects){
for(j in 1:Ntrials){
expectancy[i,j] ~ dcat(p[i, 1:Nposs, j])
mu[i,1:Nposs, j] <- tau[i] * v[i,1:Nposs, j]
for(k in 1:Nposs){
q[i,k,j] <- exp(mu[i,k,j])
p[i,k,j] <- q[i,k,j]/sum(q[i,,j])
}
pe[i,j] <- shock[i,j] - v[i,expectancy[i,j], j]
for(k in 1:Nposs){
v[i,k,j+1] <- v[i,k,j] + ifelse(expectancy[i,j] == k,
alpha[i] * pe[i,j],
0)
}
}
}
for (i in 1:Nsubjects){
tau [i] ~ dunif (0,3)
alpha [i] ~ dunif (0,1)
}
}
接下来，我们可以补一些数据。这里的数据是针对 20 个受试者，expectancy 中的 5 个可能的选择和 14 个试验。
set.seed(40120)
# v has to be a Nsubjects x Nchoices x Ntrials+1 matrix
v <- array(NA, dim=c(20, 5, 15))
# The first trial of v is initialized to 0 fo all subjects and choices
v[,,1] <- matrix(0, nrow=20, ncol=5)
# expectancy and shock are both Nsubjects x Ntrials matrices
expectancy <- matrix(sample(1:5, 20*14, replace=TRUE), ncol=14)
shock <- matrix(sample(c(0,1), 20*14, replace=TRUE), ncol=14)
dlist <- list(
Nsubjects = 20,
Ntrials = 14,
Nposs = 5,
expectancy = expectancy,
shock = shock,
v=v
)
最后，我们可以指定初始值和 运行 模型。
myinits <- list(list(tau = runif (nrow(expectancy),0,3),
alpha = runif (nrow(expectancy)0,1)))
parameters <- c("tau",'alpha')
library(R2jags)
# jags sampling
samples <- jags(dlist, inits=myinits, parameters,
model.file = RW_model,
n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
您可以从输出中看到它生成了样本，并根据受试者估计了不同的 tau 和 alpha 参数。
samples$BUGSoutput
Inference for Bugs model at "/var/folders/55/q9y1hbcx13b5g50kks_p0mb00000gn/T//Rtmp32zeG4/model1726bdee0c99.txt", fit using jags,
1 chains, each with 1000 iterations (first 500 discarded)
n.sims = 500 iterations saved
mean sd 2.5% 25% 50% 75% 97.5%
alpha[1] 0.6 0.3 0.1 0.4 0.7 0.9 1.0
alpha[2] 0.4 0.3 0.0 0.1 0.3 0.7 0.9
alpha[3] 0.3 0.2 0.0 0.1 0.2 0.4 0.9
alpha[4] 0.4 0.3 0.0 0.1 0.3 0.6 0.9
alpha[5] 0.4 0.3 0.0 0.1 0.3 0.6 1.0
alpha[6] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[7] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[8] 0.3 0.2 0.0 0.1 0.2 0.5 0.9
alpha[9] 0.3 0.3 0.0 0.1 0.3 0.5 0.9
alpha[10] 0.4 0.2 0.0 0.2 0.4 0.5 0.9
alpha[11] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[12] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[13] 0.3 0.3 0.0 0.1 0.3 0.5 0.9
alpha[14] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[15] 0.5 0.3 0.0 0.2 0.5 0.7 1.0
alpha[16] 0.4 0.3 0.0 0.1 0.3 0.7 0.9
alpha[17] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[18] 0.4 0.3 0.0 0.2 0.4 0.7 1.0
alpha[19] 0.4 0.3 0.0 0.2 0.4 0.7 1.0
alpha[20] 0.3 0.3 0.0 0.0 0.2 0.5 0.9
deviance 909.3 5.1 901.1 905.8 908.8 912.5 919.9
tau[1] 1.1 0.6 0.2 0.7 1.1 1.5 2.3
tau[2] 0.8 0.7 0.0 0.3 0.6 1.2 2.7
tau[3] 1.0 0.8 0.0 0.3 0.8 1.6 2.8
tau[4] 1.1 0.7 0.0 0.4 0.9 1.6 2.7
tau[5] 0.9 0.7 0.0 0.3 0.7 1.3 2.7
tau[6] 1.0 0.8 0.0 0.3 0.9 1.6 2.8
tau[7] 1.0 0.7 0.1 0.5 0.9 1.5 2.7
tau[8] 1.1 0.8 0.0 0.4 0.9 1.7 2.9
tau[9] 1.0 0.8 0.0 0.3 0.8 1.6 2.7
tau[10] 1.6 0.8 0.1 0.9 1.7 2.3 2.9
tau[11] 1.1 0.9 0.0 0.3 0.9 1.8 2.8
tau[12] 1.0 0.8 0.0 0.4 0.8 1.6 2.9
tau[13] 0.9 0.8 0.0 0.3 0.6 1.4 2.7
tau[14] 1.1 0.8 0.0 0.5 0.9 1.6 2.7
tau[15] 1.2 0.7 0.1 0.7 1.1 1.6 2.7
tau[16] 1.0 0.8 0.1 0.4 0.9 1.5 2.8
tau[17] 1.1 0.8 0.1 0.4 0.8 1.6 2.9
tau[18] 1.2 0.8 0.1 0.5 1.1 1.7 2.7
tau[19] 1.2 0.8 0.1 0.5 1.0 1.8 2.8
tau[20] 1.0 0.9 0.0 0.3 0.8 1.5 2.9
DIC info (using the rule, pD = var(deviance)/2)
pD = 12.9 and DIC = 922.1
DIC is an estimate of expected predictive error (lower deviance is better).
由于我不知道该模型会带来什么，您必须验证它是否确实按照您的预期进行，但看起来这是一个很好的起点。