我正在研究易感-感染-恢复模型,该模型定义在元群体上,并通过基于距离的评分将细胞成对联系起来。在这样做的时候,我的状态会降到零以下,(S)。我不明白为什么。
事实证明,这可能发生,因为不知何故,没有“约束”的变化。这是我还是有解决方案?
我试图创造一个更小更直接的例子。这一个使用随机值,但它只是为了表明,当它相对较小时,速率是好的,但如果它发生变化,由于许多附近的细胞被感染(这里是由于随机数),那么估计失败。
library(deSolve)
ode(
y = c(
S1 = 10,
I1 = 0.001),
times = 0:5000,
parms = list(# beta = 0.05
beta = 0.01),
func = function(times, state, parms) {
with(parms, {
# S <- state[c(1, 3)]
# I <- state[c(2, 4)]
S <- state[1]
I <- state[2]
beta <- beta + 2 * rbinom(1, size = 1, prob = 0.4)
list(c(dS = -beta * S * I,
dI = beta * S * I))
})
},
jacfunc = function(t, y, p) {
with(p, {
S <- y[1]
I <- y[2]
(matrix(
c(beta * I, beta * S, -beta * I, -beta * S),
nrow = 2,
ncol = 2,
byrow = TRUE
))
})
},
jactype = "fullusr",
verbose = TRUE
) %>%
# diagnostics() %>%
print() %>%
# zapsmall()
identity()
- 指定
hmin
并不能挽救这一点 - 我也试过提供一个雅可比矩阵,但这也无济于事。
1条答案
按热度按时间pbossiut1#
诊断
该守则包含两个问题:
1.模型函数需要是严格确定的,并且不能包含随机分量,因为求解器不能收敛。在ode系统中,时间根据定义是连续的,因此随机分量在内部是没有意义的。从技术上讲,求解器将为每个迭代步骤绘制不同的随机值,因此模型本质上是未定义的。如果意图是提供可变参数
beta
,则必须预先用明确定义的时间步长来定义。1.代码行
beta + 2 * rbinom(1, size = 1, prob = 0.4)
导致beta
的值在0.01和2.01之间切换。对于这种类型的模型,2.01的比率是不切实际的大,而不仅仅是“适度”。从技术上讲,它贯穿并给出了一个合理的数字,S
立即下降到零。然而,由于求解器具有有限的精度,因此它可以达到(小的)负值。为了得到一个合理的数字,可以将
beta
直接乘以beta_t
,或者只将随机部分乘以<< 2的值。解决方案
将随机分量实现为强制
将随机分量存储在输入向量(例如
beta_r
),其长度与times
相同。模型函数中的time
参数始终包含模拟的当前时间。其旨在访问时变输入。因为R中的向量以索引1开始,所以使用floor(time) + 1
。输入向量(通常称为强制)可以作为全局变量或作为附加的命名参数传递给func
。为了简化示例,我省略了雅可比矩阵。
一系列模拟
要运行一系列模拟,可以使用循环或
lapply
。deSolve::plot
-函数允许直接绘制这样的列表。流水线方法也是可能的。针对大梯度调整求解器参数
如前所述,
beta = 2.0
对于这样的系统来说非常大。但是,如果真的打算这样做,可以减少时间步长,外部使用times
或内部使用hmax
。当然,它会变得非常缓慢。