优化函数中的起始值问题
problem with starting values in optim function
我正在使用 optim 函数处理优化问题。要最大化的函数是似然函数。我正在尝试评估一长串要评估的数据集,在某些情况下它会变得混乱,因为 lik.function 由于起始值而不会收敛。我提供的示例是函数未找到解决方案的示例。所以,我想知道一种方法,使最佳函数成为 select 网格中的起始值,找到解决方案,否则继续前进。这是我的代码,我尽量让它尽可能短。抱歉,我使用了很多 trycatch。
#data set
thisdata<-matrix(c(0.3014754, -1.8827312, 0.03221715, 0.08229814,
1.7730673, -0.9852836, 0.12997904, 0.04904762,
4.8520303, -1.2527630, 1.00781250, 0.12857143,
1.9582560, -3.0834379, 0.04961323, 0.17430025,
2.2284771, -2.5530445, 0.15824176, 0.08291110,
3.3672958, -1.6218604, 0.25862069, 0.07484568,
3.2358734, -1.3581235, 0.14847512, 0.06984127,
0.5930637, -3.3499041, 0.03696742, 0.51754386,
1.1451323, -3.0012725, 0.09415584, 0.11663597,
1.7147984, -3.3843903, 0.04370370, 0.17231638), nrow = 10, ncol=4, byrow = T)
colnames(thisdata)<-c('eta.obs', 'xi.obs', 'var.eta', 'var.xi')
#likelihood function
lik.to.optim <- function(theta, data){
mu.alpha <- theta[1]
beta <- theta[2]
mu.xi <- theta[3]
sigma2.xi <- theta[4]
sigma2.alpha<- theta[5]
if(sigma2.xi <= 0 | sigma2.alpha <=0 | beta^2*sigma2.xi-sigma2.alpha<0)
{
return(NA)
}
else{
Sigma<-matrix(c(beta^2*sigma2.xi-sigma2.alpha, beta*sigma2.xi-sigma2.alpha/beta,
beta*sigma2.xi-sigma2.alpha/beta, sigma2.xi), 2,2)
ris<-sum(dmvnorm(data[,1:2],c(mu.alpha+beta*mu.xi, mu.xi), Sigma, log=T))
}
return(ris)
}
#another function calling the previous lik. function
fun_adicional<-function(base){
NA.matrix<- matrix(NA, nrow=5, ncol=5)
unos<-c(1,1,1,1,1)
themle<-tryCatch(optim(unos, lik.to.optim, data=base, control=list(fnscale=-1))$par,
error=function(e) paste= c(NA,NA,NA,NA,NA),
warning=function(w) paste=c(NA,NA,NA,NA,NA))
hessiano0 <- tryCatch(optim(unos, lik.to.optim, data=base, control= list(fnscale=-1),
hessian=T)$hessian,
error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
lavar<-tryCatch(solve(-hessiano0), error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
se_actual <- sqrt(diag(lavar))
#double check
if(any(is.na(se_actual))){
#using the last mle estimate
init <- themle
ris.par <- tryCatch(optim(init, lik.to.optim, data=base, control=list(fnscale=-1))$par,
error=function(e) paste= c(NA,NA,NA,NA,NA),
warning=function(w) paste=c(NA,NA,NA,NA,NA))
hessiano0 <- tryCatch(optim(ris.par, lik.to.optim, data=base, control= list(fnscale=-1), hessian=T)$hessian,
error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
naive_var_a<-tryCatch(solve(-hessiano0), error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
V_CV_list<-naive_var_a
estimas<-ris.par
st_err<-sqrt(diag(V_CV_list))
} else{
estimas<-mle.par
V_CV_list<-lavar
st_err<-se_actual
}
here<-list(themle,lavar)
return(here)
}
fun_adicional(thisdata)
根据您的评论,我猜您想要 运行 网格搜索。我将使用 ?optim
:
中的一个简单示例
fr <- function(x) { ## Rosenbrock Banana function
x1 <- x[1]
x2 <- x[2]
100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
此函数的解是一个长度为 2 的向量 x
。
所以,运行一个gridSearch
.
library("NMOF")
gridSearch(function(x, fr) optim(x, fn = fr)$value,
fr = fr,
levels = list(c(-1, -2), ## the starting values
c(1, 2, 3)))
## 2 variables with 2, 3 levels: 6 function evaluations required.
## $minfun
## [1] 4.731118e-08
##
## $minlevels
## [1] -1 1
所以最好的起始值是 c(-1, 1)
,这导致
objective 函数值 4.731118e-08
.
另一个的objective函数值(values
)
起点 (levels
) 也被返回。
## $values
## [1] 4.731118e-08 1.605096e-06 4.568377e-07 1.006949e-06
## [5] 1.622874e-06 2.821143e-07
##
## $levels
## $levels[[1]]
## [1] -1 1
##
## $levels[[2]]
## [1] -2 1
##
## $levels[[3]]
## [1] -1 2
##
## $levels[[4]]
## [1] -2 2
##
## $levels[[5]]
## [1] -1 3
##
## $levels[[6]]
## [1] -2 3
我正在使用 optim 函数处理优化问题。要最大化的函数是似然函数。我正在尝试评估一长串要评估的数据集,在某些情况下它会变得混乱,因为 lik.function 由于起始值而不会收敛。我提供的示例是函数未找到解决方案的示例。所以,我想知道一种方法,使最佳函数成为 select 网格中的起始值,找到解决方案,否则继续前进。这是我的代码,我尽量让它尽可能短。抱歉,我使用了很多 trycatch。
#data set
thisdata<-matrix(c(0.3014754, -1.8827312, 0.03221715, 0.08229814,
1.7730673, -0.9852836, 0.12997904, 0.04904762,
4.8520303, -1.2527630, 1.00781250, 0.12857143,
1.9582560, -3.0834379, 0.04961323, 0.17430025,
2.2284771, -2.5530445, 0.15824176, 0.08291110,
3.3672958, -1.6218604, 0.25862069, 0.07484568,
3.2358734, -1.3581235, 0.14847512, 0.06984127,
0.5930637, -3.3499041, 0.03696742, 0.51754386,
1.1451323, -3.0012725, 0.09415584, 0.11663597,
1.7147984, -3.3843903, 0.04370370, 0.17231638), nrow = 10, ncol=4, byrow = T)
colnames(thisdata)<-c('eta.obs', 'xi.obs', 'var.eta', 'var.xi')
#likelihood function
lik.to.optim <- function(theta, data){
mu.alpha <- theta[1]
beta <- theta[2]
mu.xi <- theta[3]
sigma2.xi <- theta[4]
sigma2.alpha<- theta[5]
if(sigma2.xi <= 0 | sigma2.alpha <=0 | beta^2*sigma2.xi-sigma2.alpha<0)
{
return(NA)
}
else{
Sigma<-matrix(c(beta^2*sigma2.xi-sigma2.alpha, beta*sigma2.xi-sigma2.alpha/beta,
beta*sigma2.xi-sigma2.alpha/beta, sigma2.xi), 2,2)
ris<-sum(dmvnorm(data[,1:2],c(mu.alpha+beta*mu.xi, mu.xi), Sigma, log=T))
}
return(ris)
}
#another function calling the previous lik. function
fun_adicional<-function(base){
NA.matrix<- matrix(NA, nrow=5, ncol=5)
unos<-c(1,1,1,1,1)
themle<-tryCatch(optim(unos, lik.to.optim, data=base, control=list(fnscale=-1))$par,
error=function(e) paste= c(NA,NA,NA,NA,NA),
warning=function(w) paste=c(NA,NA,NA,NA,NA))
hessiano0 <- tryCatch(optim(unos, lik.to.optim, data=base, control= list(fnscale=-1),
hessian=T)$hessian,
error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
lavar<-tryCatch(solve(-hessiano0), error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
se_actual <- sqrt(diag(lavar))
#double check
if(any(is.na(se_actual))){
#using the last mle estimate
init <- themle
ris.par <- tryCatch(optim(init, lik.to.optim, data=base, control=list(fnscale=-1))$par,
error=function(e) paste= c(NA,NA,NA,NA,NA),
warning=function(w) paste=c(NA,NA,NA,NA,NA))
hessiano0 <- tryCatch(optim(ris.par, lik.to.optim, data=base, control= list(fnscale=-1), hessian=T)$hessian,
error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
naive_var_a<-tryCatch(solve(-hessiano0), error=function(e) paste= NA.matrix,
warning=function(w) paste=NA.matrix)
V_CV_list<-naive_var_a
estimas<-ris.par
st_err<-sqrt(diag(V_CV_list))
} else{
estimas<-mle.par
V_CV_list<-lavar
st_err<-se_actual
}
here<-list(themle,lavar)
return(here)
}
fun_adicional(thisdata)
根据您的评论,我猜您想要 运行 网格搜索。我将使用 ?optim
:
fr <- function(x) { ## Rosenbrock Banana function
x1 <- x[1]
x2 <- x[2]
100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
此函数的解是一个长度为 2 的向量 x
。
所以,运行一个gridSearch
.
library("NMOF")
gridSearch(function(x, fr) optim(x, fn = fr)$value,
fr = fr,
levels = list(c(-1, -2), ## the starting values
c(1, 2, 3)))
## 2 variables with 2, 3 levels: 6 function evaluations required.
## $minfun
## [1] 4.731118e-08
##
## $minlevels
## [1] -1 1
所以最好的起始值是 c(-1, 1)
,这导致
objective 函数值 4.731118e-08
.
另一个的objective函数值(values
)
起点 (levels
) 也被返回。
## $values
## [1] 4.731118e-08 1.605096e-06 4.568377e-07 1.006949e-06
## [5] 1.622874e-06 2.821143e-07
##
## $levels
## $levels[[1]]
## [1] -1 1
##
## $levels[[2]]
## [1] -2 1
##
## $levels[[3]]
## [1] -1 2
##
## $levels[[4]]
## [1] -2 2
##
## $levels[[5]]
## [1] -1 3
##
## $levels[[6]]
## [1] -2 3