0

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

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)情节。因此，当我使用ggplotand时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')
4

1 回答 1

2

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()withtmpf将结果添加到绘图中的原因。

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

在此处输入图像描述

于 2015-12-08T17:16:17.087 回答