R语言使用nls拟合,为什么总说循环次数大于50

Python019

R语言使用nls拟合,为什么总说循环次数大于50,第1张

nls的数据源必须有误差。不能精确等于公式返回值(零残差)。循环次数大于50通常是使用 函数精确返回值 作为数据源去拟合函数。必须给y值加上随机误差。

z=function(x,a,b){a*sin(x)+b*cos(x)}

x=seq(1,10,9/500)

y=z(x,1,1) # a=1 b=1 是期望拟合出的结果。

cor=data.frame(x=x,y=y)

cor$res=runif(length(cor$x),min=-0.005,max=0.005)

cor$yres=cor$y+cor$res

#yres =y加上随机误差,y是精确返回值

> nls(cor$yres~z(cor$x,a,b),data=cor,start=list(a=0.8,b=1.3))

Nonlinear regression model

  model: cor$yres ~ z(cor$x, a, b)

   data: cor

     a      b 

0.9999 1.0002 

 residual sum-of-squares: 0.004213

Number of iterations to convergence: 1 

Achieved convergence tolerance: 2.554e-07

#使用精确返回值拟合就会出错。

> nls(cor$y~z(cor$x,a,b),data=cor,start=list(a=1,b=1))

Error in nls(cor$y ~ z(cor$x, a, b), data = cor, start = list(a = 1, b = 1)) : 

  循环次数超过了50这个最大值

模型拟合

对于人口模型可以采用Logistic增长函数形式,它考虑了初期的指数增长以及总资源的限制。其函数形式如下。

首先载入car包以便读取数据,然后使用nls函数进行建模,其中theta1、theta2、theta3表示三个待估计参数,start设置了参数初始值,设定trace为真以显示迭代过程。nls函数默认采用Gauss-Newton方法寻找极值,迭代过程中第一列为RSS值,后面三列是各参数估计值。然后用summary返回回归结果。

library(car)

pop.mod1 <- nls(population ~ theta1/(1+exp(-(theta2+theta3*year))),start=list(theta1 = 400, theta2 = -49, theta3 = 0.025), data=USPop, trace=T)

summary(pop.mod)

在上面的回归过程中我们直接指定参数初始值,另一种方法是采用搜索策略,首先确定参数取值范围,然后利用nls2包的暴力方法来得到最优参数。但这种方法相当费时。

还有一种更为简便的方法就是采用内置自启动模型(self-starting Models),此时我们只需要指定函数形式,而不需要指定参数初始值。本例的logistic函数所对应的selfstarting函数名为SSlogis

pop.mod2 <- nls(population ~ SSlogis(year,phi1,phi2,phi3),data=USPop)

二、判断拟合效果

非线性回归模型建立后需要判断拟合效果,因为有时候参数最优化过程会捕捉到局部极值点而非全局极值点。最直观的方法是在原始数据点上绘制拟合曲线。

library(ggplot2)

p <- ggplot(USPop,aes(year, population))