Начальный список параметров для nls и рациональной функции - PullRequest
Новая
       20

Начальный список параметров для nls и рациональной функции

0 голосов
/ 16 ноября 2018

Я хотел бы добавить рациональную функцию к моим данным:

data: 

 [1] 2.000000 3.000000 2.333333 1.750000 2.000000 1.833333 2.416667  1.916667
 [9] 1.750000 2.166667 2.116667 1.916667 1.944444 1.611111 1.722222 1.777778
[17] 1.877778 1.944444 1.958333 1.833333 2.041667 2.020833 1.908333 1.916667
[25] 1.733333 1.833333 1.800000 1.933333 1.893333 1.866667 1.888889 1.805556
[33] 1.833333 1.847222 1.822222 1.805556 1.833333 1.904762 1.880952 1.833333
[41] 1.804762 1.809524 1.708333 1.708333 1.750000 1.708333 1.683333 1.687500
[49] 1.611111 1.666667 1.648148 1.611111 1.611111 1.611111 1.650000 1.600000
[57] 1.650000 1.625000 1.630000 1.616667 1.469697 1.560606 1.590909 1.651515
[65] 1.651515 1.651515 1.513889 1.555556 1.625000 1.638889 1.647222 1.652778
[73] 1.679487 1.717949 1.705128 1.698718.

Модель, которую я хотел бы подогнать, следующая: 

Model <- function(t, a, b, c, d) {
                 (a + b*t)/(1 + c*t + d*t^2)
 }

Я знаю, что сначала я должен дать стартовый список a, b, c ... для nls, но я действительно не знаю, как установить параметры.Поскольку я не эксперт, я нашел в этом документе http://www.css.cornell.edu/faculty/dgr2/teach/R/R_rat.pdf полезное руководство.Но в какой-то момент он говорит:

"При заданном наборе упорядоченных пар (ti, yi), где в общем случае проводятся повторные измерения при каждом значении t, параметры рациональной функции можно подобратьс помощью нелинейной оценки наименьших квадратов, например, с помощью метода nls в R. Один у нас есть четыре параметра, мы можем вычислить значение t, при котором это максимизируется, путем вычисления первой производной .... ".

Хотя я не сообщаю дополнительные данные, у меня есть еще один столбец, который представляет время (целые числа от 1:76, представляющие годы).

Кто-нибудь может мне помочь?

Best

E.

Ответы [ 2 ]

0 голосов
/ 16 ноября 2018

Модель не указана полностью в вопросе, но предполагается, что модель приведена в приведенном ниже коде, и данные, воспроизводимые в примечании 2 ниже, если мы устанавливаем c = d = 0, то это линейная модель, поэтому мы можем использовать коэффициенты из линейная модель подходит в качестве начальных значений: 

fm1 <- lm(y ~ t)
st2 <- list(a = coef(fm1)[[1]], b = coef(fm1)[[2]], c = 0, d = 0)
fm2 <- nls(y ~ Model(t, a, b, c, d), start = st2)

дает: 

> fm2
Nonlinear regression model
  model: y ~ Model(t, a, b, c, d)
   data: parent.frame()
        a         b         c         d 
2.5097712 0.6038808 0.3205409 0.0008663 
 residual sum-of-squares: 1.684

Number of iterations to convergence: 16 
Achieved convergence tolerance: 8.029e-06

Глядя на графически подгонку: 

# model is shown in red. See Note 1 for fm4 (blue) and fm5 (green) models.
plot(y ~ t)
lines(fitted(fm2) ~ t, col = "red")
lines(fitted(fm4) ~ t, col = "blue")
lines(fitted(fm5) ~ t, col = "green")
legend("topright", c("fm2", "fm4", "fm5"), col = c("red", "blue", "green"), lty = 1)

screenshot

Примечание 1

Ниже приведена другая модель, которая подходит почти так же, но использует только 3 параметра. См. Синюю линию на графике выше. 

fm3 <- lm(log(y) ~ log(t))
st4 <- list(a = coef(fm3)[[1]], b = 0, c = coef(fm3)[[2]])
fm4 <- nls(y ~ exp(a + b/t + c*log(t)), start = st4)

> fm4
Nonlinear regression model
  model: y ~ exp(a + b/t + c * log(t))
   data: parent.frame()
      a       b       c 
 0.9845 -0.1767 -0.1157 
 residual sum-of-squares: 1.685

Number of iterations to convergence: 4 
Achieved convergence tolerance: 2.625e-06

Также эта модель не так уж плоха. Он использует только два параметра, он линейный в них и имеет остаточную сумму квадратов 1.728837 по сравнению с 1.684 для модели fm2 и 1.685 для модели fm4. См. Зеленую линию на графике выше. 

fm5 <- lm(y ~ log(t))

> fm5

Call:
lm(formula = y ~ log(t))

Coefficients:
(Intercept)       log(t)  
     2.4029      -0.1793  

> deviance(fm5)
[1] 1.728837

Примечание 2 

y <- c(2, 3, 2.333333, 1.75, 2, 1.833333, 2.416667, 1.916667, 1.75, 
2.166667, 2.116667, 1.916667, 1.944444, 1.611111, 1.722222, 1.777778, 
1.877778, 1.944444, 1.958333, 1.833333, 2.041667, 2.020833, 1.908333, 
1.916667, 1.733333, 1.833333, 1.8, 1.933333, 1.893333, 1.866667, 
1.888889, 1.805556, 1.833333, 1.847222, 1.822222, 1.805556, 1.833333, 
1.904762, 1.880952, 1.833333, 1.804762, 1.809524, 1.708333, 1.708333, 
1.75, 1.708333, 1.683333, 1.6875, 1.611111, 1.666667, 1.648148, 
1.611111, 1.611111, 1.611111, 1.65, 1.6, 1.65, 1.625, 1.63, 1.616667, 
1.469697, 1.560606, 1.590909, 1.651515, 1.651515, 1.651515, 1.513889, 
1.555556, 1.625, 1.638889, 1.647222, 1.652778, 1.679487, 1.717949, 
1.705128, 1.698718)

t <- seq_along(y)
0 голосов
/ 16 ноября 2018

Полагаю, вы можете сделать что-то вроде 

nls(y ~ (a1 + b1*times) / (1 + c1*times + d1*times^2))

, где y - данные, которые вы указали выше, а times=1:76.Я добавил 1 к параметрам, потому что nls не распознал c как параметр, а скорее как функцию c().

Однако, когда я запускаю это, я получаю ошибкуsingular gradient и предлагает инициализировать начальные значения чем-то отличным от 1 (по умолчанию).Вы можете указать начальные значения с аргументом start = list("a1"=0.1, "b1"=0.1, "c1"=0.1, "d1"=0.1), но это, похоже, не помогло.Возможно, у вас есть лучшее представление о том, какими должны быть начальные значения?

...