Экспоненциальный распад вписывается в г - PullRequest
0 голосов
/ 24 мая 2018

Я хотел бы подогнать функцию экспоненциального затухания в R к следующим данным:

data <- structure(list(x = 0:38, y = c(0.991744340878828, 0.512512332368168, 
0.41102449265681, 0.356621905557202, 0.320851602373477, 0.29499198506227, 
0.275037747162642, 0.25938850981822, 0.245263623938863, 0.233655093612007, 
0.224041426946405, 0.214152907133301, 0.207475138903635, 0.203270738895484, 
0.194942528735632, 0.188107106969046, 0.180926819430008, 0.177028560207711, 
0.172595416846822, 0.166729221891201, 0.163502461048814, 0.159286528409165, 
0.156110097827889, 0.152655498715612, 0.148684858095915, 0.14733605355542, 
0.144691873223729, 0.143118852619617, 0.139542186417186, 0.137730138713745, 
0.134353615271572, 0.132197800438632, 0.128369567159113, 0.124971834736476, 
0.120027536018095, 0.117678812415655, 0.115720611113327, 0.112491329844252, 
0.109219168085624)), class = "data.frame", row.names = c(NA, 
-39L), .Names = c("x", "y"))

Я пробовал подгонять по nls, но сгенерированная кривая не близка к фактическим данным.

введите описание изображения здесь

Было бы очень полезно, если бы кто-нибудь мог объяснить, как работать с такими нелинейными данными и найти функцию наилучшего соответствия.

Ответы [ 3 ]

0 голосов
/ 24 мая 2018

Двухэкспоненциальная модель подошла бы намного лучше, хотя все еще не идеально.Это будет означать, что у вас может быть два процесса распада одновременно.

fit <- nls(y ~ SSbiexp(x, A1, lrc1, A2, lrc2), data = data)
#A1*exp(-exp(lrc1)*x)+A2*exp(-exp(lrc2)*x)

plot(y ~x, data = data)
curve(predict(fit, newdata = data.frame(x)), add = TRUE)

resulting plot

Если погрешность измерения зависит от величины, вы можете использовать ее дляweighting.

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

0 голосов
/ 26 мая 2018

Попробуйте y ~ .lin / (b + x^c).Обратите внимание, что при использовании "plinear" линейный параметр .lin пропускается при указании формулы для nls, а также для него не указывается начальное значение.

Также обратите внимание, что параметры .lin и b приблизительно равны 1 при оптимальном, поэтому мы также можем попробовать модель с одним параметром y ~ 1 / (1 + x^c).Это форма однопараметрической логистической кривой выживания.AIC для этой модели с одним параметром хуже, чем для модели с тремя параметрами (сравните AIC(fm1) и AIC(fm3)), но модель с одним параметром все еще может быть предпочтительнее из-за ее экономии и того факта, что подбор визуально неотличим от модели 3.модель параметров.

opar <- par(mfcol = 2:1, mar = c(3, 3, 3, 1), family = "mono")

# data = data.frame with x & y col names; fm = model fit; main = string shown above plot
Plot <- function(data, fm, main) {
  plot(y ~ x, data, pch = 20)
  lines(fitted(fm) ~ x, data, col = "red")
  legend("topright", bty = "n", cex = 0.7, legend = capture.output(fm))
  title(main = paste(main, "- AIC:", round(AIC(fm), 2)))
}  

# 3 parameter model
fo3 <- y ~ 1/(b + x^c) # omit .lin parameter; plinear will add it automatically
fm3 <- nls(fo3, data = data, start = list(b = 1, c = 1), alg = "plinear")
Plot(data, fm3, "3 parameters")

# one parameter model
fo1 <- y ~ 1 / (1 + x^c)
fm1 <- nls(fo1, data, start = list(c = 1))
Plot(data, fm1, "1 parameter")

par(read.only = opar)

screenshot

AIC

Добавляя решения в другие ответы, мы можем сравнить значения AIC.Мы пометили каждое решение числом используемых им параметров (степени свободы были бы на единицу больше, чем это) и переработали решение log-log, чтобы использовать nls вместо lm и иметь LHS из y, так как одинне может сравнивать значения AIC моделей, имеющих разные левые части или использующих разные подпрограммы оптимизации, поскольку используемые константы логарифмического правдоподобия могут различаться:

    df     AIC
fm3  4 -329.35
fm1  2 -307.69
fm4  5 -215.96
fm2  3 -167.33
0 голосов
/ 24 мая 2018
data <- structure(list(x = 0:38, y = c(0.991744340878828, 0.512512332368168, 
0.41102449265681, 0.356621905557202, 0.320851602373477, 0.29499198506227, 
0.275037747162642, 0.25938850981822, 0.245263623938863, 0.233655093612007, 
0.224041426946405, 0.214152907133301, 0.207475138903635, 0.203270738895484, 
0.194942528735632, 0.188107106969046, 0.180926819430008, 0.177028560207711, 
0.172595416846822, 0.166729221891201, 0.163502461048814, 0.159286528409165, 
0.156110097827889, 0.152655498715612, 0.148684858095915, 0.14733605355542, 
0.144691873223729, 0.143118852619617, 0.139542186417186, 0.137730138713745, 
0.134353615271572, 0.132197800438632, 0.128369567159113, 0.124971834736476, 
0.120027536018095, 0.117678812415655, 0.115720611113327, 0.112491329844252, 
0.109219168085624)), class = "data.frame", row.names = c(NA, 
-39L), .Names = c("x", "y"))

# Do this because the log of 0 is not possible to calculate
data$x = data$x +1

fit = lm(log(y) ~ log(x), data = data)
plot(data$x, data$y)
lines(data$x, data$x ^ fit$coefficients[2], col = "red")

Это намного лучше, чем использование nls forumla.И при составлении графика подгонка кажется довольно неплохой.

...