R nls奇异梯度

Question

R nls奇异梯度

9

我尝试搜索了其他关于这个问题的帖子，但没有一个解决办法适用于我。我有一个自然实验的结果，并且想展示一个事件连续发生的次数符合指数分布。我的R shell代码如下:

f <- function(x,a,b) {a * exp(b * x)}
> x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27
> y
 [1] 1880  813  376  161  100   61   31    9    8    2    7    4    3    2    0
[16]    1    0    0    0    0    0    1    0    0    0    0    1
> dat2
    x    y
1   1 1880
2   2  813
3   3  376
4   4  161
5   5  100
6   6   61
7   7   31
8   8    9
9   9    8
10 10    2
11 11    7
12 12    4
13 13    3
14 14    2
> fm <- nls(y ~ f(x,a,b), data = dat2, start = c(a=1, b=1)) 
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model
> fm <- nls(y ~ f(x,a,b), data = dat2, start = c(a=7, b=-.5)) 
Error in nls(y ~ f(x, a, b), data = dat2, start = c(a = 7, b = -0.5)) : 
  singular gradient
> fm <- nls(y ~ f(x,a,b), data = dat2, start = c(a=7,b=-.5),control=nls.control(maxiter=1000,warnOnly=TRUE,minFactor=1e-5,tol=1e-10),trace=TRUE) 
4355798 :   7.0 -0.5
Warning message:
In nls(y ~ f(x, a, b), data = dat2, start = c(a = 7, b = -0.5),  :
  singular gradient

请原谅格式不太好，这是我在此处的第一篇文章。x包含直方图的箱子，y包含每个箱子在该直方图中出现的次数。由于0计数箱子会扰乱指数回归，dat2截止到14，而且我只需要适应那前14个。那些计数超过14的箱子，我有生物学上的理由认为它们很特殊。最初遇到的问题是无穷大，但我并没有得到0值。按照另一篇帖子建议给出合理的初始值后，我遇到了奇异梯度错误。我看到的唯一其他帖子都有更多的变量，我尝试增加迭代次数，但没有成功。感谢您的任何帮助。

- sessmurda

2个回答

5

请查看minpack.lm包中的nlsLM函数。这是nls的更健壮版本，可以处理残差平方和为零的数据。

https://www.r-bloggers.com/a-better-nls/

- user3634351

请解释一下你提供的解决方案是如何解决这个问题的。外部链接的缺点是未来并不总是可访问的。 - Lepidopteron

1

针对我试图解决的问题，nlsLM()函数能够返回一个合理的解决方案，而nls()则会出现奇异梯度错误。请尝试下面的示例。抱歉格式有些混乱，我还在学习这种格式。nls()失败了，但nlsLM正确地恢复了参数。#模拟一些数据 set.seed(20160227) x<-seq(0,50,1) y<-((runif(1,10,20)x)/(runif(1,0,10)+x))+rnorm(51,0,1) y<-1(1-exp( -1*(x/2)^0.5) ) f <- function(x,a,b,c) {a*(1-exp( -1*(x/b)^c) )} nls(y~~f(x,a,b,c),start=list(a=1,b=1, c=1) ) nlsLM(y~~f(x,a,b,c),start=list(a=1,b=1, c=1) ) - user3634351

最佳答案，至少对于强非线性拟合来说。 - AlainD

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- G. Grothendieck · Accepted Answer

1) 线性化以获得起始值 您需要更好的起始值：

# starting values
fm0 <- nls(log(y) ~ log(f(x, a, b)), dat2, start = c(a = 1, b = 1))

nls(y ~ f(x, a, b), dat2, start = coef(fm0))

提供：

Nonlinear regression model
  model: y ~ f(x, a, b)
   data: x
        a         b 
4214.4228   -0.8106 
 residual sum-of-squares: 2388

Number of iterations to convergence: 6 
Achieved convergence tolerance: 3.363e-06

1a) 同样地，我们可以使用 lm 通过编写以下代码来获得初始值：

y ~ a * exp(b * x)

如

y ~ exp(log(a) + b * x)

然后取两者的对数，得到一个以log(a)和b线性相关的模型：

log(y) ~ log(a) + b * x

可以使用 lm 来解决：

fm_lm <- lm(log(y) ~ x, dat2)
st <- list(a = exp(coef(fm_lm)[1]), b = coef(fm_lm)[2])
nls(y ~ f(x, a, b), dat2, start = st)

提供：

Nonlinear regression model
  model: y ~ f(x, a, b)
   data: dat2
       a        b 
4214.423   -0.811 
 residual sum-of-squares: 2388

Number of iterations to convergence: 6 
Achieved convergence tolerance: 3.36e-06

1b) 我们也可以通过重新参数化来使其工作。在这种情况下，只要我们按照参数转换的方式转换初始值，a=1和b=1就可以起作用。

nls(y ~ exp(loga + b * x), dat2, start = list(loga = log(1), b = 1))

提供：

Nonlinear regression model
  model: y ~ exp(loga + b * x)
   data: dat2
  loga      b 
 8.346 -0.811 
 residual sum-of-squares: 2388

Number of iterations to convergence: 20 
Achieved convergence tolerance: 3.82e-07

因此b如所示，而a = exp(loga) = exp(8.346) = 4213.3

2) plinear 另一个更简单的可能性是使用alg="plinear"，在这种情况下不需要参数线性输入的起始值。在这种情况下，问题中b=1的起始值似乎已足够。

nls(y ~ exp(b * x), dat2, start = c(b = 1), alg = "plinear")

提供：

Nonlinear regression model
  model: y ~ exp(b * x)
   data: dat2
        b      .lin 
  -0.8106 4214.4234 
 residual sum-of-squares: 2388

Number of iterations to convergence: 11 
Achieved convergence tolerance: 2.153e-06