R: optim和nlme错误
R: optim和nlme错误
提问于 2017-08-09 15:59:50
library(nlme)
mydat <-
structure(list(class = c(91L, 91L, 91L, 91L, 91L, 91L, 92L, 92L,
92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L,
92L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L,
94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L,
95L, 95L, 95L), days = c(5.7, 11.1, 13.9, 15.3, 18.3, 18.9, 1.9,
2.1, 2.9, 3.4, 4.4, 5, 6.9, 10.4, 11.6, 13, 13.4, 15.7, 15.9,
17.3, 17.7, 19.4, 2.3, 12.6, 15.4, 2, 2.4, 4.9, 5.4, 5.7, 7,
7.1, 7.7, 9.1, 9.6, 12.6, 16.6, 17.1, 1, 2, 4.4, 5.6, 5.7, 10.4,
12.1, 12.9, 13, 15.6, 16.1, 18.6), outcome = c(3.31586988521325,
3.26226473964254, 3.26098046236747, 3.26086828734987, 3.26081582971007,
3.2608136610639, 1.54175217273846, 1.74336277564818, 2.48010804039677,
2.73455940271066, 2.86602619542132, 2.85753365781511, 2.81667209739959,
2.80430238913247, 2.80395479988755, 2.80383291979961, 2.8038189189449,
2.80379174103878, 2.80379113213262, 2.80378896928776, 2.80378871890839,
2.80378827755366, 1.96537220180574, 3.00124636046136, 3.00096700482166,
2.05608815148142, 2.44248026102198, 4.03918455971327, 4.08570704450138,
4.09781416453829, 4.09869791687544, 4.09759058744364, 4.09084045815843,
4.07921200433542, 4.07668896637335, 4.07081047825795, 4.06993724272757,
4.06991925387706, 1.07225026715462, 3.72090308875724, 4.93988448353623,
4.7209971449984, 4.70681751285687, 4.49435510591282, 4.48824648431355,
4.4870870783591, 4.48698191392076, 4.48574152736339, 4.48566703246688,
4.48551396890595), s = c(0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693)), .Names = c("class",
"days", "outcome", "s"), row.names = 1001:1050, class = "data.frame")
library(nlme)
obj_NM <- function(arg){
model = nlme(outcome ~ exp(beta1) * s * days / (s * days + exp(1 - beta3 * s * days)) +
exp(beta2) * exp(1 - beta4 * s * days) / (s * days + exp(1 - beta4 * s * days)),
data = mydat,
fixed = list(beta1 ~ 1, beta2 ~ 1, beta3 ~ 1, beta4 ~ 1),
random = list(class = pdDiag(list(beta1 ~ 1, beta2 ~ 1, beta4 ~ 1))),
start = list(fixed = c(beta1 = arg[1], beta2 = arg[2], beta3 = arg[3], beta4 = arg[4])), verbose = FALSE)
return(model$logLik)
}
control = list(fnscale = -1)
optim(par = c(1.310239, -4.668217, 17.01345, 2.402943), fn = obj_NM, hessian = TRUE, control = control)运行上述代码将给出错误和警告：
Error in chol.default((value + t(value))/2) :
the leading minor of order 2 is not positive definite In addition: Warning messages:
1: In nlme.formula(outcome ~ exp(beta1) * s * days/(s * days + exp(1 - :
Singular precision matrix in level -1, block 1
2: In nlme.formula(outcome ~ exp(beta1) * s * days/(s * days + exp(1 - :
Singular precision matrix in level -1, block 1我在这里的目标是获得恒河矩阵，并研究为什么我的nlme模型可能不适合。我试图最大化我的目标函数，因此我设置了fnscale = -1 (文档说，为了使optim执行最大化，它应该是负的)。但是，我不知道如何处理错误信息。有没有一种方法可以让optim输出恒河矩阵？nlme的错误似乎阻止了它这样做。
回答 1
Stack Overflow用户
发布于 2022-10-12 01:30:38
我已经能够使用DEoptim进行初始的“预优化步骤”来优化这个函数：
mydat <- structure(list(class = c(91L, 91L, 91L, 91L, 91L, 91L, 92L, 92L,
92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L,
92L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L,
94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L,
95L, 95L, 95L),
days = c(5.7, 11.1, 13.9, 15.3, 18.3, 18.9, 1.9,
2.1, 2.9, 3.4, 4.4, 5, 6.9, 10.4, 11.6, 13, 13.4, 15.7, 15.9,
17.3, 17.7, 19.4, 2.3, 12.6, 15.4, 2, 2.4, 4.9, 5.4, 5.7, 7,
7.1, 7.7, 9.1, 9.6, 12.6, 16.6, 17.1, 1, 2, 4.4, 5.6, 5.7, 10.4,
12.1, 12.9, 13, 15.6, 16.1, 18.6),
outcome = c(3.31586988521325, 3.26226473964254, 3.26098046236747, 3.26086828734987, 3.26081582971007,
3.2608136610639, 1.54175217273846, 1.74336277564818, 2.48010804039677,
2.73455940271066, 2.86602619542132, 2.85753365781511, 2.81667209739959,
2.80430238913247, 2.80395479988755, 2.80383291979961, 2.8038189189449,
2.80379174103878, 2.80379113213262, 2.80378896928776, 2.80378871890839,
2.80378827755366, 1.96537220180574, 3.00124636046136, 3.00096700482166,
2.05608815148142, 2.44248026102198, 4.03918455971327, 4.08570704450138,
4.09781416453829, 4.09869791687544, 4.09759058744364, 4.09084045815843,
4.07921200433542, 4.07668896637335, 4.07081047825795, 4.06993724272757,
4.06991925387706, 1.07225026715462, 3.72090308875724, 4.93988448353623,
4.7209971449984, 4.70681751285687, 4.49435510591282, 4.48824648431355,
4.4870870783591, 4.48698191392076, 4.48574152736339, 4.48566703246688,
4.48551396890595),
s = c(0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693)),
.Names = c("class", "days", "outcome", "s"),
row.names = 1001 : 1050,
class = "data.frame")
library(nlme)
library(DEoptim)
obj_NM <- function(arg, bool_Print = FALSE)
{
logLik <- tryCatch(nlme(outcome ~ exp(beta1) * s * days / (s * days + exp(1 - beta3 * s * days)) +
exp(beta2) * exp(1 - beta4 * s * days) / (s * days + exp(1 - beta4 * s * days)), data = mydat,
fixed = list(beta1 ~ 1, beta2 ~ 1, beta3 ~ 1, beta4 ~ 1),
random = list(class = pdDiag(list(beta1 ~ 1, beta2 ~ 1, beta4 ~ 1))),
start = list(fixed = c(beta1 = arg[1], beta2 = arg[2], beta3 = arg[3], beta4 = arg[4])), verbose = FALSE)$logLik,
error = function(e) NA)
if(bool_Print == TRUE)
{
print(logLik)
}
if(is.na(logLik))
{
return(10 ^ 30)
}else
{
return(-logLik)
}
}
obj_Iter <- DEoptim(fn = obj_NM, lower = rep(-1000, 4), upper = rep(1000, 4))
optim(par = obj_Iter$optim$bestmem, fn = obj_NM, hessian = TRUE, method = "BFGS")
$par
par1 par2 par3 par4
184.59178 178.77184 179.57188 90.01865
$value
[1] 2.010262
$counts
function gradient
68 1
$convergence
[1] 0
$message
NULL页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/45595549
复制相关文章
相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号