首页 / 技术问答 / R问答列表 / 问答详情

R 使用 FME：：modFit 或 optim 将模型拟合约束添加到参数子集

我试图通过最小化最小二乘来拟合一个微分方程模型，其中要优化的四个参数中有三个不能低于零.除了这一限制外，我不希望施加任何进一步的限制.

以下是数据、模型和帮助器函数:

library(deSolve)
library(FME)

# Data
y.data <- data.frame(P = c(57537.4049311277, 58610.5565091595, 59380.7326528789, 59831.1412677501, 59956.9401326381, 59770.9438648549, 59339.7636243282, 58743.8966682831, 58062.6867570216, 57372.6802888231, 56729.6957034719, 56118.5748350006, 55507.8890166606, 54867.5694355373, 54168.9851576162, 53385.1758149806, 52500.082193378, 51534.2686505716, 50516.106686327), 
SSL = c(5918.49121715316, 6223.93819671156, 6522.34271426757, 6817.47760344158, 7115.49425740418, 7423.50644959141, 7748.62985898112, 8097.51802623853, 8472.32658437055, 8872.9327092582, 9294.62369710465, 9730.87922007949, 10175.6848086032, 10623.7661461289, 11075.7535333951, 11538.2076285642, 12018.4500914707, 12518.681172246, 13039.7328357312),
S = c(61.8032786885246, 69.8469945355191, 71.1475409836066, 75.0163934426229, 65.6393442622951, 63.7049180327869, 60.0437158469945, 67.9344262295082, 66.5573770491803, 67.0327868852459, 57.6010928961749, 43.7158469945355, 57.224043715847, 65.1366120218579, 49.2131147540984, 30.7158469945355, 34.3715846994535, 12.3661202185792, 27.4972677595628))

# Initial state conditions
y0 <- c(P=y.data[1,]$P, SSL=y.data[1,]$SSL, S=y.data[1,]$S)

# Initial free parameter values
a <- 0.00075
b <- 0.004464286
rs <- 0.01
ks <- 100
p0 <- c(a=a,b=b,rs=rs,ks=ks)

# Series along which to integrate
t <- 0:18
l <- length(t)

# Differential equation model
LVfn_dc <- function(time, y, param) {
  P = y[1]
  SSL = y[2]
  S = y[3]
  with(as.list(param), {
    dPdt <- ((-0.4101*4)*time^3) + ((18.058*3)*time^2) - ((297.9*2)*time) + 1511.2
    dSSLdt <- ((7.2468*2)*time) + 265.67      
    dSdt <- (rs*S)*((ks-S-(a*P)-(b*SSL))/ks)
    return(list(c(dPdt,dSSLdt,dSdt)))
  })
}

# Function to solve differential equation
sde <- function(param, time, f, y) {
  sol <- deSolve::ode(y = y, 
                      times = time, 
                      func = f, 
                      parms = param,
                      method = "rk4",
                      rtol = 1e-15, atol = 1e-15, maxsteps = 1e9)
  return(sol)
}

# Least squares function for optimization using FME::modFit
LS_modFit <- function(param, time, falala, y, y.values){
  # call sde function to compute y.theta(t) after initializing
  y.theta.all <- c()
  y.theta <- c()
  call.sde <- sde(param,time,falala,y)
  y.theta.all <- call.sde[,-1] # disregard t column
  # only keep y.theta(t.i) values for each data point t.i
  for(i in 1:length(time)){
    if(time[i]%%1 == 0){
      y.theta <- rbind(y.theta,y.theta.all[i,])
    }
  }
  error <- sum((y.values - y.theta)^2)
  return(error)
}

我try 使用FME::modFit和stats::OpTim，只限制受约束参数的下限，但这会产生"非有限差值"错误消息.

# Optimizing with limited lower bounds
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     upper=c(Inf,Inf,Inf,Inf), lower=c(0,0,-Inf,1),
                     method = "L-BFGS-B")

当我对所有四个参数使用非无限但宽泛的上下限时(在运行sde函数时不会产生无限值)，我仍然收到非有限类型的错误.

# Optimizing with upper and lower bounds for all parameters
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     upper=c(0.02,0.02,1,400), lower=c(0,0,-2,65),
                     method = "L-BFGS-B")

你能洞察到我哪里错了吗？

推荐答案

并非所有优化器本身都支持约束，因此FME对其他优化器使用转换.我认为有两种不同的方法可以让它发挥作用:

  1. 使用有限约束，例如100或1e6
  2. 使用本地支持约束的优化器，例如Port或bobyqa.

我还修改了以下内容:

  • 我不明白get(i)的目的，把它换成了p0.
  • 我没有使用循环，而是使用了预定义的函数modCost，并将time作为自变量添加到y.data.请注意，如果状态变量具有不同的标度，则定义加权方法非常重要
  • 将rk4替换为lsoda.
  • 我还建议减少rtol和rtol，因为这会使模拟变慢并且帮助不大，因为优化器有他们自己的容差.
library(deSolve)
library(FME)

