R优化多个参数

Question

我正在使用 R optim() 函数来估计优化用户定义函数的参数集，如下所示。但是 optim() 输出是：

Error in optim(pstart, llAgedepfn, method = "L-BFGS-B", upper = up, lower = lo) : L-BFGS-B needs finite values of 'fn'

请帮忙。完整的脚本如下所示：

dataM<-cbind(c(1.91,0.29,0.08,0.02,0.01,0.28,0.45,0.36,0.42,0.17,0.16,0.06,0.17,0.17,0.12),
               c(0.27,4.54,0.59,0.05,0.04,0.13,0.48,0.68,0.66,0.18,0.11,0.06,0.08,0.08,0.08),
               c(0.07,0.57,4.48,0.48,0.02,0.05,0.09,0.43,0.78,0.52,0.17,0.10,0.05,0.05,0.14),
               c(0.02,0.04,0.44,4.34,0.36,0.09,0.07,0.11,0.41,0.77,0.43,0.10,0.03,0.04,0.14),
               c(0.01,0.04,0.01,0.36,2.20,0.46,0.19,0.15,0.19,0.34,0.62,0.30,0.09,0.03,0.22),
               c(0.22,0.11,0.05,0.09,0.45,0.91,0.61,0.43,0.37,0.26,0.41,0.63,0.29,0.16,0.15),
               c(0.31,0.35,0.07,0.05,0.16,0.54,0.81,0.59,0.48,0.36,0.33,0.43,0.47,0.26,0.20),
               c(0.22,0.45,0.29,0.08,0.11,0.34,0.53,0.85,0.71,0.39,0.27,0.26,0.26,0.28,0.38),
               c(0.22,0.36,0.44,0.26,0.12,0.24,0.36,0.59,0.91,0.61,0.35,0.28,0.20,0.22,0.29),
               c(0.09,0.10,0.30,0.49,0.22,0.17,0.28,0.33,0.62,0.80,0.52,0.29,0.20,0.11,0.46),
               c(0.10,0.07,0.12,0.32,0.48,0.32,0.30,0.27,0.42,0.61,0.78,0.47,0.33,0.23,0.49),
               c(0.04,0.04,0.06,0.08,0.24,0.53,0.41,0.28,0.36,0.36,0.50,0.67,0.51,0.19,0.47),
               c(0.10,0.05,0.04,0.02,0.07,0.23,0.43,0.26,0.23,0.23,0.33,0.48,0.75,0.51,0.49),
               c(0.05,0.04,0.03,0.05,0.02,0.10,0.19,0.22,0.21,0.10,0.18,0.14,0.40,0.79,0.82),
               c(0.03,0.02,0.03,0.03,0.06,0.04,0.06,0.12,0.11,0.18,0.16,0.14,0.16,0.34,1.26)
)

NormCM <- dataM/eigen(CMWkday)$values[1] #Normalizing the contact mtrix - divide by the largest eigen value

w <- c(495,528,548,603,617,634,720,801,957,937,798,755,795,1016,2469) 

g2 <- c(770,622,726,559,410,547,564,472,399,397,340,308,337,91,84) 

h2 <- c(269,426,556,430,271,284,303,207,194,181,126,106,74,24,23) 

z2 <- h2/g2

g1 <- c(774,527,665,508,459,539,543,492,402,412,365,342,213,146,152) 

h1 <- c(56,31,84,173,103,85,123,70,71,80,55,25,18,12,26) 
z1 <- h1/g1

#### Normal loglikelihood #########

llnormfn <- function(q) {  

  tol <- 1e-9
  final.size.start <- 0.8
  zeta <- rep(final.size.start, nrow(NormCM))
  last.zeta <- rep(0, nrow(NormCM))
  first.run <- T
  current.diff <- tol+1
  loglik <- 0

  while (current.diff > tol) {

    zeta <- 1-exp(-(q*(zeta%*%NormCM)))
    current.diff <- sum(abs(last.zeta-zeta))
    last.zeta <-zeta

  }
  mu <- c(zeta)

  zigma <- z1*(1-z1)/g1 + (z1+mu)*(1-(z1+mu))/g2

  logliknorm <- -sum((((z2-z1)-mu)**2)/2*zigma + 0.5*log(2*pi*zigma))

  return(logliknorm)

} 

pstart <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
up <- c(5,5,5,5,5,5,5,5,5,5,5,5,5,5,5)
lo <- c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1)
estm <- optim(pstart, llnormfn, method = "L-BFGS-B", upper = up, lower = lo )

Answer 1

您的 llnormfn 不是 return 范围内所有参数值的有限值。例如在上限：

> llnormfn(up)
[1] NaN
Warning message:
In log(2 * pi * zigma) : NaNs produced

因为zigma这里必须小于零。

如果你稍微限制一下范围，你最终可以找到它起作用的地方...

> llnormfn(up-2)
[1] NaN
Warning message:
In log(2 * pi * zigma) : NaNs produced
> llnormfn(up-3)
[1] 42.96818

让我们检查它在较低范围内的工作情况：

> llnormfn(lo)
[1] 41.92578

看起来不错。因此，要么您将上限设置在函数的计算有效范围之外，要么您的 llnormfn 函数中存在错误，或两者兼而有之，或者其他原因。

如果您执行运行减少上限的优化，您会收敛：

> estm <- optim(pstart, llnormfn, method = "L-BFGS-B", upper = up-3, lower = lo )
> estm
$par
 [1] 1.9042672 1.0891264 0.9916916 0.6208685 1.2413983 1.4822433 1.1243878
 [8] 1.5224263 1.3686933 1.4876350 1.6231518 2.0000000 2.0000000 2.0000000
[15] 2.0000000

$value
[1] 38.32182

$counts
function gradient 
      23       23 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

虽然您可能会注意到其中一些参数在上限值 (2.0) 这是一个警钟。

检查您的函数对其输入值的行为是否合理 - 尝试修复除一个以外的所有值并绘制 llnormfn 在改变一个时的行为方式。我只是快速浏览了一下，函数看起来一点也不流畅，有很多不连续性，所以我怀疑 BFGS 是否是一个好的优化方法。

例如，在 0.1 和 2 之间改变第五个参数：

> s = seq(0.1,2,len=300)
> ss = sapply(1:length(s),function(i){ll=lo;ll[5]=s[i];llnormfn(ll)})
> plot(s,ss)

给出：

R优化多个参数

R optimize multiple parameters

r

mathematical-optimization