使用nls对含integrate的函数做参数估计时遇报错求助

问题描述 #

调用intsol1传入单个参数值时计算正常，但在nls_multstart或ggplot的geom_function中调用时，触发如下报错：

Error in integrate(aint1, lower = 0, upper = t) :
length(upper) == 1 is not TRUE

问题原因 #

nls系列函数和geom_function会向量化传入参数（即t会以向量形式传入），但integrate的upper参数仅接受单个数值，无法直接处理向量输入。虽然你对内部的aint1做了向量化，但外层的积分逻辑仍未适配向量t的场景。

解决方案 #

对整个intsol1函数做向量化处理，让它能逐个处理向量中的每个t值。以下提供两种实现方式：

修改后的代码 #

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)
}

# 方案1：用Vectorize包装整个intsol1函数（简洁版）
intsol1 <- Vectorize(function(x=1, t, diff=1, cl=1) {
  # 积分变量改用t_val避免和外部t冲突
  aint1 <- function(t_val) {exp(cl*t_val)*diffusion1d(x, t_val, diff)}
  bint1 <- exp(-cl*t)*integrate(aint1, lower=0, upper=t)$value
  return(bint1)
}, vectorize.args = "t")

# 方案2：在函数内部用sapply遍历t向量（直观可控版）
# intsol1 <- function(x=1, t, diff=1, cl=1) {
#   sapply(t, function(t_val) {
#     aint1 <- function(t_int) {exp(cl*t_int)*diffusion1d(x, t_int, diff)}
#     exp(-cl*t_val)*integrate(aint1, lower=0, upper=t_val)$value
#   })
# }

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))

# 重新执行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)
print(nlssummary1)

关键说明 #

• 两种方案均能让函数同时适配单个t值和向量t的输入场景，完美兼容nls和ggplot的调用需求。
• 修改时注意将积分内部的变量名与外部的t区分开（如改用t_val/t_int），避免变量作用域冲突。

内容的提问来源于stack exchange，提问作者Murali