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

lndjwyie  于 2023-04-03  发布在  其他
关注(0)|答案(2)|浏览(263)

我试图通过最小化最小二乘法来拟合微分方程模型，其中四个参数中的三个参数不能低于零。除了这个限制，我不希望应用任何进一步的约束。
下面是数据、模型和辅助函数：

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)
}

我尝试过同时使用FME：：modFit和stats：：optim，只限制约束参数的下限，但这会产生“non-finite finite-difference value”错误消息。

# Optimizing with limited lower bounds
fit <- FME::modFit(f = LS_modFit, p = get(i),
                     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 = get(i),
                     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")

你知道我哪里做错了吗？

r

2条答案

按热度按时间
zd287kbt

zd287kbt1#

你可以考虑使用DEoptim R包。我试过了，你的代码没有错误：

library(DEoptim)
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)
}

obj_DEoptim <- DEoptim(fn = LS_modFit, lower = c(0,0,-2,65), upper = c(0.02,0.02,1,400), 
                       time = t, falala = LVfn_dc, y = y0, y.values = y.data)
赞(0）分享  回复(0）举报  2023-04-03
db2dz4w8

db2dz4w82#

并非所有优化器都支持约束，所以FME对其他优化器使用转换。我看到两种不同的方法来实现它：
1.使用有限约束，例如100或1 e6
1.使用原生支持约束的优化器，例如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"))

创建于2023-04-01使用reprex v2.0.2

赞(0）分享  回复(0）举报  2023-04-03
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备2022004421号-2 NULL123 https://www.null123.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com