Ошибка в нелинейных наименьших сквалах в кривых R - Logistic и Gompertz - PullRequest
Новая
       40

Ошибка в нелинейных наименьших сквалах в кривых R - Logistic и Gompertz

0 голосов
/ 09 сентября 2018

Я работаю над моделью для переменной y, в которой я намерен использовать время в качестве пояснительной переменной. Я выбрал Gompertz и логистическую кривую в качестве кандидатов, но когда я пытаюсь оценить коэффициенты (используя как nls, так и nls2), я получаю разные ошибки (сингулярность или коэффициент шага уменьшаются ниже 'minFactor'). Я был бы очень признателен за любую помощь. Вот мой код и версия информационного объекта deput.

Я выбрал начальные значения в соответствии с критериями в http://www.metla.fi/silvafennica/full/sf33/sf334327.pdf

library(nls2)

> dput(info)
structure(list(y = c(0.308, 0.279, 0.156, 0.214, 0.224, 0.222, 
0.19, 0.139, 0.111, 0.17, 0.155, 0.198, 0.811, 0.688, 0.543, 
0.536, 0.587, 0.765, 0.667, 0.811, 0.587, 0.617, 0.586, 0.633, 
2.231, 2.202, 1.396, 1.442, 1.704, 2.59, 2.304, 3.026, 2.7, 3.275, 
3.349, 3.936, 9.212, 8.773, 6.431, 6.983, 7.169, 9.756, 10.951, 
13.938, 14.378, 18.406, 24.079, 28.462, 51.461, 46.555, 39.116, 
43.982, 41.722), t = 1:53), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -53L))

summary(gomp_nls <- nls2(y ~ alpha*exp(-beta*exp(-gamma*t)), 
                 data = info, 
                 start = list(alpha = 40, beta = 4.9, gamma = 0.02), 
                 algorithm = "default")
        )

summary(logist_nls <- nls2(y ~ alpha/(1+beta*exp(-gamma*t)), 
                          data = info, 
                          start = list(alpha = 40, beta = 128, gamma = 0.02), 
                          algorithm = "default"))
        )

Буду признателен за любую помощь

1 Ответ

0 голосов
/ 12 сентября 2018

Алгоритм "default" для nls2 должен использовать nls. Вы хотите указать "brute-force" или один из других алгоритмов для поиска начального значения. Начальным значением должен быть фрейм данных из двух строк, чтобы он заполнил гиперкуб, определенный таким образом с потенциальными начальными значениями.

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

Если вы обнаружите, что результат, возвращаемый nls2, находится на границе определенной вами области, увеличьте область и попробуйте снова. (Вам может не понадобиться этот шаг, если возвращаемое начальное значение достаточно хорошее).

Наконец, запустите nls с начальными значениями, которые вы нашли. 

library(nls2)

## 1

fo1 <- y ~ alpha*exp(-beta*exp(-gamma*t))
st1 <- data.frame(alpha = c(10, 100), beta = c(1, 100), gamma = c(0.01, 0.20))
fm1.0 <- nls2(fo1, data = info, start = st1, algorithm = "brute-force")

fm1 <- nls(fo1, data = info, start = coef(fm1.0))

## 2

fo2 <- y ~ alpha/(1+beta*exp(-gamma*t))
st2 <- data.frame(alpha = c(10, 1000), beta = c(1, 10000), gamma = c(0.01, 0.20))
fm2.0 <- nls2(fo2, data = info, start = st2, algorithm = "brute-force")

fm2 <- nls(fo2, data = info, start = coef(fm2.0))

# plot both fits

plot(y ~ t, info)
lines(fitted(fm1) ~ t, info, col = "blue")
lines(fitted(fm2) ~ t, info, col = "red")

screenshot

Примечание

Обратите внимание, что для показанных данных эти две двухпараметрические экспоненциальные модели подходят достаточно хорошо, поэтому, если вас интересует только диапазон, где он экспоненциально растет, то это могут быть альтернативы для рассмотрения. (Первый ниже лучше, потому что коэффициенты больше похожи друг на друга. Второй может иметь проблемы с масштабированием.) 

fm3 <- nls(y ~ a * exp(b/t), info, start = c(a = 1, b = 1))
fm4 <- nls(y ~ a * t^b, info, start = c(a = .001, b = 6))
...