【问题标题】:Exponential regression with nls in R [duplicate]R中nls的指数回归[重复]
【发布时间】:2014-10-25 09:07:16
【问题描述】:

我正在尝试解决以下问题:

x <- c(0.11557577149788574,2.1552479877306925,2.5505873377321175,1.0995198836006757,3.710225290286669,2.386870541964232,0.11557577149788574,0.11557577149788574,2.1552479877306925,2.5505873377321175,1.0995198836006757,3.710225290286669,2.386870541964232,0.11557577149788574)
y <- c(16500,11500,11500,13630,7000,11995,13490,16500,11500,11500,13630,7000,11995,13490)

df <- data.frame(x, y)

m <- nls(y ~ I(a*exp(-b*x)+c), data=df, start=list(a=14000, b=1, c=100), control=nls.control(maxiter=10000, minFactor=1e-7, tol=1e-5, printEval=F, warnOnly=F))

但是,即使我尝试更改起始值并且 nls 控件也不会返回任何值。 我做错了什么?我需要更多积分来解决这个问题吗?

谢谢!

【问题讨论】:

    标签: r regression exponential


    【解决方案1】:

    您需要更好的起始值。

    首先以递增的x 对数据框进行排序,以便绘制得很好。

    如果我们将c 设置为零,我们可以将更简单的模型拟合到log(y) ~ A - b*x,它在所有系数中都是线性的,因此可以通过lm 进行拟合,并且不需要起始值。

    使用该更简单模型给出的b 的起始值。另请注意,a 和 c 线性输入完整模型,因此我们可以使用 nls 的 plinear 算法,从而无需为这些系数提供起始值:

    o <- order(df$x)
    df_o <- df[o, ] # order it by increasing x
    
    fm0 <- lm(log(y) ~ I(-x), df_o) # simpler model to get better starting values
    
    st <- list(b = coef(fm0)[[2]])
    fm <- nls(y ~ cbind(exp(-b*x), 1), df_o, start = st, alg = "plinear")
    
    plot(df_o, col = "red", pch = 20)
    lines(fitted(fm) ~ x, df_o)
    

    结果如下,其中.lin1 是a,.lin2 是c:

    > fm
    Nonlinear regression model
      model: y ~ cbind(exp(-b * x), 1)
       data: df_o
             b      .lin1      .lin2 
       -0.4903 -1529.0253 16509.4422 
     residual sum-of-squares: 10555038
    
    Number of iterations to convergence: 5 
    Achieved convergence tolerance: 8.504e-07
    

    这里我们将输入数据显示为红点,并通过拟合结果画线:

    【讨论】:

    • 真的很好的答案,我想知道你为什么放I(-x),是因为数据吗?例如,如果我有一个从 0 开始的值会发生什么,我们会遇到 log(0) 的问题,有没有办法解决它?谢谢!
    • 在 lm 公式中,-x 表示删除变量 x,但 I(-x) 表示负 x。关于 log(0),lm 有一个子集参数,可用于删除 y=0 行。
    • 但是去掉 0 值会影响拟合吗?
    • 日志仅用于获取起始值。
    • 使用 df_o
    【解决方案2】:

    问题似乎出在您的数据中。

    如果您尝试plot( y ~ x ),那么它看起来并不是指数级的:实际上它似乎比线性略快,特别是如果您发现您的数据中实际上有两个点具有x y 值3.710225 7000 所以您正在尝试将凸曲线拟合到表明凹曲线的数据。两个建议:

    1. 去掉这两点:

      df2 &lt;- df[df$y!=7000,]

      m2 &lt;- nls(y ~ I(a*exp(-b*x)+c), data=df2, start=list(a=14000, b=1, c=100), control=nls.control(maxiter=10000, minFactor=1e-7, tol=1e-5, printEval=F, warnOnly=F))

    给予

    Nonlinear regression model
      model: y ~ I(a * exp(-b * x) + c)
       data: df2
            a         b         c 
    1.418e+04 1.202e-01 1.030e+03 
     residual sum-of-squares: 9781328
    
    Number of iterations to convergence: 14 
    Achieved convergence tolerance: 2.573e-06
    
    1. 通过颠倒a 和b 的符号,让您的起点建议一个凹曲线

      m &lt;- nls(y ~ I(a*exp(-b*x)+c), data=df, start=list(a=-14000, b=-1, c=100), control=nls.control(maxiter=10000, minFactor=1e-7, tol=1e-5, printEval=F, warnOnly=F))

    给予

    Nonlinear regression model
      model: y ~ I(a * exp(-b * x) + c)
       data: df
             a          b          c 
    -1529.0204    -0.4903 16509.4360 
     residual sum-of-squares: 10555038
    
    Number of iterations to convergence: 7 
    Achieved convergence tolerance: 7.471e-07
    

    【讨论】:

    • “y ~ I(a * exp(-b * x) + c)”中的I代表什么?
    • @William 我从原始问题中得到它。在 R 中,它表示AsIs,在函数formula 中用于禁止将“+”、“-”、“*”和“^”等运算符解释为公式运算符,因此它们用作算术运算符
    • 感谢您的解释! :D
    【解决方案3】:

    您的数据集中没有足够的信息来可靠地估计此非线性模型的系数。您的数据中几乎没有曲率。如果您使用 Levenberg-Marquardt,它会遇到或多或少等于线性拟合的参数值:

    plot(y~x)
    
    library(minpack.lm)
    m <- nlsLM(y ~ a*exp(-b*x)+c, data=df, start=list(a=14000, b=1, c=100),  
               control = nls.lm.control(maxiter = 1000))
    summary(m)
    #Formula: y ~ a * exp(-b * x) + c
    #
    #Parameters:
    #    Estimate Std. Error t value Pr(>|t|)
    #a  2.308e+06  6.822e+08   0.003    0.997
    #b  8.506e-04  2.518e-01   0.003    0.997
    #c -2.292e+06  6.822e+08  -0.003    0.997
    #
    #Residual standard error: 1250 on 11 degrees of freedom
    
    #Number of iterations till stop: 94 
    #Achieved convergence tolerance: 1.49e-08
    #Reason stopped: Number of calls to `fcn' has reached or exceeded `maxfev' == 400.
    
    lines(x, predict(m), col="green", lwd=2)
    abline(lm(y~x), col="red", lty=2, lwd=2)
    

    如果模型是基于科学的,那么您测量的 x 值的范围可能太小了。

    【讨论】:

      猜你喜欢
      • 2021-07-13
      • 2012-07-24
      • 2014-12-16
      • 2020-09-19
      • 2016-02-23
      • 2020-05-02
      • 2014-11-28
      • 1970-01-01
      • 2016-09-14
      相关资源
      最近更新 更多