У меня есть некоторые трудности с получением конкретной кривой, подходящей для R, в то время как она отлично работает в коммерческой программе подбора кривой.
Формула, которой должны соответствовать данные:
y(t) = A * exp(-a*(t)) + B * exp(-b*(t)) - (A+B) * exp(-c*(t))
Поэтому для этого я хочу использовать нелинейную регрессию, встроенную в R. Я занимаюсь этим в течение дня, и просто не могу заставить его работать. Проблема полностью связана с начальными значениями, поэтому я использую NLS2 для поиска начальных значений методом перебора.
y <- c(0,0.01377,0.01400875,0.0119175,0.00759375,0.00512125,0.004175,0.00355375,
0.00308875,0.0028925,0.00266375)
t <- c(0,3,6,12,24,48,72,96,120,144,168)
df <- data.frame(t,y)
plot(t,y);
#Our model:
fo <- y ~ f1*exp(-k1*t)+f2*exp(-k2*t)-(f1+f2)*exp(-k3*t);
#Define the outer boundaries to search for initial values
grd <- data.frame(f1=c(0,1),
f2=c(0,1),
k1=c(0,2),
k2=c(0,2),
k3=c(0,0.7));
#Do the brute-force
fit <- nls2(fo,
data=df,
start = grd,
algorithm = "brute-force",
control=list(maxiter=20000))
fit
coef(fit)
final <- nls(fo, data=df, start=as.list(coef(fit)))
Он должен дать следующие значения:
f1 0.013866
f2 0.005364
k1 0.063641
k2 0.004297
k3 0.615125
Хотя даже при довольно высоких значениях итераций я получаю только ерунду. Я явно делаю что-то не так, но не вижу этого
Изменить на основе комментария @Roland:
Метод, который вы предлагаете с приближением k1-3 с линейным подходом, кажется, работает с некоторыми наборами данных, но не со всеми из них. Ниже приведен код, который я использую сейчас на основе ваших данных.
#Oral example:
y <- c(0,0.0045375,0.0066325,0.00511375,0.00395875,0.003265,0.00276,
0.002495,0.00231875);
t <- c(0,12,24,48,72,96,120,144,168);
#IV example:
#y <- c(0,0.01377,0.01400875,0.0119175,0.00759375,0.00512125,0.004175,
#0.00355375,0.00308875,0.0028925,0.00266375)
#t <- c(0,3,6,12,24,48,72,96,120,144,168)
DF <- data.frame(y, t)
fit1 <- nls(y ~ cbind(exp(-k1*t), exp(-k2*t), exp(-k3*t)), algorithm = "plinear", data = DF,
start = list(k1 = 0.002, k2 = 0.02, k3= 0.2))
k1_predict <-summary(fit1)$coefficients[1,1]
k2_predict <-summary(fit1)$coefficients[2,1]
k3_predict <-summary(fit1)$coefficients[3,1]
fo <- y ~ f1*exp(-k1*t)+f2*exp(-k2*t)-(f1+f2)*exp(-k3*t);
fit2 <- nls(fo, data = DF,
start = list(k1 = k1_predict, k2 = k2_predict, k3 = k3_predict, f1 = 0.01, f2 = 0.01))
summary(fit2);
plot(t,y);
curve(predict(fit2, newdata = data.frame(t = x)), 0, 200, add = TRUE, col = "red")
@ГРАММ. Гротендик. Подобно предложению Роланда, ваше предложение также превосходно в том смысле, что оно способно подогнать одни наборы данных, но затрудняется с другими. Приведенный ниже код основан на вашем вводе и завершается с сингулярной матрицей градиента.
#Oral example:
y <- c(0,0.0045375,0.0066325,0.00511375,0.00395875,0.003265,0.00276,
0.002495,0.00231875);
t <- c(0,12,24,48,72,96,120,144,168);
#IV example:
#y <- c(0,0.01377,0.01400875,0.0119175,0.00759375,0.00512125,0.004175,
#0.00355375,0.00308875,0.0028925,0.00266375)
#t <- c(0,3,6,12,24,48,72,96,120,144,168)
df <- data.frame(y, t)
grd <- data.frame(f1=c(0,1),
f2=c(0,1),
k1=c(0,2),
k2=c(0,2),
k3=c(0,0.7));
set.seed(123)
fit <- nls2(fo,
data=df,
start = grd,
algorithm = "random",
control = nls.control(maxiter = 100000))
nls(fo, df, start = coef(fit), alg = "port", lower = 0)
plot(t,y);
curve(predict(nls, newdata = data.frame(t = x)), 0, 200, add = TRUE, col = "red")
