私の問題（およびそれを解決するのに役立つと思われるもの）は、「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
}