# Data
y.data <- data.frame(
  time = 0:18,
  P = c(57537.4049311277, 58610.5565091595, 59380.7326528789, 59831.1412677501, 59956.9401326381, 59770.9438648549, 59339.7636243282, 58743.8966682831, 58062.6867570216, 57372.6802888231, 56729.6957034719, 56118.5748350006, 55507.8890166606, 54867.5694355373, 54168.9851576162, 53385.1758149806, 52500.082193378, 51534.2686505716, 50516.106686327),
  SSL = c(5918.49121715316, 6223.93819671156, 6522.34271426757, 6817.47760344158, 7115.49425740418, 7423.50644959141, 7748.62985898112, 8097.51802623853, 8472.32658437055, 8872.9327092582, 9294.62369710465, 9730.87922007949, 10175.6848086032, 10623.7661461289, 11075.7535333951, 11538.2076285642, 12018.4500914707, 12518.681172246, 13039.7328357312),
  S = c(61.8032786885246, 69.8469945355191, 71.1475409836066, 75.0163934426229, 65.6393442622951, 63.7049180327869, 60.0437158469945, 67.9344262295082, 66.5573770491803, 67.0327868852459, 57.6010928961749, 43.7158469945355, 57.224043715847, 65.1366120218579, 49.2131147540984, 30.7158469945355, 34.3715846994535, 12.3661202185792, 27.4972677595628)
)

# Initial state conditions
y0 <- c(P=y.data[1,]$P, SSL=y.data[1,]$SSL, S=y.data[1,]$S)

# Initial free parameter values
a <- 0.00075
b <- 0.004464286
rs <- 0.01
ks <- 100
p0 <- c(a=a,b=b,rs=rs,ks=ks)

# Series along which to integrate
t <- 0:18
l <- length(t)

# Differential equation model
LVfn_dc <- function(time, y, param) {
  P = y[1]
  SSL = y[2]
  S = y[3]
  with(as.list(param), {
    dPdt <- ((-0.4101*4)*time^3) + ((18.058*3)*time^2) - ((297.9*2)*time) + 1511.2
    dSSLdt <- ((7.2468*2)*time) + 265.67
    dSdt <- (rs*S)*((ks-S-(a*P)-(b*SSL))/ks)
    return(list(c(dPdt,dSSLdt,dSdt)))
  })
}

# Function to solve differential equation
sde <- function(param, time, f, y) {
  sol <- deSolve::ode(y = y,
                      times = time,
                      func = f,
                      parms = param,
                      method = "lsoda", # lsoda is more stable and faster than rk4
                      rtol = 1e-6, atol = 1e-6)
  return(sol)
}

# Least squares function for optimization using FME::modFit
LS_modFit <- function(param, time, falala, y, y.values){
  call.sde <- sde(param, time, falala, y)
  ## important! choose appropriate weighting method,
  ## depending on scale of variables
  ## options:  one of "none", "std", "mean"
  modCost(call.sde, y.values, weight="std")
}


# Optimizing with limited lower bounds
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     ## infinite "constraints" work only
                     ## for a subset of solvers, e.g. Port and bobyqa
                     upper=c(Inf,Inf,Inf,Inf), lower=c(0,0,-Inf,1),
                     ## workaround for the other optimizers
                     ## recommended: don't make limits too big
                     #upper=c(1e6,1e6,1e6,1e6), lower=c(0,0,-1e6,1),
                     method = "Port")

guess <- sde(p0, t, LVfn_dc, y0)
fitted <- sde(fit$par, t, LVfn_dc, y0)

plot(guess, fitted, obs=y.data, mfrow=c(1, 3))
legend("bottomleft", col=1:2, lty=1:2, legend=c("guess", "fitted"))

solution of the ODE system with fitted parameters

创建于2023-04-01年第reprex v2.0.2页

R相关问答推荐

有没有方法将paste 0功能与列表结合起来？

使用ggplot 2根据R中的类别排列Likert比例gplot

ggplot的轴标签保存在officer中时被剪切

修改用R编写的用户定义函数

如何在Chart_Series()中更改轴值的 colored颜色 ？

在R中无法读入具有Readxl和lApply的数据集

在ggplot中为不同几何体使用不同的 colored颜色 比例

如何写商，水平线，在一个单元格的表在R

如何用书面利率绘制geom_bar图

从R中的对数正态分布生成随机数的正确方法

在另一个包中设置断点&S R函数

将摘要图添加到facet_WRAP gglot的末尾

如何在反曲线图中更改X标签

数值型数据与字符混合时如何进行绑定

如何在R中创建条形图，使条形图在y轴上围绕0.5而不是0构建条形图？

TidyVerse中长度不等的列结合向量

我已经运行了几个月的代码的`Palette()`中出现了新的gglot错误

如何用不同长度的向量填充列表？

在shiny 表格中输入的文本在第一次后未更新

GgHighlight找不到它创建的列：`Highlight..1`->；`Highlight.....`

实用课程推荐
相关教程推荐
休闲君评测 |  TOOLFK工具网 |  古诗文网 |  良许Linux教程网 |  Python之禅 |  红色石头 |  入党申请书 | 

© 2020 Hi LearnFk 意见反馈 免责声明 关于我们 桂ICP备11002319号-3 友链&广告位+QQ: 1963612630