r - 用 nlme 和 lsoda 拟合一阶方程
问题描述
nlme
我试图使用和来拟合一阶微分模型lsoda
。这是基本思想:我首先定义允许生成微分方程解的函数:
library(deSolve)
ODE1 <- function(time, x, parms) {with(as.list(c(parms, x)), {
import <- excfunc(time)
dS <- import*k/tau - (S-yo)/tau
res <- c(dS)
list(res)})}
solution_ODE1 = function(tau1,k1,yo1,excitation,time){
excfunc <- approxfun(time, excitation, rule = 2)
parms <- c(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
xstart = c(S = yo1)
out <- lsoda(xstart, time, ODE1, parms)
return(out[,2])
}
然后,我根据两个 ID 的等式生成数据:
time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(signal = c(solution_ODE1(3,2,0.1,excitation,time)+rnorm(length(time),0,0.1),
solution_ODE1(3.2,1.5,0.3,excitation,time)+rnorm(length(time),0,0.1)),
time = rep(time,2),
excitation = rep(excitation,2),
ID = rep(c("A","B"),each = length(time)))
这是它的样子:
library(ggplot2)
ggplot(simu_data)+
geom_point(aes(time,signal,color = "signal"),size = 2)+
geom_line(aes(time,excitation,color = "excitation"))+
facet_wrap(~ID)
然后我尝试使用 nlme 来适应:
fit1 <- nlme(signal ~ solution_ODE1(damping,gain,eq,excitation,time),
data = simu_data,
fixed = damping + gain + eq ~1,
random = damping ~ 1 ,
groups = ~ ID,
start = c(damping = 5, gain = 1,eq = 0))
我得到了这个错误,我没有得到:
eval(替代(expr),数据,enclos = parent.frame())中的错误:找不到对象'k'
这traceback
表明错误来自 ODE1 模型,该模型在生成值时起作用。
16. eval(substitute(expr), data, enclos = parent.frame())
15. eval(substitute(expr), data, enclos = parent.frame())
14. with.default(as.list(c(parms, x)), {
import <- excfunc(time)
dS <- import * k/tau - (S - yo)/tau
res <- c(dS) ...
13. with(as.list(c(parms, x)), {
import <- excfunc(time)
dS <- import * k/tau - (S - yo)/tau
res <- c(dS) ...
12. func(time, state, parms, ...)
11. Func2(times[1], y)
10. eval(Func2(times[1], y), rho)
9. checkFunc(Func2, times, y, rho)
8. lsoda(xstart, time, ODE1, parms)
7. solution_ODE1(damping, gain, eq, excitation, time)
6. eval(model, data.frame(data, pars))
5. eval(model, data.frame(data, pars))
4. eval(modelExpression[[2]], envir = nlEnv)
3. eval(modelExpression[[2]], envir = nlEnv)
2. nlme.formula(signal ~ solution_ODE1(damping, gain, eq, excitation,
time), data = simu_data, fixed = damping + gain + eq ~ 1,
random = damping ~ 1, groups = ~ID, start = c(damping = 5,
gain = 1, eq = 0))
1. nlme(signal ~ solution_ODE1(damping, gain, eq, excitation, time),
data = simu_data, fixed = damping + gain + eq ~ 1, random = damping ~
1, groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))
有谁知道我应该如何进行?
编辑
我尝试按照 mikeck 的建议进行修改:
ODE1 <- function(time, x, parms) {
import <- parms$excfunc(time)
dS <- import*parms$k/parms$tau - (x["S"]-parms$yo)/parms$tau
res <- c(dS)
list(res)}
生成数据没有问题。但是nlme
现在使用给定:
checkFunc(Func2, times, y, rho) 中的错误:func() (0) 返回的导数数量必须等于初始条件向量的长度 (100)
具有以下回溯:
> traceback()
10: stop(paste("The number of derivatives returned by func() (",
length(tmp[[1]]), ") must equal the length of the initial conditions vector (",
length(y), ")", sep = ""))
9: checkFunc(Func2, times, y, rho)
8: lsoda(xstart, time, ODE1, parms) at #5
7: solution_ODE1(damping, gain, eq, excitation, time)
6: eval(model, data.frame(data, pars))
5: eval(model, data.frame(data, pars))
4: eval(modelExpression[[2]], envir = nlEnv)
3: eval(modelExpression[[2]], envir = nlEnv)
2: nlme.formula(signal ~ solution_ODE1(damping, gain, eq, excitation,
time), data = simu_data, fixed = damping + gain + eq ~ 1,
random = damping ~ 1, groups = ~ID, start = c(damping = 5,
gain = 1, eq = 0))
1: nlme(signal ~ solution_ODE1(damping, gain, eq, excitation, time),
data = simu_data, fixed = damping + gain + eq ~ 1, random = damping ~
1, groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))
解决方案
在您的示例中,您的times
向量不会单调运行。我认为这与lsoda
. 时间在这里运作的方式的背景/意义是什么?用两组拟合随机效应模型并没有什么意义。您是否试图将相同的曲线拟合到两个独立的时间序列?
这是一个精简的示例,进行了一些调整(并非所有内容都可以折叠为数字向量而不会丢失必要的结构):
library(deSolve)
ODE1 <- function(time, x, parms) {
with(as.list(parms), {
import <- excfunc(time)
dS <- import*k/tau - (x-yo)/tau
res <- c(dS)
list(res)
})
}
solution_ODE1 = function(tau1,k1,yo1,excitation,time){
excfunc <- approxfun(time, excitation, rule = 2)
parms <- list(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
xstart = yo1
out <- lsoda(xstart, time, ODE1, parms)
return(out[,2])
}
time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(time = rep(time,2),
excitation = rep(excitation,2))
svec <- c(damping = 3, gain = 1.75, eq = 0.2)
这有效:
with(c(simu_data, as.list(svec)),
solution_ODE1(damping,gain,eq,excitation[1:50],time[1:50]))
但是,如果我们再包含一个步骤(以便时间重置为 0),它将失败:
with(c(simu_data, as.list(svec)),
solution_ODE1(damping,gain,eq,excitation[1:51],time[1:51]))
lsoda(xstart, time, ODE1, parms) 中的错误:在采取任何集成步骤之前检测到非法输入 - 请参阅书面消息
推荐阅读
- angular - 具有左右间距和固定页脚的 Angular Flex 布局列
- laravel - laravel - 混合清单不存在
- java - 如何从另一个位置添加 textView 的第二行
- python - 使用 pandas 计算一天中的事件持续时间
- python - 绘制象限图以根据 X 和 Y 的平均值区分 4 组中的人口并找到最终计数
- python - 试图通过 selenium 抓取但不能遍历动态内容
- spring-boot - Spring Security:使用域逻辑进行授权
- bash - Cannot write pipe character | in IntelliJ Terminal
- flutter - 当我使用带有标头授权令牌的颤振 http 库时,它在控制台中显示 500。但在邮递员中一切都好
- amazon-web-services - Terraform + GithubActions