【发布时间】：2019-11-13 01:16:53
【问题描述】：
我想使用 optim 优化 R 中的值 (numcaseswk0)。然而，这个值不是 ODE 的参数 - 只是一个初始值。
下面显示的代码尝试执行此操作，但优化过程不断失败，只产生用户提供的上限。我怀疑这可能是因为numcaseswk0 不是 ODE 的参数之一。如果有人能指出我如何解决这个问题，我会很高兴。谢谢。
library(deSolve)
### ODE FUNCTION
HAVODE <- function(t, states, parameters){
with(as.list(c(states,parameters)),
{
N <- S + L + Z + I + R
dS <- -beta * S * (I/N)
dL <- beta * S * (I/N) - (1/durL)*L
dI <- (1/durL)*L +(1/durRel)*R - (1/durI)*I
dZ <- (1-propRelapse)*(1/durI)*I
dR <- propRelapse*(1/durI)*I - (1/durRel)*R
return(list(c(dS, dL, dI, dZ, dR)))
})
}
### COST FUNCTION
calib_function <- function(x, parameters,observed.){
## Variable to be optimized
numcaseswk0 <- x
initpop = parameters[1]
durL = parameters[2]
durI = parameters[3]
fracImmune = parameters[4]
durRel = parameters[5]
propRelapse = parameters[6]
probdetec = parameters[7]
beta = parameters[8]
## Starting values for states
S. = (1-fracImmune)*initpop
L. = numcaseswk0 # *** I want this value to be optimized
I. = 0
Z. = fracImmune*initpop
R. = 0
states = c(S=S., L= L. , I=I., Z=Z., R=R.)
## Parameters to be fed into ODE solver
parameters1 = c(durL = durL, durI = durI,durRel = durRel, propRelapse = propRelapse, beta = beta )
tspan = seq(0, length(observed.)+10);
# Run the ODE solver
result <- data.frame(ode(y = states, times = tspan, func = HAVODE, parms = parameters1))
# Calculating model response (number of detected incident cases)
IncDetec <- probdetec *((1/durL)*result[, 3] + (1/durRel)*result[, 6])
model_response <- IncDetec[-1][1:length(observed.)] # exclude initial week
# Calculate negative log likelihood of model responses
NLLK <- -sum(dpois(x = floor(model_response), lambda = observed., log = TRUE ))
if (NLLK == Inf){
NLLK = 999999 # if NLLK is infinity, replace by a large number
}
return(NLLK)
}
## vector of starting values
x0 <- 2
## set lower and upper bounds for these variables
upper <- 10
lower <- 1
## Call the cost function with optim
calib_parameters <- c(135722, 9.2088, 2.6047, 0.47, 3.930, 7.21, 0.094, 0.517)
optimization_results <- optim(par=x0, lower = lower, upper = upper, method = 'Brent', fn = calib_function, parameters = calib_parameters, observed. = abs(rnorm(50, mean=6, sd=3)))
运行上面的代码给出：
> optimization_results
$par
[1] 1.000001
$value
[1] 113463174144
$counts
function gradient
NA NA
$convergence
[1] 0
$message
NULL
optim 产生的估计值是提供的下限值 (lower=1)。您可能还会注意到没有函数评估。为什么优化不适用于numcaseswk0？
【问题讨论】：
-
您好像没有给我们minimal reproducible example？光看代码很难解决这个问题……
-
嗨@BenBolker，感谢您指出需要重新构建问题的事实。我已经做到了这一点，确保我提供了重现问题所需的尽可能多的内容。如果您能再看一下，我将不胜感激。
-
我又试了一次，但代码似乎仍然缺少一些部分。当我运行它时，我得到：
Error in seq.default(length(observed.), et, length.out = numsteps + 1) : object 'numsteps' not found -
@tpetzoldt，你是对的；我已经将
numsteps存储在 R 的内存中。但是，该代码现在已经过编辑，因此您在重新运行时不会遇到错误。感谢您的耐心等待。 -
@tpetzoldt 好吧，你又是对的，对此我很抱歉。已修复（现在是真实的）。
标签： r optimization