My problem (and what I think can help solve it) is explained before the line “FOR PLAYBACK”. After that, I just posted my code, just in case playback can help solve it.
I use optim and constrOptim.nl to solve the optimization problem with constraints inside the g function (see below).
I know that the initial values ​​used below are not ideal, but I chose them because they cause a problem that I encounter in a shorter program. I use this program to calibrate model parameters for data, and this problem also occurs for better values, higher tolerances, etc.
Error
I call the get_par function, which I wrote with:
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)
Then i get
Error in optim: inital value 'vmmin' is not finite.
What I watched so far
So, I started debugging my code to find out exactly where this error occurs. The error occurs in the function g (see below) in the line (with values ​​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)
Function tracking constrOptim.nl displays
Outer iteration: 18 Min(hin): 1.026858e-19 Max(abs(heq)): 0 par: 10 9.99998 1.02686e-19 fval = 6399
for the last iteration. I assume that there is some numerical problem with 1.02686e-19 appearing on the last iteration.
I looked at the constrOptim.nl and albama functions (using debug ()), and the error exactly occurs on the line
theta.old <- theta atemp <- optim(par = theta, fn = fun, gr = grad, control = control.optim, method = "BFGS", hessian = TRUE, ...)
where theta = theta.old matters
Browse[2]> theta.old [1] 1.000002e+01 9.999985e+00 -3.349452e-20
Therefore, it has a record that is slightly below zero (its absolute value is even less than the accuracy of the machine, right?).
When you look at the fun function, you realize that it calls the function
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, ...) returns a negative-written vector, so this function returns NaN. I believe this should cause an error message: "Error in optimization: the value inital" vmmin "is not finite." My question is:
How can i fix this? I was thinking about making the program somehow stop when such small values ​​happen, but I haven't done it yet. What are you offering?
Thank you very much in advance,
FOR PLAYBACK:
Here is my program:
The hinv, hinv.jac, heq, and heq.jac functions are for limitations only. The function in which I am optimized is 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<-(xm)/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 }