我的问题（以及我认为可以帮助解决它的问题）在“FOR REPRODUCTION”行之前得到了解释。之后我刚刚发布了我的代码，以防万一复制可能有助于解决它。
我使用 optim 和 constrOptim.nl 来解决我编写的函数 g （见下文）内的约束优化问题。
我知道下面使用的初始值并不理想，但我选择它们是因为它们会导致我在较短的程序中遇到问题。我使用该程序将模型参数校准为数据，并且对于更好的初始值、更高的公差等也会出现此问题。
错误
我调用我编写的函数 get_par ：
v<-c(0.12504710,0.09329359,0.06778733, 0.04883216, 0.04187344,0.02886261,0.02332951,0.02178576,0.02282214,0.02956336,0.03478598)
Ti=1/12
x<-log(cbind(0.8,0.85,0.9,0.95,0.975,1,1.025,1.05,1.1,1.15,1.2))
g(par2=c(-5,5),v=v,Ti=Ti,x=x)
然后我得到
Error in optim: inital value 'vmmin' is not finite.
到目前为止我所观察到的
所以我开始调试我的代码来找出这个错误到底发生在哪里。错误发生在函数 g（见下文）中（值 sigma=5,m=-5,y=(xm)/sigma,vtilde=v/12）
#print(paste("vW: sigma: ",sigma,"mv:",mv))
argmin<-constrOptim.nl(par=c(3*sigma,sigma,mv/2),fn=f,hin.jac=hinv.jac,
hin=hinv,heq.jac=heqv.jac,heq=heqv,control.outerlist(trace=T),
control.optim=list(abstol=10^(-10)),y=y,vtilde=vtilde,sigma=sigma)
函数 constrOptim.nl 的轨迹显示
Outer iteration: 18
Min(hin): 1.026858e-19 Max(abs(heq)): 0
par: 10 9.99998 1.02686e-19
fval = 6399
最后一次迭代。我猜想在最后一次迭代中出现 1.02686e-19 存在某种数值问题。
我查看了函数 constrOptim.nl 和 albama （使用 debug() ），并且错误恰好发生在该行中
theta.old <- theta
atemp <- optim(par = theta, fn = fun, gr = grad, control = control.optim,
method = "BFGS", hessian = TRUE, ...)
其中 theta=theta.old 具有值
Browse[2]> theta.old
[1] 1.000002e+01 9.999985e+00 -3.349452e-20
因此它有一个略低于零的条目（它的绝对值甚至小于机器精度，不是吗？）。
当您查看函数 fun 时，您会意识到它调用了函数
R:
function (theta, theta.old, ...)
{
gi <- hin(theta, ...)
if (any(gi < 0))
return(NaN)
gi.old <- hin(theta.old, ...)
hjac <- hin.jac(theta.old, ...)
bar <- sum(gi.old * log(gi) - hjac %*% theta)
if (!is.finite(bar))
bar <- -Inf
fn(theta, ...) - mu * bar
}
hin(theta,...)=hinv(theta,...) 返回一个带有负数的向量，因此该函数返回 NaN。我想这应该会导致错误消息：“优化错误：初始值'vmmin'不是有限的”。我现在的问题是：
我该如何解决？我想过在出现如此小的值时强制程序以某种方式终止，但我还没有设法做到这一点。你有什么建议？
提前谢谢了，
繁殖：
这是我的程序：
函数 hinv、hinv.jac、heq 和 heq.jac 仅用于约束。我优化的函数是g。
library(alabama)
library(dfoptim)
#function f, par = (c,d,atilde)
f<-function(par3,y,vtilde,sigma){
sum((par3[3]+par3[2]*y+par3[1]*sqrt(y^2+1)-vtilde)^2)
}
#Equality/Inequality constraints
heqv<-function(par3,y,vtilde,sigma){
J1<-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2)
J2<-matrix(0,nrow=3,ncol=3)
J2[1:2,1:2]<-J1
J2[3,3]<-1
j<-J2%*%par3
j[2]-2*sqrt(2)*sigma
}
#Jacobian-matrix
hinv.jac<-function(par3,y,vtilde,sigma){
#J1, J2: Drehungen für die constraints
J1<-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2)
J2<-matrix(0,nrow=3,ncol=3)
J2[1:2,1:2]<-J1
J2[3,3]<-1
hjac<-matrix(cbind(1,-1,0,0,0,0,0,0,0,0,1,-1),nrow=4)%*%J2
hjac
}
hinv<-function(par3,y,vtilde,sigma){
#J1, J2: Drehungen für die constraints
J1<-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2)
J2<-matrix(0,nrow=3,ncol=3)
J2[1:2,1:2]<-J1
J2[3,3]<-1
j<-J2%*%par3
h<-rep(NA,4)
h[1]<- j[1]
h[2]<- sqrt(2)*2*sigma-j[1]
h[3]<-j[3]
h[4]<-max(vtilde)-j[3]
h
}
#Jacobian-matrix
heqv.jac<-function(par3,y,vtilde,sigma){
#J1, J2: Drehungen für die constraints
J1<-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2)
J2<-matrix(0,nrow=3,ncol=3)
J2[1:2,1:2]<-J1
J2[3,3]<-1
cbind(J2[2,1],J2[2,2],0)
}
#function g input: par2= (m,sigma): optimization of function f
g<-function(par2,v,Ti,x){
#definition of parameters being used
m<-par2[1]
sigma<-par2[2]
y<-(x-m)/sigma #Transformation von x zu y gemäß paper
vtilde<-Ti*v
mv<-max(vtilde)
#print(paste("vW: sigma: ",sigma,"mv:",mv))
argmin<-constrOptim.nl(par=c(3*sigma,sigma,mv/2),fn=f,hin.jac=hinv.jac,hin=hinv,heq.jac=heqv.jac,heq=heqv,control.outer=list(trace=F),control.optim=list(abstol=10^(-10)),y=y,vtilde=vtilde,sigma=sigma)
argmin$par
}