r 中的优化：非有限有限差分误差答案

【问题标题】：optim in r :non finite finite difference errorr 中的优化：非有限有限差分误差
【发布时间】：2016-03-25 12:41:35
【问题描述】：

我有一个想要最大化的简单似然函数（来自均值 = 0 的正常分布）。 optim一直给我这个错误：优化错误（par = phi，fn = loglike，估计 = 估计，NULL，hessian = TRUE，：非有限有限差分值 [1]

这是我的数据和似然函数：

y =  [ -0.01472  0.03942  0.03592 0.02776 -0.00090 ]

C = a varcov matrix:

  1.66e-03 -0.000120 -6.780e-06  0.000102 -4.000e-05

 -1.20e-04  0.001387  7.900e-05 -0.000140 -8.000e-05

 -6.78e-06  0.000079  1.416e-03 -0.000070  8.761e-06

  1.02e-04 -0.000140 -7.000e-05  0.001339 -6.000e-05

 -4.00e-05 -0.000080  8.761e-06 -0.000060  1.291e-03

我的对数似然函数是： lglkl = -.5*(log(det(v)) + (t(y)%%vi%%y))` 其中 v = phi*I + C 和 vi=inverse(v ) 和 I= 5*5 单位矩阵。

我正在尝试获取“phi”的 mle 估计值。我认为这将是一个简单的优化问题，但我正在苦苦挣扎。非常感谢任何帮助。提前致谢。我的代码如下：

loglike <- function(phi,y) {

v = phi*I + C
vi = solve(v)
loglike = -.5*(log(det(v)) + (t(y)%*%vi%*%y))
return(-loglike)
}

phi = 0
parm <- optim(par=phi,fn=loglike,y=y,NULL,hessian = TRUE, method="L-BFGS-B",lower=0,upper=1000)

【问题讨论】：

如果您使用有效的 R 语法共享您的数据（y 和 V），那就太棒了。使用dput() 非常适合。
你是这个意思吗？ dput(y) c(-0.01472, 0.03942, 0.03592, 0.02776, -9e-04) > dput(v) 结构(c(0.00166, -0.00012, -6.78e-06, 0.000102, -4e-05, -0.00012, 0.001387、7.9e-05、-0.00014、-8e-05、-6.78e-06、7.9e-05、0.001416、-7e-05、8.761e-06、0.000102、-0.00014、-7e-05、0.001339、 -6e-05, -4e-05, -8e-05, 8.761e-06, -6e-05, 0.001291), .Dim = c(5L, 5L ))
没错。我可以将其复制/粘贴到 R 中以获取 y 和 V 对象。理想情况下，您会将其编辑到您的问题中。
刚刚对我的帖子进行了更正。我之前的帖子中几乎没有错误。

标签： r optimization

【解决方案1】：

您遇到的错误是因为 φ 在一定数量的迭代之后变为负数（这表明算法没有正确应用约束）。此外，解决方案不会收敛到单个值，而是会在几个小值之间跳转，然后才能达到更新后的协方差矩阵不再是正定的情况。在那个阶段你得到 det(v) optim 算法在那个阶段退出。

要查看发生了什么，请在下面的代码中使用maxit 和ndeps 参数。

require("matrixcalc")

#-------------------------------------------------
# Log-likelihood function
#-------------------------------------------------
loglike <- function(phi, y) {

  # Shift the covariance matrix
  print(paste("phi = ", phi))
  #v = phi*I + (1 - phi)*C
  v = phi*I + C
  stopifnot(is.positive.definite(v))

  # Invert shifted matrix
  vi = solve(v)

  # Compute log likelihood
  loglike = -.5*(log(det(v)) + (t(y) %*% vi %*% y))
  print(paste("L = ", loglike))

  return(-loglike)
}

#-------------------------------------------------
# Data
#-------------------------------------------------
y = c(-0.01472, 0.03942, 0.03592, 0.02776, -9e-04) 
C =  structure(c(0.00166, -0.00012, -6.78e-06, 0.000102, -4e-05, -0.00012, 
                 0.001387, 7.9e-05, -0.00014, -8e-05, -6.78e-06, 7.9e-05, 
                 0.001416, -7e-05, 8.761e-06, 0.000102, -0.00014, -7e-05, 
                 0.001339, -6e-05, -4e-05, -8e-05, 8.761e-06, -6e-05, 0.001291), 
                 .Dim = c(5L, 5L ))

#--------
# Initial parameter
#--------
I = diag(5)
phi = 50

#--------
# Minimize
#--------
parm <- optim(par = phi, fn = loglike, y = y, NULL, hessian = TRUE, 
              method = "L-BFGS-B", lower = 0.0001, upper = 1000,
              control = list(trace = 3,
                             maxit = 1000,
                             ndeps = 1e-4) )

【讨论】：