0

我正在使用 optim 函数解决优化问题。最大化的函数是似然函数。我正在尝试一个很长的要评估的数据集列表，在某些情况下它会变得混乱，因为 lik.function 由于起始值而不会收敛。我提供的示例是函数找不到解决方案的示例。所以，我想知道一种方法来制作最佳函数以从它们的网格中选择起始值，以找到解决方案，否则继续前进。这是我的代码，我尽量让它最短。对不起，我使用了很多 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)
4

1 回答 1

1

从您的评论中，我猜您想运行网格搜索。我将使用一个简单的示例?optim：

fr <- function(x) {   ## Rosenbrock Banana function
    x1 <- x[1]
    x2 <- x[2]
    100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}

这个函数的解是一个x长度为 2 的向量。所以，运行一个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)，这导致目标函数值为4.731118e-08。

还返回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
于 2020-06-12T14:50:18.400 回答