Блог

Я пытаюсь найти подходящие для моих данных. Но пока что не везло. Попробовал логарифмические, отличные от пакета drc .. но я уверен, что должен быть лучший, я просто не знаю, какой тип. На другой ноте — я был бы благодарен за совет о том, как вообще заниматься охотой на кривых.

 library(drc)  dflt;-structure(list(x = c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19,   20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,   36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,   52, 53, 54, 55), y = c(0.1066, -0.6204, -0.2028, 0.2621, 0.4083,   0.4497, 0.6343, 0.7762, 0.8809, 1.0029, 0.8089, 0.7845, 0.8009,   0.9319, 0.9414, 0.9505, 0.9323, 1.0321, 0.9381, 0.8975, 1.0929,   1.0236, 0.9589, 1.0644, 1.0411, 1.0763, 0.9679, 1.003, 1.142,   1.1049, 1.2868, 1.1569, 1.1952, 1.0802, 1.2125, 1.3765, 1.263,   1.2507, 1.2125, 1.2207, 1.2836, 1.3352, 1.1311, 1.2321, 1.4277,   1.1645), w = c(898, 20566, 3011, 1364, 1520, 2376, 1923, 1934,   1366, 1010, 380, 421, 283, 262, 227, 173, 118, 113, 95, 69, 123,   70, 80, 82, 68, 83, 76, 94, 101, 97, 115, 79, 98, 84, 92, 121,   97, 102, 93, 92, 101, 74, 124, 64, 52, 63)), row.names = c(NA,   -46L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(`47` = 47L), class = "omit"))    fit lt;- drm(data = df,y ~ x,fct=LL.4(), weights = w)  plot(fit)  

1) Если мы проигнорируем веса, то y = a b * x c/x^2, по-видимому, подходит и является линейным по коэффициентам, поэтому его легко подогнать. Это кажется наклонным вверх, поэтому мы начали с линии, но затем нам нужно было смягчить ее, поэтому мы добавили взаимный термин. Обратная квадратичная функция работала немного лучше, чем обычная обратная, основанная на остаточной сумме квадратов, поэтому мы переключились на это.

 fm lt;- lm(y ~ x   I(1 / x^2), df) coef(summary(fm)) ## Estimate Std. Error t value Pr(gt;|t|) ## (Intercept) 1.053856e 00 0.116960752 9.010341 1.849238e-11 ## x 4.863077e-03 0.002718613 1.788808 8.069195e-02 ## I(1/x^2) -1.460443e 02 16.518887452 -8.841049 3.160306e-11  

Коэффициент x-члена не имеет значения на уровне 5% — значение p составляет 8% в таблице выше, — поэтому мы можем удалить его, и он будет соответствовать почти так же хорошо, давая модель только с двумя параметрами. На графике ниже соответствие fm с 3 параметрами является сплошным, а соответствие fm2 с 2 параметрами-пунктирным.

 fm2 lt;- lm(y ~ I(1 / x^2), df)  plot(y ~ x, df) lines(fitted(fm) ~ x, df) lines(fitted(fm2) ~ x, df, lty = 2)  

2) Другой подход заключается в использовании двух прямых линий. Это все еще непрерывно, но имеет одну недифференцируемую точку в точке перехода. Модель имеет 4 параметра, перехваты и наклоны каждой линии. Ниже мы используем веса. Его преимущество заключается в очевидной мотивации, основанной на появлении данных. Точка разрыва на пересечении двух линий может иметь значение как точка перехода между начальным ростом с более высоким уклоном и последующим ростом с более низким уклоном.

 # starting values use lines fitted to 1st ten and last 10 points fm_1 lt;- lm(y ~ x, df, subset = 1:10) fm_2 lt;- lm(y ~ x, df, subset = seq(to = nrow(df), length = 10)) st lt;- list(a = coef(fm_1)[[1]], b = coef(fm_1)[[2]],  c = coef(fm_2)[[1]], d = coef(fm_2)[[2]])  fm3 lt;- nls(y ~ pmin(a   b * x, c   d * x), df, start = st, weights = w)  # point of transition X lt;- with(as.list(coef(fm3)), (a - c) / (d - b)); X ## [1] 16.38465 Y lt;- with(as.list(coef(fm3)), a   b * X); Y ## [1] 0.8262229  plot(y ~ x, df) lines(fitted(fm3) ~ x, df)  

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

Ее быстрый и грязный пример определяемой пользователем функции (с ростом пакетов). Это, безусловно, можно сделать аналогичным образом с drc.

 library("growthrates")  grow_userdefined lt;- function (time, parms) {  with(as.list(parms), {  y lt;- (K * y0)/(y0   (K - y0) * exp(-mumax * time))   shift   return(as.matrix(data.frame(time = time, y = y)))  }) }  fit lt;- fit_growthmodel(FUN=grow_userdefined,   p = c(y0 = -1, K = 1, mumax = 0.1, shift = 1),   time = df$x, y = df$y)  plot(fit) summary(fit)  

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