R：return 产品总和的有效方法（与 JAGS 兼容）
R: Efficient way to return sums of products (compatible with JAGS)
我正在尝试计算数组中由乘积之和产生的单元格概率。这是我将在 R 中的 JAGS 中 运行 模型的一部分。在此示例中，对象可能处于 3 种状态之一。在每个时间间隔，它可能保持相同状态或转到其他状态之一。使用下面的示例代码，起始状态为 1，最终状态为 2。数组 X 是状态*状态转换矩阵（例如，X[1,2] 是状态 1 中的对象在时间 t 处于状态的概率2 在时间 t+1)。请注意，这些只是用来测试 test1 和 test2（见下文）给出相同结果的数字，因为我确信手算（test1）是正确的。
X <- array((0),dim=c(3,3))
X[1,1] <- 0.1
X[1,2] <- 0.2
X[1,3] <- 0.3
X[2,1] <- 0.4
X[2,2] <- 0.5
X[2,3] <- 0.6
X[3,1] <- 0.7
X[3,2] <- 0.8
X[3,3] <- 0.9
示例已简化。我最终将为每个时间间隔有一个单独的状态*状态转换矩阵，因为转换是时间相关的。为简单起见，我还从示例中删除了观察概率。考虑到观察的概率（因为它最终需要），每个单元格描述了在该单元格的状态和时间中第一次观察到对象的概率 。 X 的每个元素都包含一个生存参数（对于时间 t 的状态 i）和转换参数（对于状态 i 在 t-1 到状态 j 在时间 t）。我将在模型中分别导出这些组件，我认为这部分不会太难。
回到下面给出的例子。当释放状态（样本场合 1 之前的时间间隔）为 1 时，我已经手动（测试 1）计算了状态 2 在样本场合 4 的概率。此计算包括到达状态 2 的所有可能方式的乘积之和在时间间隔内从状态 1 开始。
test1 <-
(X[1,1]*X[1,1]*X[1,1]*X[1,2])+
(X[1,1]*X[1,1]*X[1,2]*X[2,2])+
(X[1,1]*X[1,1]*X[1,3]*X[3,2])+
(X[1,1]*X[1,2]*X[2,1]*X[1,2])+
(X[1,1]*X[1,2]*X[2,2]*X[2,2])+
(X[1,1]*X[1,2]*X[2,3]*X[3,2])+
(X[1,1]*X[1,3]*X[3,1]*X[1,2])+
(X[1,1]*X[1,3]*X[3,2]*X[2,2])+
(X[1,1]*X[1,3]*X[3,3]*X[3,2])+
(X[1,2]*X[2,1]*X[1,1]*X[1,2])+
(X[1,2]*X[2,1]*X[1,2]*X[2,2])+
(X[1,2]*X[2,1]*X[1,3]*X[3,2])+
(X[1,2]*X[2,2]*X[2,1]*X[1,2])+
(X[1,2]*X[2,2]*X[2,2]*X[2,2])+
(X[1,2]*X[2,2]*X[2,3]*X[3,2])+
(X[1,2]*X[2,3]*X[3,1]*X[1,2])+
(X[1,2]*X[2,3]*X[3,2]*X[2,2])+
(X[1,2]*X[2,3]*X[3,3]*X[3,2])+
(X[1,3]*X[3,1]*X[1,1]*X[1,2])+
(X[1,3]*X[3,1]*X[1,2]*X[2,2])+
(X[1,3]*X[3,1]*X[1,3]*X[3,2])+
(X[1,3]*X[3,2]*X[2,1]*X[1,2])+
(X[1,3]*X[3,2]*X[2,2]*X[2,2])+
(X[1,3]*X[3,2]*X[2,3]*X[3,2])+
(X[1,3]*X[3,3]*X[3,1]*X[1,2])+
(X[1,3]*X[3,3]*X[3,2]*X[2,2])+
(X[1,3]*X[3,3]*X[3,3]*X[3,2])
test1 # = 0.9288
我写了一个 for 循环来简化每个释放状态和重新捕获状态的计算 (test2)。问题是，在更多的时间间隔内，我的方法将最终计算一个具有很多维度的数组的乘积总和，我担心它在 JAGS 中的计算量会非常大。
test2 <- array((0),dim=c(3,3,3,3,3))
for (m in 1:3){
for (i in 1:3){
for (j in 1:3){
for (k in 1:3){
for (s in 1:3){
test2[m,i,j,k,s] <- X[m,i]*X[i,j]*X[j,k]*X[k,s]
}
}
}
}
}
sumtest2 <- sum(test2[1,,,,2]) # probability of state 2 at observation
# time 4 after release in state 1 at
# release time 1.
test2
sumtest2 # = 0.9288 (the same as test1 and therefore the correct result)
我的问题是实现相同目标的更有效方法是什么？
谢谢。
--------------------------------
感谢 Dex 提供的非常有用的答案。我有一个跟进问题，关于时间相关的转移矩阵。上面的代码效果很好，但我无法将它用于不同次数的迭代，其中每次迭代的转换矩阵都不同。
对于每个释放时机 t 有 T-t 重获时机（其中 T 是总学习年数）。我想用矩阵乘以匹配第一次可能夺回时间到最后一次夺回时间的矩阵。
这是我想要实现的长版的一部分：
# This long version does work, with lines for release occasion 1
# and recapture occasions 1:4 only:
pr2 <- array((0), dim=c(3,3,4,4))
pr2[,,1,1] <- diag(1, 3, 3) %*% (X[,,1])
pr2[,,1,2] <- diag(1, 3, 3) %*% (X[,,1] %*% X[,,2])
pr2[,,1,3] <- diag(1, 3, 3) %*% (X[,,1] %*% X[,,2] %*% X[,,3])
pr2[,,1,4] <- diag(1, 3, 3) %*% (X[,,1] %*% X[,,2] %*% X[,,3] %*% X[,,4])
我试过使用额外的时间索引，然后使用带有矩阵乘法的 for 循环：
# With time dependence:
X <- array((0),dim=c(3,3,4))
X[1,1,1] <- 0.1
X[1,2,1] <- 0.2
X[1,3,1] <- 0.7
X[2,1,1] <- 0.25
X[2,2,1] <- 0.35
X[2,3,1] <- 0.4
X[3,1,1] <- 0.15
X[3,2,1] <- 0.25
X[3,3,1] <- 0.6
X[,,2] <- X[,,1]+0.05
X[,,3] <- X[,,1]-0.05
X[,,4] <- X[,,1]+0.1
# not working:
library(expm)
for (t in 1:4){
for (j in t:4){
pr1[,,t,j] <- diag(1, 3, 3) %*% (X[,,j]%^%j)
}
}
# This also doesn't work:
library(expm)
for (t in 1:4){
for (j in t:4){
pr1[,,t,j] <- diag(1, 3, 3) %*% (X[,,j]%^%1:j)
}
}
我也尝试过列表方法：
# List approach
# With time dependence:
X1 <- array((0),dim=c(3,3))
X1[1,1] <- 0.1
X1[1,2] <- 0.2
X1[1,3] <- 0.7
X1[2,1] <- 0.25
X1[2,2] <- 0.35
X1[2,3] <- 0.4
X1[3,1] <- 0.15
X1[3,2] <- 0.25
X1[3,3] <- 0.6
X2 <- X1+0.05
X3 <- X1-0.05
X4 <- X1+0.1
XT <- list(X1,X2,X3,X4)
pr1 <- array((0), dim=c(3,3,4,4))
# This also doesn't work:
library(expm)
for (t in 1:4){
for (j in t:4){
pr1[,,t,j] <- diag(1, 3, 3) %*% (XT[[1:j]]%^%j)
}
}
以下代码显示了计算时间相关数组中概率的很长的方法。它给出了正确的值。完整地写出来（有 23 次释放事件和 23 次重新捕获事件，对于 3 个 release/recapture 状态），虽然可能看起来效率不高，特别是如果分析需要稍后调整。这仅显示了在第一次释放事件中首次释放的个人的前四次重新捕获事件。每个单元格给出了一个人在那个场合 第一次 被重新捕获的概率。重新捕获概率由 p 给出，因此未重新捕获的概率由 1-p 给出。在我的分析中，T 和 p 是要估计的参数。此处定义了它们的值，以确保机制按预期工作。
# Make a 3D transition matrix T where z-dimension is time
T_states <- 3
T_terms <- 4
Tz <- 4
T <- array(runif(T_states * T_states * Tz)/1.5, dim = c(T_states, T_states, Tz))
# Normalise by rows - this step ensures rows sum to 1
for (z in seq(Tz)) {
T[, , z] <- T[, , z] / apply(T[, , z], 1, sum)
}
# state and time dependent recapture probability array [,recapture state,occasion]
# the first dimension is needed to repeat p across release states
p <- array((0),dim=c(3,3,4))
p[,1,1] <- 0.5
p[,1,2] <- 0.55
p[,1,3] <- 0.6
p[,1,4] <- 0.65
p[,2,1] <- 0.6
p[,2,2] <- 0.65
p[,2,3] <- 0.7
p[,2,4] <- 0.75
p[,3,1] <- 0.7
p[,3,2] <- 0.75
p[,3,3] <- 0.8
p[,3,4] <- 0.85
# Put them together (with 1-p for earlier recapture occasions):
pr2 <- array((0), dim=c(3,3,4,4))
pr2[,,1,1] <- diag(1, 3, 3) %*% (T[,,1]*p[,,1])
pr2[,,1,2] <- diag(1, 3, 3) %*% ((T[,,1]*(1-p[,,1])) %*% (T[,,2]*p[,,2]))
pr2[,,1,3] <- diag(1, 3, 3) %*% ((T[,,1]*(1-p[,,1])) %*% (T[,,2]*(1-p[,,2])) %*% (T[,,3]*p[,,1]))
pr2[,,1,4] <- diag(1, 3, 3) %*% ((T[,,1]*(1-p[,,1])) %*% (T[,,2]*(1-p[,,2])) %*% (T[,,3]*(1-p[,,3])) %*% (T[,,4]*p[,,4]))
pr2[,,1,]
我一直没能找到一种方法来为此编写更高效的代码来提供相同的值（很多方法都提供不正确的值）。
你可以用矩阵乘法来解决这个问题。转移矩阵可以矩阵相乘在一起，以在多次迭代中给出一个转移矩阵。对于四次迭代，我们执行四次：
(X %*% X %*% X %*% X)
您还可以使用 expm 包中的 %^% 运算符来进行矩阵求幂，而无需将其写出。
要找出对象在四次迭代后从状态 1 开始到状态 2 结束的可能性，我们可以将其表示为 1x3 矩阵并将其乘以转移矩阵四次。
matrix(c(1, 0, 0), nrow = 1, ncol = 3) %*% (X %*% X %*% X %*% X)
这给出了它处于三种状态的概率：
[,1] [,2] [,3]
[1,] 0.756 0.9288 1.1016
时机
我根据您的代码制作了另一个示例，其中的转换矩阵具有合理的概率。矩阵乘法比 for 循环快大约 100 倍并且给出相同的答案。
library("microbenchmark")
X <- array((0),dim=c(3,3))
X[1,1] <- 0.1
X[1,2] <- 0.2
X[1,3] <- 0.7
X[2,1] <- 0.25
X[2,2] <- 0.35
X[2,3] <- 0.4
X[3,1] <- 0.15
X[3,2] <- 0.25
X[3,3] <- 0.6
microbenchmark(
{test2 <- array((0),dim=c(3,3,3,3,3))
for (m in 1:3){
for (i in 1:3){
for (j in 1:3){
for (k in 1:3){
for (s in 1:3){
test2[m,i,j,k,s] <- X[m,i]*X[i,j]*X[j,k]*X[k,s]
}
}
}
}
}
sumtest2 <- sum(test2[1,,,,2])},
matrix(c(1, 0, 0), nrow = 1, ncol = 3) %*% X %*% X %*% X %*% X
)
我正在尝试计算数组中由乘积之和产生的单元格概率。这是我将在 R 中的 JAGS 中 运行 模型的一部分。在此示例中，对象可能处于 3 种状态之一。在每个时间间隔，它可能保持相同状态或转到其他状态之一。使用下面的示例代码，起始状态为 1，最终状态为 2。数组 X 是状态*状态转换矩阵（例如，X[1,2] 是状态 1 中的对象在时间 t 处于状态的概率2 在时间 t+1)。请注意，这些只是用来测试 test1 和 test2（见下文）给出相同结果的数字，因为我确信手算（test1）是正确的。
X <- array((0),dim=c(3,3))
X[1,1] <- 0.1
X[1,2] <- 0.2
X[1,3] <- 0.3
X[2,1] <- 0.4
X[2,2] <- 0.5
X[2,3] <- 0.6
X[3,1] <- 0.7
X[3,2] <- 0.8
X[3,3] <- 0.9
示例已简化。我最终将为每个时间间隔有一个单独的状态*状态转换矩阵，因为转换是时间相关的。为简单起见，我还从示例中删除了观察概率。考虑到观察的概率（因为它最终需要），每个单元格描述了在该单元格的状态和时间中第一次观察到对象的概率 。 X 的每个元素都包含一个生存参数（对于时间 t 的状态 i）和转换参数（对于状态 i 在 t-1 到状态 j 在时间 t）。我将在模型中分别导出这些组件，我认为这部分不会太难。
回到下面给出的例子。当释放状态（样本场合 1 之前的时间间隔）为 1 时，我已经手动（测试 1）计算了状态 2 在样本场合 4 的概率。此计算包括到达状态 2 的所有可能方式的乘积之和在时间间隔内从状态 1 开始。
test1 <-
(X[1,1]*X[1,1]*X[1,1]*X[1,2])+
(X[1,1]*X[1,1]*X[1,2]*X[2,2])+
(X[1,1]*X[1,1]*X[1,3]*X[3,2])+
(X[1,1]*X[1,2]*X[2,1]*X[1,2])+
(X[1,1]*X[1,2]*X[2,2]*X[2,2])+
(X[1,1]*X[1,2]*X[2,3]*X[3,2])+
(X[1,1]*X[1,3]*X[3,1]*X[1,2])+
(X[1,1]*X[1,3]*X[3,2]*X[2,2])+
(X[1,1]*X[1,3]*X[3,3]*X[3,2])+
(X[1,2]*X[2,1]*X[1,1]*X[1,2])+
(X[1,2]*X[2,1]*X[1,2]*X[2,2])+
(X[1,2]*X[2,1]*X[1,3]*X[3,2])+
(X[1,2]*X[2,2]*X[2,1]*X[1,2])+
(X[1,2]*X[2,2]*X[2,2]*X[2,2])+
(X[1,2]*X[2,2]*X[2,3]*X[3,2])+
(X[1,2]*X[2,3]*X[3,1]*X[1,2])+
(X[1,2]*X[2,3]*X[3,2]*X[2,2])+
(X[1,2]*X[2,3]*X[3,3]*X[3,2])+
(X[1,3]*X[3,1]*X[1,1]*X[1,2])+
(X[1,3]*X[3,1]*X[1,2]*X[2,2])+
(X[1,3]*X[3,1]*X[1,3]*X[3,2])+
(X[1,3]*X[3,2]*X[2,1]*X[1,2])+
(X[1,3]*X[3,2]*X[2,2]*X[2,2])+
(X[1,3]*X[3,2]*X[2,3]*X[3,2])+
(X[1,3]*X[3,3]*X[3,1]*X[1,2])+
(X[1,3]*X[3,3]*X[3,2]*X[2,2])+
(X[1,3]*X[3,3]*X[3,3]*X[3,2])
test1 # = 0.9288
我写了一个 for 循环来简化每个释放状态和重新捕获状态的计算 (test2)。问题是，在更多的时间间隔内，我的方法将最终计算一个具有很多维度的数组的乘积总和，我担心它在 JAGS 中的计算量会非常大。
test2 <- array((0),dim=c(3,3,3,3,3))
for (m in 1:3){
for (i in 1:3){
for (j in 1:3){
for (k in 1:3){
for (s in 1:3){
test2[m,i,j,k,s] <- X[m,i]*X[i,j]*X[j,k]*X[k,s]
}
}
}
}
}
sumtest2 <- sum(test2[1,,,,2]) # probability of state 2 at observation
# time 4 after release in state 1 at
# release time 1.
test2
sumtest2 # = 0.9288 (the same as test1 and therefore the correct result)
我的问题是实现相同目标的更有效方法是什么？
谢谢。
--------------------------------
感谢 Dex 提供的非常有用的答案。我有一个跟进问题，关于时间相关的转移矩阵。上面的代码效果很好，但我无法将它用于不同次数的迭代，其中每次迭代的转换矩阵都不同。
对于每个释放时机 t 有 T-t 重获时机（其中 T 是总学习年数）。我想用矩阵乘以匹配第一次可能夺回时间到最后一次夺回时间的矩阵。
这是我想要实现的长版的一部分：
# This long version does work, with lines for release occasion 1
# and recapture occasions 1:4 only:
pr2 <- array((0), dim=c(3,3,4,4))
pr2[,,1,1] <- diag(1, 3, 3) %*% (X[,,1])
pr2[,,1,2] <- diag(1, 3, 3) %*% (X[,,1] %*% X[,,2])
pr2[,,1,3] <- diag(1, 3, 3) %*% (X[,,1] %*% X[,,2] %*% X[,,3])
pr2[,,1,4] <- diag(1, 3, 3) %*% (X[,,1] %*% X[,,2] %*% X[,,3] %*% X[,,4])
我试过使用额外的时间索引，然后使用带有矩阵乘法的 for 循环：
# With time dependence:
X <- array((0),dim=c(3,3,4))
X[1,1,1] <- 0.1
X[1,2,1] <- 0.2
X[1,3,1] <- 0.7
X[2,1,1] <- 0.25
X[2,2,1] <- 0.35
X[2,3,1] <- 0.4
X[3,1,1] <- 0.15
X[3,2,1] <- 0.25
X[3,3,1] <- 0.6
X[,,2] <- X[,,1]+0.05
X[,,3] <- X[,,1]-0.05
X[,,4] <- X[,,1]+0.1
# not working:
library(expm)
for (t in 1:4){
for (j in t:4){
pr1[,,t,j] <- diag(1, 3, 3) %*% (X[,,j]%^%j)
}
}
# This also doesn't work:
library(expm)
for (t in 1:4){
for (j in t:4){
pr1[,,t,j] <- diag(1, 3, 3) %*% (X[,,j]%^%1:j)
}
}
我也尝试过列表方法：
# List approach
# With time dependence:
X1 <- array((0),dim=c(3,3))
X1[1,1] <- 0.1
X1[1,2] <- 0.2
X1[1,3] <- 0.7
X1[2,1] <- 0.25
X1[2,2] <- 0.35
X1[2,3] <- 0.4
X1[3,1] <- 0.15
X1[3,2] <- 0.25
X1[3,3] <- 0.6
X2 <- X1+0.05
X3 <- X1-0.05
X4 <- X1+0.1
XT <- list(X1,X2,X3,X4)
pr1 <- array((0), dim=c(3,3,4,4))
# This also doesn't work:
library(expm)
for (t in 1:4){
for (j in t:4){
pr1[,,t,j] <- diag(1, 3, 3) %*% (XT[[1:j]]%^%j)
}
}
以下代码显示了计算时间相关数组中概率的很长的方法。它给出了正确的值。完整地写出来（有 23 次释放事件和 23 次重新捕获事件，对于 3 个 release/recapture 状态），虽然可能看起来效率不高，特别是如果分析需要稍后调整。这仅显示了在第一次释放事件中首次释放的个人的前四次重新捕获事件。每个单元格给出了一个人在那个场合 第一次 被重新捕获的概率。重新捕获概率由 p 给出，因此未重新捕获的概率由 1-p 给出。在我的分析中，T 和 p 是要估计的参数。此处定义了它们的值，以确保机制按预期工作。
# Make a 3D transition matrix T where z-dimension is time
T_states <- 3
T_terms <- 4
Tz <- 4
T <- array(runif(T_states * T_states * Tz)/1.5, dim = c(T_states, T_states, Tz))
# Normalise by rows - this step ensures rows sum to 1
for (z in seq(Tz)) {
T[, , z] <- T[, , z] / apply(T[, , z], 1, sum)
}
# state and time dependent recapture probability array [,recapture state,occasion]
# the first dimension is needed to repeat p across release states
p <- array((0),dim=c(3,3,4))
p[,1,1] <- 0.5
p[,1,2] <- 0.55
p[,1,3] <- 0.6
p[,1,4] <- 0.65
p[,2,1] <- 0.6
p[,2,2] <- 0.65
p[,2,3] <- 0.7
p[,2,4] <- 0.75
p[,3,1] <- 0.7
p[,3,2] <- 0.75
p[,3,3] <- 0.8
p[,3,4] <- 0.85
# Put them together (with 1-p for earlier recapture occasions):
pr2 <- array((0), dim=c(3,3,4,4))
pr2[,,1,1] <- diag(1, 3, 3) %*% (T[,,1]*p[,,1])
pr2[,,1,2] <- diag(1, 3, 3) %*% ((T[,,1]*(1-p[,,1])) %*% (T[,,2]*p[,,2]))
pr2[,,1,3] <- diag(1, 3, 3) %*% ((T[,,1]*(1-p[,,1])) %*% (T[,,2]*(1-p[,,2])) %*% (T[,,3]*p[,,1]))
pr2[,,1,4] <- diag(1, 3, 3) %*% ((T[,,1]*(1-p[,,1])) %*% (T[,,2]*(1-p[,,2])) %*% (T[,,3]*(1-p[,,3])) %*% (T[,,4]*p[,,4]))
pr2[,,1,]
我一直没能找到一种方法来为此编写更高效的代码来提供相同的值（很多方法都提供不正确的值）。
你可以用矩阵乘法来解决这个问题。转移矩阵可以矩阵相乘在一起，以在多次迭代中给出一个转移矩阵。对于四次迭代，我们执行四次：
(X %*% X %*% X %*% X)
您还可以使用 expm 包中的 %^% 运算符来进行矩阵求幂，而无需将其写出。
要找出对象在四次迭代后从状态 1 开始到状态 2 结束的可能性，我们可以将其表示为 1x3 矩阵并将其乘以转移矩阵四次。
matrix(c(1, 0, 0), nrow = 1, ncol = 3) %*% (X %*% X %*% X %*% X)
这给出了它处于三种状态的概率：
[,1] [,2] [,3]
[1,] 0.756 0.9288 1.1016
时机
我根据您的代码制作了另一个示例，其中的转换矩阵具有合理的概率。矩阵乘法比 for 循环快大约 100 倍并且给出相同的答案。
library("microbenchmark")
X <- array((0),dim=c(3,3))
X[1,1] <- 0.1
X[1,2] <- 0.2
X[1,3] <- 0.7
X[2,1] <- 0.25
X[2,2] <- 0.35
X[2,3] <- 0.4
X[3,1] <- 0.15
X[3,2] <- 0.25
X[3,3] <- 0.6
microbenchmark(
{test2 <- array((0),dim=c(3,3,3,3,3))
for (m in 1:3){
for (i in 1:3){
for (j in 1:3){
for (k in 1:3){
for (s in 1:3){
test2[m,i,j,k,s] <- X[m,i]*X[i,j]*X[j,k]*X[k,s]
}
}
}
}
}
sumtest2 <- sum(test2[1,,,,2])},
matrix(c(1, 0, 0), nrow = 1, ncol = 3) %*% X %*% X %*% X %*% X
)