8

我的问题（以及我认为可以帮助解决它的问题）在“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
}
4

0 回答 0