r - 有没有办法在 R 中使用 optim 优化 ODE 的初始值？

Question

我想在使用时优化一个值（ numcaseswk0）。然而，这个值不是 ODE 的参数 - 只是一个初始值。Roptim

下面显示的代码尝试执行此操作，但优化过程不断失败，只产生用户提供的上限。我怀疑这可能是因为它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？

标签： roptimization

Answer 1

我认为您的模型规范有问题。优化器运行正确（在技术意义上），最优参数为零。

请看评论，有什么改变：

• 更大的惩罚值
• 将成本值重新调整到可行范围
• 对您的成本函数进行一些手动测试调用
• 一些可能有助于调试的 cat() 函数
• 另请注意，“点”语法（L.、R. 等）不是必需的

以下可能仍然不是您想要的，但希望有助于使其运行。祝你好运！

    library(deSolve)

    ### ODE FUNCTION

    HAVODE <- function(t, states, parameters){
      with(as.list(c(states,parameters)),
           {
             N    <- S + L + Z + I + R
             #cat("N=", N, "\n")
             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
      #cat("x=", x, "\n")

      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.)
      #cat(states, "\n")

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

      ## tpe: set it to a really large value
      if (!is.finite(NLLK)){
        NLLK = 0.1 * .Machine$double.xmax    # if NLLK is infinity, replace by a large number
      }

      ## tpe: re-scale return value to a numerically feasible range
      return(NLLK * 1e-10)  
    }

    ## vector of starting values
    x0 <- 2
    ## set lower and upper bounds for these variables
    upper <- 10
    lower <- 0

    ## Call the cost function with optim

    calib_parameters <- c(135722, 9.2088, 2.6047, 0.47, 3.930, 7.21, 0.094, 0.517)

    ## tpe: reproducible comparison data
    set.seed(42)
    observed <- abs(rnorm(50, mean=6, sd=3))

    ## test manually
    calib_function(1, calib_parameters, observed)
    calib_function(0, calib_parameters, observed)
    calib_function(10, calib_parameters, observed)

    ## tpe: we see that zero *is* the best among these

    ## optimize automatically
    optimization_results <- optim(par=x0, lower = lower,  upper = upper, 
                                  method = 'L-BFGS-B', fn = calib_function, 
                                  parameters = calib_parameters,  
                                  observed. =  observed,
                                  control=list(trace=TRUE))


    optimization_results
    ## tpe: optimized par is again zero, that confirms the manual test