Как найти хорошие начальные значения для функции nls? - PullRequest
Новая
       10

Как найти хорошие начальные значения для функции nls?

7 голосов
/ 14 марта 2012

Я не понимаю, почему у меня не может быть функции nls для этих данных.Я пробовал с разными начальными значениями, и у меня всегда одна и та же ошибка.

Вот что я делал: 

expFct2 = function (x, a, b,c)
{
  a*(1-exp(-x/b)) + c  
}
vec_x <- c(77.87,87.76,68.6,66.29)
vec_y <- c(1,1,0.8,0.6)
dt <- data.frame(vec_x=vec_x,vec_y=vec_y)
ggplot(data = dt,aes(x = vec_x, y = vec_y)) +  geom_point() + 
     geom_smooth(data=dt, method="nls", formula=y~expFct2(x, a, b, c),
       se=F, start=list(a=1, b=75, c=-5)

У меня всегда есть эта ошибка: 

Error in method(formula, data = data, weights = weight, ...) : 
  singular gradient

Ответы [ 3 ]

9 голосов
/ 14 марта 2012

Подгонка трехпараметрической нелинейной модели к четырем точкам данных будет в любом случае довольно сложной задачей, хотя в этом случае данные ведут себя хорошо. Точка # 1 состоит в том, что ваше начальное значение для вашего параметра c (-5) было слишком далеко. Построение изображения кривой, соответствующей вашим начальным параметрам (см. Ниже), поможет вам понять это (так, если бы вы поняли, что полученная вами кривая будет варьироваться от минимума c до максимума c+a и диапазона ваши данные от 0,6 до 1 ...)

Однако, даже с лучшим начальным предположением, я обнаружил, что суетюсь с параметрами управления (то есть control=nls.control(maxiter=200)), за которыми следуют другие предупреждения - nls не известен своей надежностью. Поэтому я попробовал модель SSasympOff, которая реализует самозапускающуюся версию кривой, которую вы хотите подогнать. 

start1 <- list(a=1, b=75, c=-5)
start2 <- list(a=0.5, b=75, c=0.5)  ## a better guess

pfun <- function(params) {
  data.frame(vec_x=60:90,
             vec_y=do.call(expFct2,c(list(x=60:90),params)))
}
library(ggplot2)
ggplot(data = dt,aes(x = vec_x, y = vec_y)) +  geom_point() +
  geom_line(data=pfun(start1))+
  geom_line(data=pfun(start2),colour="red")+
  geom_smooth(data=dt, method="nls", formula=y~SSasympOff(x, a, b, c),
              se=FALSE)

Мой совет в целом состоит в том, что проще выяснить, что происходит, и исправить проблемы, если вы поместите nls вне из geom_smooth и построите кривую, которую хотите добавить, используя predict.nls. ..

В более общем смысле, способ получить хорошие начальные параметры состоит в том, чтобы понять геометрию подходящей функции и какие параметры определяют, какие аспекты кривой. Как я уже упоминал выше, c - это минимальное значение смещенной экспоненциальной кривой с насыщением, a - диапазон, а b - параметр масштаба (вы можете видеть, что когда x=b, кривая равна 1-exp(-1) или примерно 2/3 пути от минимума до максимума). Либо немного алгебры и исчисления (то есть с ограничениями), либо игра с функцией curve(), являются хорошими способами сбора этой информации.

8 голосов
/ 14 марта 2012

Это может быть записано с двумя линейными параметрами (.lin1 и .lin2) и одним нелинейным параметром (b), например: 

a*(1-exp(-x/b)) + c  
= (a+c) - a * exp(-x/b)
= .lin1 + .lin2 * exp(-x/b)

, где .lin1 = a+c и .lin2 = -a (поэтому a = - .lin2 и c = .lin1 + .lin2) Это позволяет нам использовать "plinear", который требует только указания начального значения для одного нелинейного параметра (устраняя проблему того, как установить начальные значения для других параметров) и который сходится, несмотря наначальное значение b=75, далекое от значения решения: 

nls(y ~ cbind(1, exp(-x/b)), start = list(b = 75), alg = "plinear")

Вот результат прогона, из которого мы можем видеть из размера .lin2, что проблема плохо масштабируется: 

> x <- c(77.87,87.76,68.6,66.29)
> y <- c(1,1,0.8,0.6)
> nls(y ~ cbind(1, exp(-x/b)), start = list(b = 75), alg = "plinear")
Nonlinear regression model
  model:  y ~ cbind(1, exp(-x/b)) 
   data:  parent.frame() 
         b      .lin1      .lin2 
 3.351e+00  1.006e+00 -1.589e+08 
 residual sum-of-squares: 7.909e-05

Number of iterations to convergence: 9 
Achieved convergence tolerance: 9.887e-07 
> R.version.string
[1] "R version 2.14.2 Patched (2012-02-29 r58660)"
> win.version()
[1] "Windows Vista (build 6002) Service Pack 2"

РЕДАКТИРОВАТЬ: добавлен пример запуска и комментарий по масштабированию.

2 голосов
/ 14 марта 2012

Я изо всех сил пытаюсь найти интерпретацию ваших параметров: a - наклон, b - скорость сходимости, a + c - предел, но само по себе, кажется, ничего не значит. После повторной настройки вашей функции проблема исчезает. 

f <- function (x, a,b,c) a + c * exp(-x/abs(b))
nls(y~f(x, a, b, c), data=dt, start=list(a=1, b=75, c=-5), trace=TRUE)

Однако значение c выглядит очень и очень высоким: возможно, поэтому модель изначально не сходилась. 

Nonlinear regression model
  model:  y ~ f(x, a, b, c) 
   data:  dt 
         a          b          c 
 1.006e+00  3.351e+00 -1.589e+08 
 residual sum-of-squares: 7.909e-05

Number of iterations to convergence: 9 
Achieved convergence tolerance: 2.232e-06

Вот еще одна, более разумная параметризация той же функции. 

g <- function (x, a,b,c) a * (1-exp(-(x-c)/abs(b)))
nls(y~g(x, a, b, c), data=dt, start=list(a=1, b=75, c=-5), trace=TRUE)

Nonlinear regression model
  model:  y ~ g(x, a, b, c) 
   data:  dt 
     a      b      c 
 1.006  3.351 63.257 
 residual sum-of-squares: 7.909e-05

Number of iterations to convergence: 10 
Achieved convergence tolerance: 1.782e-06
...