在R，如何做这涉及到求解微分方程非线性最小二乘优化？(In R, how to do nonline

2019-07-19 05:38发布

站内文章 / 前端开发
24 0 0

更新与重复性的例子来说明我的问题

我原来的问题是“R中信赖域反射算法的实现”。 然而，在产生重复的例子（感谢@Ben他的意见）的方式，我意识到，我的问题是，在Matlab，一个功能lsqnonlin好（这意味着没有必要选择一个很好的起点值，快）足以满足大多数情况下，我有，而在R，有没有这样的一个人参加的所有功能。 不同optmization算法在不同的情况下，效果很好。 不同的算法得出不同的解决方案。 这背后的原因可能不是R中的优化算法不如在Matlab信赖域反射算法，它也可以作为R如何处理自动分化。 这个问题从两年前开始中断的工作其实出自。 当时，教授约翰·C·纳什，包optimx的作者之一提出，出现了相当多的工作Matlab的自动分化，这可能是Matlab的lsqnonlin性能比优化功能更好的原因/ R中的算法，我无法用我的知识来弄明白。

下面的例子说明了一些问题我也遇到过（重复性更好的例子来了）。 要运行的例子中，第一次运行install_github("KineticEval","zhenglei-gao") 您需要安装包mkin及其依赖，也可能需要安装很多其他的包不同的优化算法。

基本上我试图解决如MATLAB函数中所描述的非线性最小二乘曲线拟合问题lsqnonlin的文档（ http://www.mathworks.de/de/help/optim/ug/lsqnonlin.html ）。 在我的情况下的曲线是由一组微分方程的建模。 我将解释一点与实例。 优化算法，我曾尝试包括：

  • MARQ从nls.lm的Levenburg -马夸特
  • 从港口nlm.inb
  • L-BGFS-B从optim
  • 从SPG optimx
  • solnp包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”外，其他算法没有达到optimum.Besides，“抢断”的方法（还包括其他的方法，如“bobyqa”）需要这样一个简单的例子太多功能的评价。 此外，如果我改变了初始值，使k_Parent=0.0058 （最佳为参数值），而不是随机choosen 0.1 ，“MARQ”无法找到最佳的更多！ （代码下面提供）。 我也有过的数据集，其中“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优化工具箱求解器所使用的方法是基于信任的区域。 按照CRAN任务视图页面，只有包的信赖 ，trustOptim，minqa的信任，基于区域的方法来实现。 然而， trust和trustOptim需要梯度和麻袋。 bobyqa在minqa好像不是我要找的人。 从我个人的经验，在Matlab的信赖域反射算法相比，我在R.试过算法通常执行更好，所以我试图找到该算法在河的类似implemetation

我曾问一个相关的问题在这里： 一个R函数来搜索功能

马修Plourde提供的答案给出了一个功能lsqnonlin与Matlab中相同功能的名称，但不具有信赖域反射算法来实现。 我编辑的旧这里提出一个新的问题，因为我觉得马修Plourde的回答是，一般到谁正在寻找一个函数R用户来说非常实用。

我又做了搜索，没有运气。 是否还有其实现类似功能的MATLAB的一些功能/包在那里。 如果不是这样，我不知道是否允许我直接tranlate MATLAB函数为R，并用它为我自己的目的。

Answer 1:

在一般情况下，当只有看你的问题的标题，我想简单地使用建议FME包。 但是，这不是你的问题的角度和成功可能取决于你设置你的模型的方式。

对于你在你的例子（配件与几个转化产物的降解数据）显示问题的类型，我创建了mkin包作为一个便利的包装来FME这种类型的问题。 因此，让我们来看看如何mkin 0.9-29执行在这些情况下。 随着mkin，你只能使用FME提供的算法：

实施例1 

library(mkin)

ex1_data_wide = data.frame(
  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),
  Parent = 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),
  Metab = 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))

ex1_data = mkin_wide_to_long(ex1_data_wide, time = "time")

