拟合微分方程：如何将一组数据拟合到 R 中的微分方程答案

【问题标题】：Fitting differential equations: how to fit a set of data to a differential equation in R拟合微分方程：如何将一组数据拟合到 R 中的微分方程
【发布时间】：2015-02-02 17:27:50
【问题描述】：

有一个数据集：

conc <- data.frame(time = c(0.16, 0.5, 1.0, 1.5, 2.0, 2.5, 3), concentration = c(170, 122, 74, 45, 28, 17, 10))

我想将这些数据拟合到下面的微分方程中：

dC/dt= -kC

其中 C 是数据集中的浓度和时间 t。这也将给出 k 的结果。谁能给我一个线索如何在 R 中做到这一点？谢谢。

【问题讨论】：

请在data.frame(.) 中使用=，不要 <-。
这似乎是家庭作业。
其实不是功课。数据是药代动力学中最简化的一室模型。我举一个简单例子的原因是我相信类似的方法可以用来解决多个隔室问题，因此在 PK-PD 模拟中可以使用一些链接模型。

标签： r differential-equations

【解决方案1】：

首先使用变量分离来求解微分方程。这给出了 log(C)=-k*t+C0。

绘制数据：

plot(log(concentration) ~ time,data=conc)

拟合线性模型：

fit <- lm(log(concentration) ~ time,data=conc)
summary(fit)

# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   5.299355   0.009787   541.4 4.08e-13 ***
#   time       -0.992208   0.005426  -182.9 9.28e-11 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# 
# Residual standard error: 0.01388 on 5 degrees of freedom
# Multiple R-squared: 0.9999,  Adjusted R-squared: 0.9998 
# F-statistic: 3.344e+04 on 1 and 5 DF,  p-value: 9.281e-11

绘制预测值：

lines(predict(fit)~conc$time)

提取k：

k <- -coef(fit)[2]
#0.9922081

【讨论】：

【解决方案2】：

这可能是一个解决方案：

    require('deSolve')
conc <- data.frame(time <- c(0.16, 0.5, 1.0, 1.5, 2.0, 2.5, 3), concentration <- c(170, 122, 74, 45, 28, 17, 10))

##"Model" with differential equation
model <- function(t, C, k){
  list(-k * C)
}

##Cost function with sum of squared residuals:
cost <- function(k){
  c.start <- 170
  out <- lsoda(c(C=c.start), conc$time, model, k)
  c <- sum( (conc$concentration - out[,"C"])^2)       
  c
}

##Initial value for k
k <- 3
##  Use some optimization procedure
opt <- optim(k, cost, method="Brent", lower=0, upper=10)

k.fitted <- opt$par

也许这有点天真，因为使用 lsoda 对于仅使用一个微分方程进行计算似乎有点过头了......但它肯定会优化你的 k。您可能想检查集成的 C 的起始值，我在这里将其设置为 170，您没有 t=0 的值吗？

【讨论】：

非常感谢。我刚刚发现 FME 中的 modCost 函数也可能是一个解决方案，尽管不像您提供的答案那么直接。 inside-r.org/packages/cran/FME/docs/modCost
您对如何为多个参数执行此操作有任何参考吗？或者，你能在这个问题上花几分钟时间告诉我我哪里错了吗？谢谢stackoverflow.com/questions/15730793/…
成本函数内的c <- sum((conc$concentration....c.170..122..74..45..28..17..10. - out[,"C"])^2) 行是错误的。我将其更正为：c <- sum( (conc$concentration - out[,"C"])^2) ，但我不知道它何时会经过同行评审并因此可见。