В R, как выполнить нелинейную оптимизацию по методу наименьших квадратов, которая включает решение дифференциальных уравнений?
Мой первоначальный вопрос был «Реализация алгоритма оптимизации отражающей области доверия в R». Однако, на пути создания воспроизводимого примера (спасибо @Ben за его совет), я понимаю, что моя проблема в том, что в Matlab одна функцияlsqnonlin это хорошо (то есть нет необходимости выбирать хорошее начальное значение, быстро) достаточно для большинства случаев, которые у меня есть, в то время как в R нет такой универсальной функции. Различный алгоритм оптимизации хорошо работает в разных случаях. Различные алгоритмы достигают разных решений. Причиной этого может быть не то, что алгоритмы оптимизации в R уступают алгоритму отражения области доверия в Matlab, это также может быть связано с тем, как R обрабатывает автоматическое дифференцирование. Эта проблема возникает на самом деле из-за прерванной работы два года назад. В то время, профессор Джон К Нэш, один из авторов пакетаoptimx предположил, что для Matlab было проделано достаточно много работы для автоматического дифференцирования, что может быть причиной того, что Matlab lsqnonlin работает лучше, чем функции / алгоритмы оптимизации в R. Я не могу понять это с моими знаниями.
В приведенном ниже примере показаны некоторые проблемы, с которыми я столкнулся (грядут более воспроизводимые примеры). Чтобы запустить примеры, сначала запуститеinstall_github("KineticEval","zhenglei-gao"), Вам необходимо установить пакетmkin и его зависимости, а также, возможно, потребуется установить кучу других пакетов для различных алгоритмов оптимизации.
В основном я пытаюсь решить нелинейные задачи подгонки кривой наименьших квадратов, как описано в функции Matlablsqnonlin документация (http://www.mathworks.de/de/help/optim/ug/lsqnonlin.html). Кривые в моем случае моделируются набором дифференциальных уравнений. Я объясню немного больше с примерами. Я попробовал алгоритмы оптимизации, в том числе:
nls.lmЛевенбург-МарквардтПорт отnlm.inbL-BGFS-B отoptimSpg отoptimxsolnp пакетаRsolnpЯ также попробовал несколько других, но не показывал здесь.
Краткое изложение моих вопросовЕсть ли в R надежная функция / алгоритм для использования, напримерlsqnonlin в Matlab, который может решить мой тип нелинейных задач наименьших квадратов? (Я не мог найти один.)Какова причина того, что для простого случая разные оптимизации достигают разных решений?Что делаетlsqnonlin превосходит функции в R? Алгоритм, отражающий область доверия, или другие причины?Есть ли лучшие способы решить мой тип проблем, сформулировать по-другому? Может быть, есть очень простое решение, но я просто не знаю его.Пример 1: простой случайСначала я дам коды R, а потом объясню.
ex1 <- mkinmod.full(
Parent = list(type = "SFO", to = "Metab", sink = TRUE,
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 195, fixed = 0,lower = 0,upper = Inf),
FF = list(ini = c(.1),fixed = c(0),lower = c(0),upper = c(1)),
time=c(0.0,2.8, 6.2, 12.0, 29.2, 66.8, 99.8, 127.5, 154.4, 229.9, 272.3, 288.1, 322.9),
residue = c( 157.3, 206.3, 181.4, 223.0, 163.2, 144.7,85.0, 76.5, 76.4, 51.5, 45.5, 47.3, 42.7),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)),
Metab = list(type = "SFO",
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 0, fixed = 1,lower = 0,upper = Inf),
residue =c( 0.0, 0.0, 0.0, 1.6, 4.0, 12.3, 13.5, 12.7, 11.4, 11.6, 10.9, 9.5, 7.6),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
)
ex1$diffs
Fit <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
Fit[[i]] <- mkinfit.full(ex1,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
}
names(Fit) <- alglist
kinplot(Fit[[2]])
(lapply(Fit, function(x) x$par))
unlist(lapply(Fit, function(x) x$ssr))
Вывод из последней строки:
L-BFGS-B Marq Port spg solnp
5735.744 4714.500 5780.446 5728.361 4714.499
За исключением «Marq» и «solnp», другие алгоритмы не достигли оптимального. Кроме того, для метода spg (также для других методов, таких как bobyqa) требуется слишком много вычислений функций для такого простого случая. Более того, если я изменю начальное значение и сделаюk_Parent=0.0058 (оптимальное значение для этого параметра) вместо случайного выбора0.1«Марк» уже не может найти оптимальный вариант! (Код приведен ниже). У меня также были наборы данных, где "solnp" не находит оптимального. Однако, если я используюlsqnonlin в Matlab я не столкнулся с какими-либо трудностями в таких простых случаях.
ex1_a <- mkinmod.full(
Parent = list(type = "SFO", to = "Metab", sink = TRUE,
k = list(ini = 0.0058,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 195, fixed = 0,lower = 0,upper = Inf),
FF = list(ini = c(.1),fixed = c(0),lower = c(0),upper = c(1)),
time=c(0.0,2.8, 6.2, 12.0, 29.2, 66.8, 99.8, 127.5, 154.4, 229.9, 272.3, 288.1, 322.9),
residue = c( 157.3, 206.3, 181.4, 223.0, 163.2, 144.7,85.0, 76.5, 76.4, 51.5, 45.5, 47.3, 42.7),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)),
Metab = list(type = "SFO",
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 0, fixed = 1,lower = 0,upper = Inf),
residue =c( 0.0, 0.0, 0.0, 1.6, 4.0, 12.3, 13.5, 12.7, 11.4, 11.6, 10.9, 9.5, 7.6),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
)
Fit_a <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
Fit_a[[i]] <- mkinfit.full(ex1_a,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
}
names(Fit_a) <- alglist
lapply(Fit_a, function(x) x$par)
unlist(lapply(Fit_a, function(x) x$ssr))
Теперь вывод из последней строки:
L-BFGS-B Marq Port spg solnp
5653.132 4866.961 5653.070 5635.372 4714.499
Я объясню, что я оптимизирую здесь. Если вы запустили приведенный выше сценарий и видите кривые, мы используем двухсекционную модель с реакциями первого порядка для описания кривых. Дифференциальные уравнения для выражения модели находятся вex1$diffs:
Parent
"d_Parent = - k_Parent * Parent"
Metab
"d_Metab = - k_Metab * Metab + k_Parent * f_Parent_to_Metab * Parent"
В этом простом случае из дифференциальных уравнений мы можем вывести уравнения для описания двух кривых. Оптимизируемыми параметрами являются$M_0,k_p, k_m, c=\mbox{FF_parent_to_Met} $ с ограничениями$M_0>0,k_p>0, k_m>0, 1> c >0$.
$
\begin{split}
y_{1j}&= M_0e^{-k_pt_i}+\epsilon_{1j}\\
y_{2j} &= cM_0k_p\frac{e^{-k_mt_i}-e^{-k_pt_i}}{k_p-k_m}+\epsilon_{2j}
\end{split}
$
Поэтому мы можем подогнать кривую без решения дифференциальных уравнений.
BCS1.l <- mkin_wide_to_long(BCS1)
BCS1.l <- na.omit(BCS1.l)
indi <- c(rep(1,sum(BCS1.l$name=='Parent')),rep(0,sum(BCS1.l$name=='Metab')))
sysequ.indi <- function(t,indi,M0,kp,km,C)
{
y <- indi*M0*exp(-kp*t)+(1-indi)*C*M0*kp/(kp-km)*(exp(-km*t)-exp(-kp*t));
y
}
M00 <- 100
kp0 <- 0.1
km0 <- 0.01
C0 <- 0.1
library(nlme)
result1 <- gnls(value ~ sysequ.indi(time,indi,M0,kp,km,C),data=BCS1.l,start=list(M0=M00,kp=kp0,km=km0,C=C0),control=gnlsControl())
#result3 <- gnls(value ~ sysequ.indi(time,indi,M0,kp,km,C),data=BCS1.l,start=list(M0=M00,kp=kp0,km=km0,C=C0),weights = varIdent(form=~1|name))
## Coefficients:
## M0 kp km C
## 1.946170e+02 5.800074e-03 8.404269e-03 2.208788e-01
Таким образом, прошедшее время составляет почти 0, и оптимальное значение достигается. Однако у нас не всегда есть этот простой случай. Модель может быть сложной, и необходимо решение дифференциальных уравнений. Смотрите пример 2
Пример 2, сложная модельЯ работал над этим набором данных давным-давно, и у меня нет времени, чтобы закончить сам запуск следующего сценария. (Вам может понадобиться несколько часов, чтобы закончить бег.)
data(BCS2)
ex2 <- mkinmod.full(Parent= list(type = "SFO",to = c( "Met1", "Met2","Met4", "Met5"),
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 100,fixed = 0,lower = 0,upper = Inf),
FF = list(ini = c(.1,.1,.1,.1),fixed = c(0,0,0,0),lower = c(0,0,0,0),upper = c(1,1,1,1))),
Met1 = list(type = "SFO",to = c("Met3", "Met4")),
Met2 = list(type = "SFO",to = c("Met3")),
Met3 = list(type = "SFO" ),
Met4 = list(type = "SFO", to = c("Met5")),
Met5 = list(type = "SFO"),
data=BCS2)
ex2$diffs
Fit2 <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
Fit2[[i]] <- mkinfit.full(ex2,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
}
names(Fit) <- alglist
(lapply(Fit, function(x) x$par))
unlist(lapply(Fit, function(x) x$ssr))
Это пример, где вы увидите предупреждающие сообщения, такие как:
DLSODA- At T (=R1) and step size H (=R2), the
corrector convergence failed repeatedly
or with ABS(H) = HMIN
In above message, R =
[1] 0.000000e+00 2.289412e-09
Многие из методов, используемых в решателях Matlab Optimization Toolbox, основаны на областях доверия. Согласно странице просмотра задач CRAN, только пакетыдоверять, trustOptim, minqa внедрить методы, основанные на доверительном регионе. Тем не мение,trust а такжеtrustOptim требуют градиента и гессиана.bobyqa вminqa кажется не тот, кого я ищу. Исходя из моего личного опыта, алгоритм отражения области доверия в Matlab часто работает лучше, чем алгоритмы, которые я пробовал в R. Поэтому я попытался найти аналогичную реализацию этого алгоритма в R.
Я задал связанный вопрос здесь:Функция R для поиска функции
Ответ Мэтью Плурд дает функциюlsqnonlin с тем же именем функции в Matlab, но не имеет реализованного алгоритма отражения области доверия. Я отредактировал старый и задал новый вопрос здесь, потому что я думаю, что ответ Мэтью Плурда в целом очень полезен для пользователей R, которые ищут функцию.
Я снова сделал поиск и мне не повезло. Есть ли еще какие-нибудь функции / пакеты, которые реализуют аналогичные функции Matlab. Если нет, то мне интересно, разрешено ли мне переводить функцию Matlab непосредственно в R и использовать ее для своих собственных целей.