ex1_model = mkinmod(Parent = list(type = "SFO", to = "Metab"),
                    Metab = list(type = "SFO"))

algs = c("L-BFGS-B", "Marq", "Port")

times_ex1 <- list()
fits_ex1 <- list()
for (alg in algs) {
  times_ex1[[alg]] <- system.time(fits_ex1[[alg]] <- mkinfit(ex1_model, ex1_data,
                                                             method.modFit = alg))
}

times_ex1
unlist(lapply(fits_ex1, function(x) x$ssr))

因此，文伯格 - 马夸特在nls.lm和港口算法都找到自己的最低，与LM是要快得多： 

$`L-BFGS-B`
       User      System verstrichen 
      2.036       0.000       2.051 

$Marq
       User      System verstrichen 
      0.716       0.000       0.714 

$Port
       User      System verstrichen 
      2.032       0.000       2.030 

L-BFGS-B     Marq     Port 
5742.312 4714.498 4714.498

当我告诉mkin使用形成的部分，而不仅仅是率

ex1_model = mkinmod(Parent = list(type = "SFO", to = "Metab"),
                    Metab = list(type = "SFO"), use_of_ff = "max")

并使用你的初始值， 

for (alg in algs) {
  times_ex1[[alg]] <- system.time(fits_ex1[[alg]] <- mkinfit(ex1_model, ex1_data,
    state.ini = c(195, 0),
    parms.ini = c(f_Parent_to_Metab = 0.1, k_Parent = 0.0058, k_Metab = 0.1),
    method.modFit = alg))
}

所有三种算法找到相同的解决方案，甚至更快。 但是，如果我再关掉率和分数的转型在调用mkinfit（ transform_rates = FALSE, transform_fractions = FALSE ），我得到 

L-BFGS-B     Marq     Port 
5653.132 4714.498 5653.070

所以它似乎是相关的参数在内部转换的方式（当你给界限FME做到这一点为好）。 在mkin，我做明确的内部参数的变换等等都需要使用默认设置优化参数没有界限。

实施例2 

library(mkin)
library(KineticEval) # for the dataset BCS2
data(BCS2)

ex2_data = mkin_wide_to_long(BCS2, time = "time")

ex2_model = mkinmod(Parent = list(type = "SFO", to = paste0("Met", 1:5)),
                    Met1 = list(type = "SFO", to = c("Met3", "Met4")),
                    Met2 = list(type = "SFO", to = "Met3"),
                    Met3 = list(type = "SFO"),
                    Met4 = list(type = "SFO", to = "Met5"),
                    Met5 = list(type = "SFO"))

times_ex2 <- list()
fits_ex2 <- list()

for (alg in algs) {
  times_ex2[[alg]] <- system.time(fits_ex2[[alg]] <- mkinfit(ex2_model, ex2_data,
    method.modFit = alg))
}   

times_ex2
unlist(lapply(fits_ex2, function(x) x$ssr))

此外，LM北京时间最快的，但最低最低，就是通过端口发现： 

$`L-BFGS-B`
       User      System verstrichen 
     75.728       0.004      75.653 

$Marq
       User      System verstrichen 
      6.440       0.004       6.436 

$Port
       User      System verstrichen 
     51.200       0.028      51.180 

L-BFGS-B     Marq     Port 
485.3099 572.9635 478.4379

以前我总是建议LM，但最近我还发现，有时候会陷入局部极小，这取决于不明确的参数的初始值。 一个例子是谢弗07的数据，在最后的处理单元测试mkinfit在mkin包，称为test.mkinfit.schaefer07_complex_example 。

希望这是有用的，亲切的问候，

约翰内斯

PS：我发现这个问题时，我注意到您在KineticEval包添加的信赖域反射优化，功能lsqnonlin（）一个纯的R实现在github上，我是做对反射信任区域的搜索。



文章来源: In R, how to do nonlinear least square optimization which involves solving differential equations?
标签： r optimization nonlinear-optimization
举报

收藏的人(0)

Ta的文章 更多文章
登录 后发表评论
0条评论
还没有人评论过~