通过使用 optim() 最小化残差平方和优化参数来拟合模型以观察数据
Fitting model to observed data by obtimising parameters using optim() minimising residual sum of squares
我正在尝试将这个非常简单的 4 种线性 Lotka-Volterra 竞争模型拟合到观察到的数据，但是由于某些原因，当我尝试 optim() 函数时，某些与 deSolve 有关的东西似乎失败了.
# Data
data <- data.frame(Cod = c(0.1966126, 0.1989563, 0.2567677, 0.3158896, 0.4225435, 0.7219856,
1.0570824, 0.7266830, 0.6286763, 0.6389475),
Herring = c(1.988372, 2.788014, 3.397138, 2.557245, 2.627013, 3.045617,
3.161002, 3.531306, 3.432021, 3.617174),
Sprat = c(2.030273, 3.480469, 3.009277, 1.895996, 2.457520, 1.991211, 2.350098,
2.118164, 1.693359, 1.869141),
Flounder = c(0.4758220, 0.4425532, 0.4185687, 0.4967118, 0.7102515, 0.5733075,
0.7404255, 0.5996132, 0.6235977, 0.7187621))
# Model formulation
LLV <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
db1.dt = b1*(r1+a11*b1+a12*b2+a13*b3+a14*b4)
db2.dt = b2*(r2+a22*b2+a21*b1+a23*b3+a24*b4)
db3.dt = b3*(r3+a33*b3+a31*b1+a32*b2+a34*b4)
db4.dt = b4*(r4+a44*b4+a41*b1+a42*b2+a43*b3)
list(c(db1.dt, db2.dt, db3.dt, db4.dt))
})
}
# Model input and simulation
# Model input
params <- c(r1 = -0.342085, r2 = 0.6855681, r3 = 2.757769, r4 = 0.9744113,
a11 = -1.05973762, a12 = 0.09577309, a13 = -0.01915480, a14 = 1.36098939,
a21 = 0.17533326, a22 = -0.32247342, a23 = 0.03111628, a24 = 0.30212711,
a31 = 0.5303516, a32 = -0.4869761, a33 = -0.3194882, a34 = -1.5089027,
a41 = 0.004418133, a42 = 0.163716414, a43 = -0.237873378, a44 = -1.519158802)
ini <- c(b1 = data[1,1], b2 = data[1,2], b3 = data[1,3], b4 = data[1,4])
tmax <- 10
t <- seq(1,tmax,0.1)
# Results and first parameter guess is more or less okay
results <- deSolve::ode(y = ini, times = t, func = LLV, parms = params)
matplot(data, pch = 1)
matplot(x = results[,1], y = results[,-1], type = "l", add = TRUE)
在这里，我继续编写一个函数来最小化残差平方和，当包含在 optim() 中并使用上述初始参数猜测应该会产生我正在寻找的结果。
min.RSS <- function(data, params) {
output <- deSolve::ode(y = ini, times = t, func = LLV, parms = params)
predictions <- exp(output[,-1])
observations <- data
return(sum((predictions-observations)^2))
}
result <- optim(par = params, fn = min.RSS, data = data)
fit <- deSolve::ode(y = ini, times = t, func = LLV, parms = result$par)
matplot(x = fit[,1], y = fit[,-1], type = "l", lwd = 3, add = TRUE)
任何关于如何解决这个问题的想法将不胜感激。
对于那些感兴趣的人，我设法得到了一个涉及更改 ode 集成方法的解决方案。这是工作优化器：
# Optimising parameter fit
LVmse = function(parms) {
out = as.matrix(deSolve::ode(ini, 1:10, LLV, parms, method="rk")[,-1])
RSS = sum((spp-out)^2, na.rm = TRUE) # Minimising residual sum of squares
return(RSS)
}
optimres <- optim(par = params, fn = LVmse)
你的身材比较好，但你应该非常小心这个问题。我有点疯狂并使用（开发中）fitode package 来解决这个问题。我对模型进行了拟合并得到了更好的拟合，还尝试在我的最佳拟合周围随机选择 100 个随机变化的起点进行拟合。你的残差平方和是 1.19； fitode 第一次尝试达到 0.29，100 次拟合中最好的是 RSS=0.16。 但是：这些拟合高度不稳定。该图显示了数据的拟合和未来 5 个时间步的预测（1）您的拟合（虚线）； (2) fitode初始拟合（虚线）； (3) 其他100个fitode拟合（最佳拟合在0.05 RSS以内的是实心的，比那个更差的画得很淡）。
你可以看到样本外的预测大多是疯狂的。你的拟合实际上比一些更好的拟合 更 稳定 - 它在整个社区崩溃之前到达时间步骤 13 - 但底线是在这种情况下非常适合数据决不能保证一个合理的答案。看起来 100 次拟合中有一次到达预测时间序列的末尾而没有崩溃（这似乎是基于观察到的时间序列的合理合理的“常识”预测）。
为了可靠地拟合这些数据，您需要一个参数少得多的模型，或者以先验形式提供的外部信息，或者正则化 - 某种方法惩罚暗示 'wiggly' 确定性轨迹或相互作用 parameters/growth 不合理比率的拟合。
## remotes::install_github("parksw3/fitode")
library(fitode)
## data with tags for fitode
data2 <- setNames(data,paste0(names(data),"_obs"))
data2 <- data.frame(times=seq(nrow(data2)),data2)
## Model formulation (for fitode)
LV_model <- odemodel(
name="4-species LV",
model=list(
Cod ~ Cod*(r1+a11*Cod+a12*Herring+a13*Sprat+a14*Flounder),
Herring ~ Herring*(r2+a22*Herring+a21*Cod+a23*Sprat+a24*Flounder),
Sprat ~ Sprat*(r3+a33*Sprat+a31*Cod+a32*Herring+a34*Flounder),
Flounder ~ Flounder*(r4+a44*Flounder+a41*Cod+a42*Herring+a43*Sprat)
),
observation=list(
Cod_obs ~ ols(mean=Cod),
Herring_obs ~ ols(mean=Herring),
Sprat_obs ~ ols(mean=Sprat),
Flounder_obs ~ ols(mean=Flounder)
),
initial=list(
Cod ~ data2$Cod_obs[1],
Herring ~ data2$Herring_obs[1],
Sprat ~ data2$Sprat_obs[1],
Flounder ~ data2$Flounder_obs[1]
),
link=setNames(rep("identity",length(pars)),pars),
par= pars
)
## plot results
plotres <- function(p,ODEint="rk",lty=1,
dt=0.1,
tvec=seq(1,10,by=dt),...) {
par(las=1, bty="l")
res <- deSolve::ode(ini, tvec, LLV, p, method=ODEint)
matplot(res[,1],res[,-1],type="l",lty=lty,...)
return(invisible(res[,-1]))
}
f1 <- fitode(
LV_model,
data=data2,
start=params,
control=list(maxit=1e5,trace=1000)
)
## fitode with multistart
ranfit <- function(n,fit,range=0.5) {
##
rpars <- params*runif(length(params),1-range,1+range)
newfit <- try(update(fit, start=rpars))
return(newfit)
}
cl <- makeCluster(10)
clusterSetRNGStream(cl = cl, 101)
clusterExport(cl, c("params","LV_model","data2"))
clusterEvalQ(cl,invisible(library(fitode)))
system.time(
multifit <- parLapply(cl, 1:100, ranfit, fit=f1, tvec=tvec)
)
saveRDS(multifit,file="SO65440448_multifit.rds")
ivec <- seq_along(multifit)
ivec <- ivec[sapply(multifit,function(x) !inherits(x,"try-error"))]
coef <- pred <- vector("list", length=length(ivec))
ll <- conv <- rep(NA,length(ivec))
for (i in seq_along(ivec)) {
nf <- multifit[[ivec[i]]]
coef[[i]] <- coef(nf)
pp <- predict(nf, times=1:10)
pred[[i]] <- cbind(times=pp[[1]][,1],
do.call(cbind,lapply(pp,"[",-1)))
ll[i] <- logLik(nf)
conv[i] <- nf@mle2@details$convergence
}
par(las=1,bty="l")
matplot(pred[[1]][,1],pred[[1]][,-1],
type="n",lty=1,ylim=c(0,6),
xlab="time",ylab="density")
lthresh <- 0.05
for (i in 1:length(pred)) {
good <- ll[i]>(max(ll)-lthresh)
alpha <- if (good) 0.8 else 0.1
lwd <- if (good) 2 else 1
matlines(pred[[i]][,1],pred[[i]][,-1],lty=1,
col=adjustcolor(palette()[1:4],alpha.f=alpha),
lwd=lwd)
}
matpoints(data2[,1],data2[,-1],pch=16,cex=3)
plotres(optimres$par,add=TRUE, lwd=3,lty=2,dt=1)
plotres(coef(f1),add=TRUE, lwd=3,lty=3,dt=1)
我正在尝试将这个非常简单的 4 种线性 Lotka-Volterra 竞争模型拟合到观察到的数据，但是由于某些原因，当我尝试 optim() 函数时，某些与 deSolve 有关的东西似乎失败了.
# Data
data <- data.frame(Cod = c(0.1966126, 0.1989563, 0.2567677, 0.3158896, 0.4225435, 0.7219856,
1.0570824, 0.7266830, 0.6286763, 0.6389475),
Herring = c(1.988372, 2.788014, 3.397138, 2.557245, 2.627013, 3.045617,
3.161002, 3.531306, 3.432021, 3.617174),
Sprat = c(2.030273, 3.480469, 3.009277, 1.895996, 2.457520, 1.991211, 2.350098,
2.118164, 1.693359, 1.869141),
Flounder = c(0.4758220, 0.4425532, 0.4185687, 0.4967118, 0.7102515, 0.5733075,
0.7404255, 0.5996132, 0.6235977, 0.7187621))
# Model formulation
LLV <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
db1.dt = b1*(r1+a11*b1+a12*b2+a13*b3+a14*b4)
db2.dt = b2*(r2+a22*b2+a21*b1+a23*b3+a24*b4)
db3.dt = b3*(r3+a33*b3+a31*b1+a32*b2+a34*b4)
db4.dt = b4*(r4+a44*b4+a41*b1+a42*b2+a43*b3)
list(c(db1.dt, db2.dt, db3.dt, db4.dt))
})
}
# Model input and simulation
# Model input
params <- c(r1 = -0.342085, r2 = 0.6855681, r3 = 2.757769, r4 = 0.9744113,
a11 = -1.05973762, a12 = 0.09577309, a13 = -0.01915480, a14 = 1.36098939,
a21 = 0.17533326, a22 = -0.32247342, a23 = 0.03111628, a24 = 0.30212711,
a31 = 0.5303516, a32 = -0.4869761, a33 = -0.3194882, a34 = -1.5089027,
a41 = 0.004418133, a42 = 0.163716414, a43 = -0.237873378, a44 = -1.519158802)
ini <- c(b1 = data[1,1], b2 = data[1,2], b3 = data[1,3], b4 = data[1,4])
tmax <- 10
t <- seq(1,tmax,0.1)
# Results and first parameter guess is more or less okay
results <- deSolve::ode(y = ini, times = t, func = LLV, parms = params)
matplot(data, pch = 1)
matplot(x = results[,1], y = results[,-1], type = "l", add = TRUE)
在这里，我继续编写一个函数来最小化残差平方和，当包含在 optim() 中并使用上述初始参数猜测应该会产生我正在寻找的结果。
min.RSS <- function(data, params) {
output <- deSolve::ode(y = ini, times = t, func = LLV, parms = params)
predictions <- exp(output[,-1])
observations <- data
return(sum((predictions-observations)^2))
}
result <- optim(par = params, fn = min.RSS, data = data)
fit <- deSolve::ode(y = ini, times = t, func = LLV, parms = result$par)
matplot(x = fit[,1], y = fit[,-1], type = "l", lwd = 3, add = TRUE)
任何关于如何解决这个问题的想法将不胜感激。
对于那些感兴趣的人，我设法得到了一个涉及更改 ode 集成方法的解决方案。这是工作优化器：
# Optimising parameter fit
LVmse = function(parms) {
out = as.matrix(deSolve::ode(ini, 1:10, LLV, parms, method="rk")[,-1])
RSS = sum((spp-out)^2, na.rm = TRUE) # Minimising residual sum of squares
return(RSS)
}
optimres <- optim(par = params, fn = LVmse)
你的身材比较好，但你应该非常小心这个问题。我有点疯狂并使用（开发中）fitode package 来解决这个问题。我对模型进行了拟合并得到了更好的拟合，还尝试在我的最佳拟合周围随机选择 100 个随机变化的起点进行拟合。你的残差平方和是 1.19； fitode 第一次尝试达到 0.29，100 次拟合中最好的是 RSS=0.16。 但是：这些拟合高度不稳定。该图显示了数据的拟合和未来 5 个时间步的预测（1）您的拟合（虚线）； (2) fitode初始拟合（虚线）； (3) 其他100个fitode拟合（最佳拟合在0.05 RSS以内的是实心的，比那个更差的画得很淡）。
你可以看到样本外的预测大多是疯狂的。你的拟合实际上比一些更好的拟合 更 稳定 - 它在整个社区崩溃之前到达时间步骤 13 - 但底线是在这种情况下非常适合数据决不能保证一个合理的答案。看起来 100 次拟合中有一次到达预测时间序列的末尾而没有崩溃（这似乎是基于观察到的时间序列的合理合理的“常识”预测）。
为了可靠地拟合这些数据，您需要一个参数少得多的模型，或者以先验形式提供的外部信息，或者正则化 - 某种方法惩罚暗示 'wiggly' 确定性轨迹或相互作用 parameters/growth 不合理比率的拟合。
## remotes::install_github("parksw3/fitode")
library(fitode)
## data with tags for fitode
data2 <- setNames(data,paste0(names(data),"_obs"))
data2 <- data.frame(times=seq(nrow(data2)),data2)
## Model formulation (for fitode)
LV_model <- odemodel(
name="4-species LV",
model=list(
Cod ~ Cod*(r1+a11*Cod+a12*Herring+a13*Sprat+a14*Flounder),
Herring ~ Herring*(r2+a22*Herring+a21*Cod+a23*Sprat+a24*Flounder),
Sprat ~ Sprat*(r3+a33*Sprat+a31*Cod+a32*Herring+a34*Flounder),
Flounder ~ Flounder*(r4+a44*Flounder+a41*Cod+a42*Herring+a43*Sprat)
),
observation=list(
Cod_obs ~ ols(mean=Cod),
Herring_obs ~ ols(mean=Herring),
Sprat_obs ~ ols(mean=Sprat),
Flounder_obs ~ ols(mean=Flounder)
),
initial=list(
Cod ~ data2$Cod_obs[1],
Herring ~ data2$Herring_obs[1],
Sprat ~ data2$Sprat_obs[1],
Flounder ~ data2$Flounder_obs[1]
),
link=setNames(rep("identity",length(pars)),pars),
par= pars
)
## plot results
plotres <- function(p,ODEint="rk",lty=1,
dt=0.1,
tvec=seq(1,10,by=dt),...) {
par(las=1, bty="l")
res <- deSolve::ode(ini, tvec, LLV, p, method=ODEint)
matplot(res[,1],res[,-1],type="l",lty=lty,...)
return(invisible(res[,-1]))
}
f1 <- fitode(
LV_model,
data=data2,
start=params,
control=list(maxit=1e5,trace=1000)
)
## fitode with multistart
ranfit <- function(n,fit,range=0.5) {
##
rpars <- params*runif(length(params),1-range,1+range)
newfit <- try(update(fit, start=rpars))
return(newfit)
}
cl <- makeCluster(10)
clusterSetRNGStream(cl = cl, 101)
clusterExport(cl, c("params","LV_model","data2"))
clusterEvalQ(cl,invisible(library(fitode)))
system.time(
multifit <- parLapply(cl, 1:100, ranfit, fit=f1, tvec=tvec)
)
saveRDS(multifit,file="SO65440448_multifit.rds")
ivec <- seq_along(multifit)
ivec <- ivec[sapply(multifit,function(x) !inherits(x,"try-error"))]
coef <- pred <- vector("list", length=length(ivec))
ll <- conv <- rep(NA,length(ivec))
for (i in seq_along(ivec)) {
nf <- multifit[[ivec[i]]]
coef[[i]] <- coef(nf)
pp <- predict(nf, times=1:10)
pred[[i]] <- cbind(times=pp[[1]][,1],
do.call(cbind,lapply(pp,"[",-1)))
ll[i] <- logLik(nf)
conv[i] <- nf@mle2@details$convergence
}
par(las=1,bty="l")
matplot(pred[[1]][,1],pred[[1]][,-1],
type="n",lty=1,ylim=c(0,6),
xlab="time",ylab="density")
lthresh <- 0.05
for (i in 1:length(pred)) {
good <- ll[i]>(max(ll)-lthresh)
alpha <- if (good) 0.8 else 0.1
lwd <- if (good) 2 else 1
matlines(pred[[i]][,1],pred[[i]][,-1],lty=1,
col=adjustcolor(palette()[1:4],alpha.f=alpha),
lwd=lwd)
}
matpoints(data2[,1],data2[,-1],pch=16,cex=3)
plotres(optimres$par,add=TRUE, lwd=3,lty=2,dt=1)
plotres(coef(f1),add=TRUE, lwd=3,lty=3,dt=1)