Подбор параметров нелинейного уравнения для соответствия экспериментальным данным

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

Уравнение

new_f <- function(a, t, t1, dt, m){
  1-exp(-a*((x-t1)/dt)^(m+1))
}

fit_d <- nls(y~new_f(a, t, -30, dt, m), start=list(a=1, dt=46, m=1))

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Может ли кто-нибудь помочь с этим пошаговым руководством для чайников? Любое предложение будет принята с благодарностью!.

Мои данные: t=θ и y=Y(θ)

t <- c(-30, -29, -28, -27, -26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -14, -13, -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 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, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89)

y <- c(0.0000, -0.0013, -0.0027, -0.0041, -0.0055, -0.0069, -0.0085, -0.0101, -0.0119, -0.0137, -0.0155, -0.0173, -0.0192, -0.0212, -0.0234, -0.0258, -0.0285, -0.0291, -0.0250, -0.0224, -0.0268, -0.0298, -0.0297, -0.0324, -0.0347, -0.0353, -0.0377, -0.0392, -0.0392, -0.0383, -0.0309, -0.0166, 0.0006, 0.0233, 0.0521, 0.0843, 0.1200, 0.1597, 0.1997, 0.2381, 0.2762, 0.3137, 0.3501, 0.3864, 0.4229, 0.4590, 0.4951, 0.5312, 0.5655, 0.5967, 0.6245, 0.6494, 0.6726, 0.6939, 0.7132, 0.7316, 0.7487, 0.7641, 0.7784, 0.7917, 0.8037, 0.8146, 0.8248, 0.8343, 0.8432, 0.8516, 0.8596, 0.8671, 0.8740, 0.8805, 0.8868, 0.8931, 0.8992, 0.9050, 0.9105, 0.9155, 0.9202, 0.9247, 0.9290, 0.9330, 0.9369, 0.9405, 0.9439, 0.9470, 0.9501, 0.9529, 0.9557, 0.9583, 0.9608, 0.9633, 0.9658, 0.9682, 0.9704, 0.9726, 0.9747, 0.9768, 0.9788, 0.9808, 0.9827, 0.9844, 0.9858, 0.9869, 0.9880, 0.9890, 0.9899, 0.9909, 0.9918, 0.9927, 0.9938, 0.9947, 0.9956, 0.9964, 0.9970, 0.9975, 0.9980, 0.9985, 0.9989, 0.9993, 0.9997, 1.0000)



Ответы (2)


Есть несколько проблем:

  • x в new_f должен быть t
  • проблема чрезмерно параметризована, так что параметры не поддаются идентификации. Существует бесконечное количество комбинаций параметров, удовлетворяющих этой задаче. В частности, a/dt^(m+1) можно заменить одним параметром. Чтобы исправить это, удалите dt, скажем, из начальных значений и зафиксируйте его на 46.

С этими изменениями код дает результат.

new_f <- function(a, t, t1, dt, m){
  1-exp(-a*((t-t1)/dt)^(m+1))
}

dt <- 46
fit_d <- nls(y~new_f(a, t, -30, dt, m), start = list(a = 1, m = 1))
fit_d
## Nonlinear regression model
##   model: y ~ new_f(a, t, -30, dt, m)
##    data: parent.frame()
##     a     m 
## 0.560 3.315 
## residual sum-of-squares: 0.3099
##
## Number of iterations to convergence: 13 
## Achieved convergence tolerance: 4.165e-06

plot(y ~ t)
lines(fitted(fit_d) ~ t, col = "red")

скриншот

person G. Grothendieck    schedule 22.02.2021
comment
Спасибо за выявление проблем в этом подходе. Но я не могу установить dt = 46. Мне нужно найти все три параметра a, m включая dt. Для моего подхода параметр (а) должен быть установлен между 2,3 и 6,9 (согласно литературным данным). В моем основном уравнении мне нужно найти тринадцать параметров (a1, a2, a3, m1, m2, m3, dt1, dt2, dt3, t2, t3, B1, B2), учитывая θ1 = t1 в начале первого сгорания фазы, и я могу получить B3 = 1 - (B1 + B2), что позволяет мне лучше соответствовать моим экспериментальным данным. Есть ли способ закодировать это или с чего я могу начать? - person Luis; 23.02.2021
comment
То, что вы пытаетесь сделать, невозможно. Параметры не идентифицируются. Это все равно, что сказать a+b=0 и что такое a и b. Вы можете выбрать один произвольно, а затем другой может быть определен как его отрицательный. Параметры в задаче в ответе не могут быть определены никаким алгоритмом, так как нет однозначного значения. Вы должны установить один из них, или наложить ограничение, или изменить формулу на что-то другое, чтобы их можно было идентифицировать, если вам нужны все три. Что касается границ, nls может принимать границы. См. ?nls - person G. Grothendieck; 23.02.2021
comment
Обратите внимание, что более сложная формула в вашем сообщении, не имеющая кода, может быть идентифицирована, поэтому проблема может заключаться в том, что она была слишком упрощена. - person G. Grothendieck; 23.02.2021
comment
Спасибо за расширение объяснения для лучшего понимания. Для параметра (а) могу ли я наложить ограничение между 2,3 и 6,9. Я рассмотрю (? NLS). Чтобы получить первое представление о R-кодировании, было предложено второе уравнение. Моя цель — закодировать мое первое уравнение и получить его параметры (ai, mi, ∆ti, βi, ti), чтобы понять его поведение, поскольку мое исследование сосредоточено на этих параметрах. - person Luis; 23.02.2021
comment
Связанного ограничения недостаточно. Это все равно, что сказать, что a + b = 0, а a находится между 0 и 1. Это все равно не фиксирует его. Тип ограничения, который будет работать, состоит в том, чтобы ограничить a равным, скажем, dt, а затем заменить dt на a, чтобы было 2 параметра. - person G. Grothendieck; 23.02.2021

Начиная с осмотра формы кривой можно наблюдать:

Для высоких значений x форма выглядит экспоненциальной, потому что довольно линейна в логарифмическом масштабе для y :

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

При низких значениях x форма имеет линейный вид:

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

Если мы попытаемся объединить две части кривой с помощью ступенчатой ​​функции Хевисайда H(x), результат будет не очень удовлетворительным, поскольку соединение не является гладким, как показано на следующем рисунке:

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

Чтобы получить гладкую кривую, можно заменить функцию Хевисайда логистической функцией:

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

Приведенные выше численные значения параметров не оптимизированы. Для лучшей подгонки необходимо продолжить нелинейную регрессию (это не было сделано). Основная трудность нелинейной регрессии заключается в том, чтобы установить хорошие предполагаемые значения параметров для запуска итеративного исчисления. Вот почему приведенные выше приблизительные значения могут быть очень полезными.

person JJacquelin    schedule 22.02.2021
comment
Спасибо за проверку кривой моих данных, но мне нужно использовать предложенное уравнение, так как мне нужно изучить поведение этих параметров (ai, mi, ∆θi, βi, θi) для моего исследования. Есть ли способ закодировать уравнение и получить 13 параметров (a1, a2, a3, m1, m2, m3, ∆θ1, ∆θ2, ∆θ3, θ2, θ3, β1, β2)? - person Luis; 23.02.2021