使用 nls 对具有积分功能的函数进行参数估计

Question

我正在尝试使用 nls 对使用 integrate 的函数进行参数估计。涉及 integrate 的函数 intsol1 似乎可以自行计算其参数的各个值。

但是，当我尝试在nls中调用它时，集成不起作用。

我得到的错误是： Integrate(aint1, lower = 0, upper = t) 中的错误： 长度（上）== 1 不是 TRUE

当我尝试在 ggplot 中使用 geom_function 绘制intsol1时，我遇到了同样的问题。

任何对此的意见都将受到赞赏。

library(nls.multstart)
library(aomisc)

diffusion1d <- function(x=1, t, diff=1) {
  sig <- sqrt(t*diff*2)
  df <- dnorm(x, mean=0, sd=sig)
  df <- ifelse(t == 0, 0, df)
}

intsol1 <- function(x=1, t, diff=1, cl=1) {
  x <- x
  diff <- diff
  cl <- cl
  # Integrating factor solution 1-D diffusion with reaction 
  aint1 <- function(t) {exp(cl*t)*diffusion1d(x, t, diff)}
  aint1 <- Vectorize(aint1)
  bint1 <- exp(-cl*t)*integrate(aint1, lower=0, upper=t)$value
  return(bint1)
}

nlsdata <- data.frame(Time_days=c(0.4032258, 1.6129032, 2.0161290, 2.4193548, 
                                    3.6290323, 4.8387097, 6.4516129, 10.0806452, 
                                    13.7096774, 21.3709677, 28.2258065, 
                                    35.0806452, 41.9354839, 56.0483871, 
                                    69.7580645, 83.4677419), 
                      Concentration_mcg_ml=c(1.5151515, 3.3333333, 4.5454545,
                                               5.2727273, 5.5757576, 5.8181818,
                                               5.5757576, 4.4848485, 3.8787879, 
                                               2.6060606, 1.7575758, 1.2727273, 
                                               0.9090909, 0.4242424, 0.2424242, 
                                               0.1212121))

# x, diff, cl are parameters to be estimated with nls

nlsfit1 <- nls_multstart(Concentration_mcg_ml ~ 
                           intsol1(x, Time_days, diff, cl), data=nlsdata, 
                         lower=c(x=0, diff=0, cl=0), 
                         upper=c(x=1, diff=1, cl=1), 
                         start_lower=c(x=0, diff=0, cl=0), 
                         start_upper=c(x=1, diff=1, cl=1), 
                         iter=10, 
                         supp_errors="N")

nlssummary1 <- summary(nlsfit1)

Answer 1

考虑循环

t

。请注意，您的代码可以简化为以下代码：

intsol_vec <- function(t, x = 1, diff = 1, cl = 1) {
  aint1 <- function(t) exp(cl*t)*dnorm(x, 0, sqrt(2*t*diff))
  sapply(t, \(t)exp(-cl*t)*integrate(aint1, lower=0, upper=t)$value)
}

intsol_vec(1:5)
[1] 0.1335144 0.1707453 0.1638357 0.1482900 0.1335234

---

将结果与您的函数进行比较：

diffusion1d <- function(x=1, t, diff=1) {
     sig <- sqrt(t*diff*2)
     df <- dnorm(x, mean=0, sd=sig)
     df <- ifelse(t == 0, 0, df)
 }
 
intsol1 <- function(x=1, t, diff=1, cl=1) {
     x <- x
     diff <- diff
     cl <- cl
     # Integrating factor solution 1-D diffusion with reaction 
     aint1 <- function(t) {exp(cl*t)*diffusion1d(x, t, diff)}
     aint1 <- Vectorize(aint1)
     bint1 <- exp(-cl*t)*integrate(aint1, lower=0, upper=t)$value
     return(bint1)
 }
mapply(intsol1, t=1:5)
[1] 0.1335144 0.1707453 0.1638357 0.1482900 0.1335234