problem with starting values in optim function

I'm working on an optimization problem using the optim function. The function to maximize is a likelihood function. I'm trying a very long list of data sets to be evaluated, and in some cases It gets messy because the lik.function does not converge because of the starting values. The example I provide is one in which the function does not find a solution. So, I would like to know a way to make the optimal function to select the starting values out of a grid of them, to find a solution otherwise to move on. This is my code, I try to make it the shortest I could. I'm sorry I'm using a lot of 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)

Solution

From your comment I guess you want to run a grid search. I'll use a simple example from ?optim:

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

The solution for this function is a vector x of length 2. So, run a 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

So the best starting values are c(-1, 1), which lead to an objective function value of 4.731118e-08.

The objective function values (values) for the other starting points (levels) are also returned.

## $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