以 exp(-exp(-x/c)) 为模型的非线性回归 (R,nls)答案

【问题标题】：nonlinear regression with exp(-exp(-x/c)) as model (R,nls)以 exp(-exp(-x/c)) 为模型的非线性回归 (R,nls)
【发布时间】：2015-12-07 13:57:51
【问题描述】：

我很难将以下模型拟合到后续数据：

y = a * exp( - b * exp( - x/c))

DF$x=c(200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000,4100,4200,4300,4400,4500,4600,4700,4800,4900,5000,5100,5200,5300,5400,5500,5600,5700,5800,5900,6000,6100,6200,6300,6400,6500,6600,6700,6800,6900,7000)
DF$y=c(50,150,350,550,1050,1650,2950,4750,7850,12350,18950,27250,36750,49750,63250,79050,95450,112850,134550,158050,184650,211150,237750,270150,299650,334850,373450,413050,453050,490350,534250,574050,622750,666550,707350,760250,803050,848650,893850,928250,973850,1006250,1047650,1075850,1113850,1146350,1180150,1212650,1243950,1275850,1306250,1332150,1372350,1402550,1440650,1471550,1503550,1549850,1583150,1628850,1664250,1711250,1746850,1793250,1837950,1884750,1930850,1976750,2008650)

实际上，我需要一个log(y, base = 10) 情节。因此，当我使用ggplot 和stat_smooth(method='nls',...) 时，我会得到一个可怕的配合。而且我不知道如何提高合身性。你能想出一个合理的方法吗？谢谢！

到目前为止我的代码：

ggplot(DF, aes(x = x, y = y)) + geom_point() + 
stat_smooth(method = 'nls', 
    formula = y ~ a * exp(- b * exp( - x / 50 )), 
    aes(colour = 'Exponential'), 
    se = FALSE,   start = list(a=1000,b=10)) +
    scale_y_continuous(trans='log10')

【问题讨论】：

您的代码是正确的。您需要更好的起始值或更好的模型。
是的。我已经得到了合理的起始值（参见添加的代码）。不过，我需要改善合身的方法。我对 R 拟合相当陌生，所以我需要一些深入的方法和这个明确的例子。该模型基于实验数据。此处显示的数据来自模拟。
` y ~ SSasymp(x, Asym, R0, lrc)` 在没有起始值的情况下效果很好
谢谢，我会调查的。不过，我发布的模型是必不可少的。我无法替换它。我的数据与实验数据非常相似，所以应该有一个令人满意的拟合。 so else 已经发布了对实验数据的拟合。
您能否澄清这是否是一个 ggplot 问题，或者让 nls 作为独立操作（然后可能将其添加到 ggplot 调用）是否就足够了？跨度>

标签： r ggplot2 non-linear-regression

【解决方案1】：

tl;dr 您正在尝试拟合 Gompertz 曲线，而 R 有一个 SSgompertz 函数可以解决问题。

数据：

x <- seq(200,7000,by=100)
y <- c(50,150,350,550,1050,1650,2950,4750,7850,12350,18950,27250,36750,49750,63250,79050,95450,112850,134550,158050,184650,211150,237750,270150,299650,334850,373450,413050,453050,490350,534250,574050,622750,666550,707350,760250,803050,848650,893850,928250,973850,1006250,1047650,1075850,1113850,1146350,1180150,1212650,1243950,1275850,1306250,1332150,1372350,1402550,1440650,1471550,1503550,1549850,1583150,1628850,1664250,1711250,1746850,1793250,1837950,1884750,1930850,1976750,2008650)
DF <- data.frame(x,y)

R的自启动Gompertz函数（?SSgompertz）的参数化为 a*exp(-b2*b3^x) = a*exp(-b2*exp(x*log(b3)) = a*exp(-b2*exp(x/(1/log(b3)))

library("ggplot2"); theme_set(theme_bw())
n1 <- nls(y ~ SSgompertz(x, p1,p2,p3 ), data=DF)
coef(n1)
##           p1           p2          p3
## 2.530746e+06 6.907014e+00 9.995405e-01

将其翻译成您的术语：

coef2 <- with(as.list(coef(n1)),
            c(a=p1,b=p2,c=-1/log(p3)))

##             a            b             c 
##  2.530746e+06 6.907014e+00  2.175813e+03

检查：

tmpf0 <- function(x) {
    with(as.list(coef(n1)),p1*exp(-p2*p3^x))
}

tmpf <- function(x) {
    y <- with(as.list(coef2),a*exp(-b*exp(-x/c)))
    print(range(y))
    return(y)
}
tmpf0(200)  ## 4645.226
tmpf(200)   ## ditto

我可以很容易地用预测值绘制绘图，但由于我目前无法弄清楚的原因，我不能 (1) 将拟合嵌入 ggplot2 (2) 使用 stat_function() 和 tmpf 到将结果添加到绘图中。

DF$pred <- predict(n1)
g0 <- ggplot(DF,aes(x,y))+geom_point()
g0+ geom_line(aes(y=pred),colour="red") +
        scale_y_log10()

【讨论】：