r - 尝试用 R 和 nls 在一个有条件的函数上拟合数据

Question

我正在尝试将一些数据拟合到具有有效性限制的函数中。更准确地说，如果 t<=T 和 t>T，则具有不同值的函数。

这是我试过的代码:

posExpDecay <- function(t,tau,max,toff){ 1+max*(1-exp(-(t-toff)/tau)) } 
negExpDecay <- function(t,tau,max){ 1+max*exp(-(t)/tau) } 

data<-structure(list(t = c(0.67, 1, 1.33, 1.67, 2, 4, 6, 8, 10), y = c(1.02,2.33, 3.08, 3.34, 3.41,2.50, 1.86, 1.44, 1.22)), .Names = c("t", "y"), row.names = c(13L, 17L, 21L, 25L, 29L,37L, 45L, 49L, 53L), class = "data.frame")

fit <- nls(y~ifelse(t<=tswitch,
                    posExpDecay(t,tau1,max1,toff),
                    negExpDecay(t,tau2,max2)),
                  data,
                  start=list(max1=3,tau1=0.7,max2=7,tau2=2,toff=0.1,tswitch=3))

我收到以下错误:

Error in nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter estimates

这是我的起始参数不够好(我尝试了几个)，我的问题在 R 中没有很好地翻译，还是我错过了一个基本的数学错误？

Answer 1

nls(...) 默认使用高斯牛顿法；该错误消息实际上很常见，表示雅可比矩阵无法求逆。

我认为您的问题与以下事实有关:对于其他参数的任意值，您的复合函数(公式的 RHS)在 t=tswitch 处不连续。换句话说，函数连续的要求对其他参数施加了约束——它们不是彼此独立的。此外，复合函数的导数在 t=tswitch 处永远不会连续 - 您的 posExpDecay(...) 对所有 t 都有正导数code>，而你的 negExpDecay(...) 对所有 t 都有负导数。

我不知道这种函数形式是否有理论上的原因，但这些 +/- 指数通常使用正衰减和负衰减的乘积建模，如下所示。

注意:我通常使用 minpack.lm 包中的 nlsLM(...)，它使用更强大的 Levenberg Marquardt 算法。它与基础 R 中的 nls(...) 函数具有相同的签名。

f <- function(t, max,tau1,tau2,toff) max*exp(-t/tau1)*(1-exp(-(t-toff)/tau2))
library(minpack.lm)
fit <- nlsLM(y~f(t,max,tau1,tau2,toff),data,
             start=list(max=15,tau1=0.7,tau2=2,toff=.2))
summary(fit)
# ...
# Parameters:
#      Estimate Std. Error t value Pr(>|t|)    
# max   4.72907    0.29722  15.911 1.78e-05 ***
# tau1  6.75926    0.54093  12.496 5.82e-05 ***
# tau2  0.51211    0.08209   6.238  0.00155 ** 
# toff  0.53595    0.02667  20.093 5.64e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.113 on 5 degrees of freedom
# 
# Number of iterations to convergence: 19 
# Achieved convergence tolerance: 1.49e-08

plot(y~t,data)
curve(predict(fit,data.frame(t=x)),add=T,col="blue")

如您所见，这个更简单的函数(更少的参数)非常适合。