R语言怎么为自定义的模型拟合参数？

s609078902

我想用 R 基于以下数据拟合函数 $y = \dfrac{k}{1+\dfrac{k-n}{n} \exp(-(ax+b))}$ 的参数 $$a, b, k, n$$ 该怎么做？平时都是直接用 R 的现有函数，现在遇到没有现成函数的情况就懵逼了。

试了用nls 和 nlsLM 函数都报错了，提示如下

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

查了下可能是初始值设置的不好，于是试了几组初始值，还是会报错，不知道设置初始值有没有什么技巧或原则。

data <- data.frame(
  x = 2:20,
  y = c(4.3, 13.3, 27.5, 50.4, 78.4, 109.5, 145, 193.9, 230, 280.4, 335, 377, 420, 463, 497, 522, 543, 575, 587)
)

Cloud2016

s609078902 看看函数帮助文档 ?nls 是否可以？

s609078902

Cloud2016 nls 和 nlsLM 两个函数刚刚都试了，都提示这个错误。查了一下可能是初始值选择的不好，换了几组初始值都还是报错，想请问下初始值的选择有没有什么技巧呢？

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

chuxinyuan

样本量可能有点少了。要估计的参数比较多，损失的自由度本来就多，样本量如果太少就可能出问题。可以考虑用已知的函数形式加上随机误差捏造更多的数据然后去验证。

s609078902

chuxinyuan 只有 4 个参数应该还好，之前试过用6次多项式也可以拟合出来，不确定二者之间是否有区别。其实用三次多项式就可以拟合的挺好的，现在想的是能不能用别的函数来拟合，后续比较下他们之间的稳定性（会有新的几组数据会加进来）。

数据的图形整体上是个S型，所以我找了这个函数，但是拟合不出来，请问能不能推荐些其它函数？

yihui

s609078902 数据的图形整体上是个S型

既然是 S 型，那第一反应应该是 Logistic 曲线吧。如果用 Logistic 函数，那么连初始值都不需要设置，因为它的梯度矩阵是已知的，SSlogis() 函数可以自动帮忙计算大致靠谱的初始值。

x = 2:20
y = c(4.3, 13.3, 27.5, 50.4, 78.4, 109.5, 145, 193.9, 230, 280.4, 335, 377, 420, 463, 497, 522, 543, 575, 587)
nls(y ~ SSlogis(x, Asym, xmid, scal))

Nonlinear regression model
  model: y ~ SSlogis(x, Asym, xmid, scal)
   data: parent.frame()
  Asym   xmid   scal 
607.40  11.50   2.89 
 residual sum-of-squares: 1560

Number of iterations to convergence: 0 
Achieved convergence tolerance: 2.821e-06

又及：以后帖子里发代码和数据时，请以人类可以直接复制运行的方式发出来。你这数据里的空格让我摔了好几跤才明白过来它不是普通空格。光是把数据折腾进来这一步，可能就足以劝退想帮忙的人了。

s609078902

yihui 又及：以后帖子里发代码和数据时，请以人类可以直接复制运行的方式发出来。你这数据里的空格让我摔了好几跤才明白过来它不是普通空格。光是把数据折腾进来这一步，可能就足以劝退想帮忙的人了。

不好意思，没注意到数据格式的问题，已修改。虽然意识到要附上数据，但是忽略了数据的形式问题。

yihui 既然是 S 型，那第一反应应该是 Logistic 曲线吧。

我这个其实也算是 logistic 曲线，不过应该算是变种吧。原本的 logistic 曲线的值域是在 [0, 1] 内，为了解决这个问题才加入了 k 和 n 两个参数，现在看来好像不需要。我早上刚看到您的回复，等我研究下再反馈，感谢yihui大大。

fenguoerbian

我感觉还是你选择的这个函数的形式的问题，导致了梯度矩阵出现了奇异。另外我想问下这个函数形式里，n不是指样本量吧？就是恰好你选择了用这个字母来代表这个位置上的参数是吗？

如果只是为了S形曲线的拟合的话，像生长曲线，或者其他的一些logistic类型的曲线应该是有一些现成的工具的。比方说我搜到的这个页面。这个页面里的一些S形曲线的拟合其实也是可以用nls来做的。但是里面的这个参数化的形式和你的这个总有些区别，尤其是k和n的部分。

s609078902

fenguoerbian 另外我想问下这个函数形式里，n 不是指样本量吧？就是恰好你选择了用这个字母来代表这个位置上的参数是吗？

n 不是值样本量，只是刚好用 n 代表这个参数，这个函数我是参考了这篇文章。

fenguoerbian 如果只是为了 S 形曲线的拟合的话，像生长曲线，或者其他的一些 logistic 类型的曲线应该是有一些现成的工具的。比方说我搜到的这个页面。

刚刚才看到您的回复，等我研究下再来反馈，感谢感谢。

fenguoerbian

s609078902
那应该就是你的参数化的问题，你看你的参考文献里，它的logistic曲线的exp里面是-rt，而你这边多了一个b出来。这导致了你的求解过程中会出现的奇异的问题。

确定一条logistic曲线可以靠四个参数：整条曲线的最小值、最大值、整条曲线在x轴的哪个位置时对应着纵轴上最大最小值之间的中点，以及在该位置的斜率是多少。

你的参考文献里其实是个三参数的版本，因为它给定了最小值是0，剩下的k，n，r这三个参数可以实现和logistic曲线剩余那三个参数的一一对应（换算）关系。但是你的参数化版本多了这个b之后，这就不再是一个一一对应的关系，也就是说对于同一条logistic回归曲线，你的这个参数化方式有无穷多种结果都可以得到那条曲线。所以不管你数据量有多少，你这个参数化总会得到奇异收敛方面的报错。

再直观一点来看的话，你可以看到你的参考文献里说了，它采用这样的参数化的形式，是因为他想表达一个在x = 0处，这个曲线的高度在n/k。而按你这样的写法，就成了在x = -b/a处，曲线的高度在n/k。而这个-b/a，你把分子分母同时乘以任意的非零的数，结果都不变。也就是说对同一条logistic曲线，你的参数化方法有无穷多种不同的表示方式，所以当然估不出来。

yihui

fenguoerbian 壮士，受俺一拜。分析地太清楚了。

四个参数让我想起冯诺依曼的大象：https://en.wikipedia.org/wiki/Von_Neumann%27s_elephant （即：别说一条简单的 S 曲线了，四个参数都足以拟合一头大象）

s609078902

fenguoerbian 原来如此！大神，请受小弟一拜。我对 Logistic 曲线的了解还是太肤浅了，当时增加参数 b 是为了让曲线更加灵活，没想到反而导致出现了奇异的问题。你发的文章我还在看，同时结合 R 里函数的帮助文档在学习中。

Liechi

yihui
维基百科对这个故事的介绍跟原版有些出入。