Блог

Я пытаюсь подогнать экспоненциальную модель формы y = a*e^(x*b) в R к данным, которые я запустил ниже, используя nls() функцию. Я читал здесь и в других местах, что мне нужно ввести в модель разумные параметры для a и b , но способ получения этих начальных параметров кажется очень изменчивым, и многие из предложенных способов были методами, которые я не смог найти подробностей. В одном предложении говорилось, что вы можете скопировать свои данные в таблицу Excel, подогнать свою модель к графику и настроить параметры до тех пор, пока они не будут достаточно хорошо соответствовать данным. Ну, я зашел в Google Sheets, вставил диаграмму, основанную на приведенных ниже данных, а затем выбрал Customize > Series > Trendline (exponential) , и она выдала мне формулу 5,51e ^ 0,015x . Являются ли эти допустимые значения, которые я мог бы использовать в качестве начальных параметров? Эффективно ли Google Sheets их создает, или мне нужно использовать метод tinkering или попробовать что-то еще? Я снова и снова перечитывал важность выбора правильных начальных значений, поэтому буду признателен за любую помощь по этому вопросу. Мое образование не охватывало нелинейные модели.

 x       y
19.005  5.49
18.19   6
19.59   5.885
19.93   8.96
17.615  13.85
18.795  2.72
19.11   8.09
19.885  8.11
15.76   6.66
16.48   6.27
15.805  5.375
15.825  3.06
15.985  7.795
15.755  6.255
15.485  5.925
15.475  9.925
16.45   6.055
16.285  5.24
15.92   11.15
16.775  5.57
16.075  3.275
16.475  5.635
16.825  4.72
16.28   2.035
17.26   6.07
17.245  4.9
17.98   8.06
17.35   6.94
18.22   7.8
16.27   12.2
17.555  7.335
16.98   5.76
17.415  7.51
17.5    6.18
  

В случае данных, показанных в вопросе, мы могли бы просто использовать начальное значение b = 1. Если мы используем plinear алгоритм nls , то для линейного параметра не требуется никакого начального значения, a . В этом случае a не следует указывать в формуле, поскольку это уже подразумевается. Об этом будет сообщено как .lin в выходных данных. В первой строке кода мы сортируем DF данные DFs , чтобы облегчить последующее построение графика.

 DFs <- DF[order(DF$x), ]
fo <- y ~ exp(b * x)
fm <- nls(fo, DFs, start = list(b = 1), algorithm = "plinear")
  

Однако, если nls она не работает с некоторыми другими данными, то, поскольку y она строго положительна, мы можем взять журналы обеих сторон, чтобы получить линейную модель fm0 , которую можно использовать lm для получения начального значения b . Используйте plinear алгоритм, описанный выше, чтобы избежать необходимости указывать начальное значение для a . DFs и fo находятся сверху.

 fm0 <- lm(log(y) ~ x, DFs)
fm2 <- nls(fo, DFs, start = list(b = coef(fm0)[[2]]), algorithm = "plinear")
fm2

## Nonlinear regression model
##   model: y ~ exp(b * x)
##    data: DFs
##       b    .lin 
## 0.02819 4.10908 
##  residual sum-of-squares: 205.6
##
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 2.073e-06

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

Примечание

Входные данные в воспроизводимой форме:

 Lines <- "x       y
19.005  5.49
18.19   6
19.59   5.885
19.93   8.96
17.615  13.85
18.795  2.72
19.11   8.09
19.885  8.11
15.76   6.66
16.48   6.27
15.805  5.375
15.825  3.06
15.985  7.795
15.755  6.255
15.485  5.925
15.475  9.925
16.45   6.055
16.285  5.24
15.92   11.15
16.775  5.57
16.075  3.275
16.475  5.635
16.825  4.72
16.28   2.035
17.26   6.07
17.245  4.9
17.98   8.06
17.35   6.94
18.22   7.8
16.27   12.2
17.555  7.335
16.98   5.76
17.415  7.51
17.5    6.18"
DF <- read.table(text = Lines, header = TRUE)
  

Я написал аналогичный ответ на @G.Гротендика … общий смысл в том, что вы могли бы использовать Google sheets, или Excel, или какой-либо другой инструмент, чтобы получить начальные значения nls() , но в этом нет необходимости.

Я загрузил ваши данные с помощью dd <- read.table(header=TRUE, text="x yn19.005 5.49 ...")

логарифмическая линейная модель

 m1 <- lm(log(y)~x,data=dd)
cc <- coef(m1)
## (Intercept)           x 
##   1.1404581   0.0399405 
  

Логарифмически-линейная модель на самом деле не дает такого же ответа, как nls() (хотя и достаточно хорошего для начальных значений), поскольку она делает разные неявные предположения о распределении ошибок (логарифмически-линейная предполагает, что дисперсия постоянна в логарифмическом масштабе, в то время как нелинейные подгонки предполагают, что дисперсия постоянна в логарифмическом масштабе).исходный масштаб, если не указано иное).

лог-ссылка GLM

Другой полезной альтернативой является подгонка обобщенной линейной модели с гауссовым откликом и ссылкой на журнал. На самом деле это соответствует точно такой же модели, как nls()

 m2 <- glm(y~x, family=gaussian(link="log"), data=dd)
cc2 <- coef(m2)
## (Intercept)           x 
##  1.41320439  0.02819136 
  

… что вы можете увидеть, если сравните эти значения коэффициентов с теми, которые nls() находят …

 s <- list(a=exp(cc2[[1]]), b= cc2[[2]])
m3 <- nls(y~a*exp(b*x), data=dd, start=s)
coef(m3)
  

Несколько разных способов делать прогнозы и визуализировать их…

 par(las=1,bty="l")
plot(y~x,data=dd)
lines(dd$x,exp(predict(m1)))
lines(dd$x,exp(cc[1])*exp(cc[2]*dd$x), col=2,lty=2)

lines(dd$x, predict(m2,type="response"), col=3)
lines(dd$x, exp(cc2[1])*exp(cc2[2]*dd$x), col=4,lty=2)
legend("topright",col=1:4,lty=rep(1:2,2),
       c("lm(log) predict","lm(log) formula","glm predict","glm formula"))
dev.off()