R nlsModel сингулярная матрица градиентов при ошибке оценки начальных параметров задача оценивания начальных параметров

Мне нужно сделать nls, подходящую для данных об урожае/годе (см. -a-strange-fit/64977765#64977765">Вопрос о том, что nls подходят для R - почему это так странно подходит?).

Я начинаю с этих данных и использую коэффициенты из обычной линейной модели в качестве начальных оценок параметров.

x <- c(1979L, 1980L, 1981L, 1982L, 1983L, 1984L, 1985L, 1986L, 1987L, 
1988L, 1989L, 1990L, 1991L, 1992L, 1993L, 1994L, 1997L, 1998L, 
2000L, 2001L, 2002L, 2003L, 2004L, 2012L, 2014L, 2018L)

y <- c(94.4, 100, 120.6, 97.3, 110, 110, 80, 110, 117, 115, 50, 120, 
68.4, 137, 83, 106, 124, 113, 95.8, 115.7, 60, 105, 60, 74, 95.7, 
100)

mod_lm <- lm(y~x)

Теперь я получаю эту ошибку с пакетом minpack.lm,

library(minpack.lm)

mod_nls <- nlsLM(y ~ a+x^b, start=list(a = mod_lm$coefficients[1],b=mod_lm$coefficients[2]))    

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

или этот, если я попробую пакет nls2.

library(nls2) 

mod_nls <- nls2(y ~ a+x^b, start=list(a = mod_lm$coefficients[1],b=mod_lm$coefficients[2]))

Error in numericDeriv(form[[3L]], names(ind), env) : 
 Missing value or an infinity produced when evaluating the model

Я предполагаю, что проблема заключается в плохой начальной оценке параметров. Как мне получить более точные оценки параметров на начальном этапе? Спасибо.


r nls non-linear-regression
person user8229029    schedule 24.11.2020    source источник
comment
Я получаю ту же ошибку, но не при использовании x <- 1:26. Также кажется, что он работает нормально, если не используются начальные параметры.   -  person Eric Krantz    schedule 24.11.2020
comment
Ваше решение работает нормально (с одной корректировкой). В любом случае, спасибо.   -  person user8229029    schedule 25.11.2020

Ответы (2)


arrow_upward
1
arrow_downward

Часть проблемы заключается в том, что в ваших данных нет сильной связи:

cor(x, y)
# [1] -0.219275
cor(x, y, method="spearman")
# [1] -0.2145792
cor(x, y, method="kendall")
# [1] -0.117833
plot(y~x)

Из корреляций и графика видно, что между x и y существует слабая отрицательная связь. По сути, нет убедительных доказательств того, что нелинейная подгонка будет лучше простой линейной регрессии. Линейная регрессия

Сюжет — это одна линия доказательств, а ранговые корреляции — другая. Если бы отношение было простым полиномом, то корреляции Спирмена и Кендалла были бы больше, чем корреляция Пирсона. Как указывает @Erik Krantz, позволяя функциям угадывать начальные значения (1, 1), работает лучше, чем использование коэффициентов линейной регрессии.

Вот график линии, описываемой вашими начальными условиями, по сравнению с данными:

plot(x, coef(mod_lm)[1] + x^coef(mod_lm)[2], ylim=c(0, 1000), type="l")
points(x, y)

Начало подгонки

Просто взглянув на график, вы увидите, что первое значение слишком велико и что начальное значение около 100 было бы лучше:

mod_nls <- nlsLM(y ~ a+x^b, start=list(a = 100, b = coef(mod_lm)[2]))
mod_nls
# Nonlinear regression model
#   model: y ~ a + x^b
#    data: parent.frame()
#       a     b.x 
# 98.1638 -0.1306 
#  residual sum-of-squares: 12142
# 
# Number of iterations to convergence: 13 
# Achieved convergence tolerance: 1.49e-08
person dcarlson    schedule 24.11.2020
comment
Да, я понимаю, что в этой ситуации подойдет простая линейная подгонка, но я не могу посмотреть на 1000 этих графиков, чтобы определить это, поэтому мне нужен способ точно автоматизировать это. - person user8229029; 25.11.2020
comment
Первой проверкой может быть сравнение корреляции Спирмена/Кендалла с корреляцией Пирсона. - person dcarlson; 25.11.2020

arrow_upward
0
arrow_downward

Мы можем использовать nls2, чтобы получить приблизительное соответствие, и мы видим, что это почти горизонтальная линия. Использование коэффициента для b, показанного ниже, дает приблизительно x ^ b = 1, поэтому второй член является постоянным и может быть включен в a, и мы точно так же можем использовать модель y ~ a, и в этой новой сокращенной модели a оценивается как mean(y). Дело в том, что любое достаточно большое отрицательное значение b даст этот результат, поэтому модель чрезмерно параметризована, и мы можем сократить ее до члена a.

library(nls2)

fm <- nls2(y ~ a + x^b, start = list(a = range(y), b = c(-10, 10)), alg = "brute")
fm
## Nonlinear regression model
##   model: y ~ a + x^b
##    data: parent.frame()
##      a      b 
## 99.714 -4.286 
##  residual sum-of-squares: 12179
##
## Number of iterations to convergence: 64 
## Achieved convergence tolerance: NA

plot(y ~ x)
lines(fitted(fm) ~ x, col = "red")

скриншот

person G. Grothendieck    schedule 25.11.2020