проблема с начальными значениями в функции оптимизации

Я работаю над проблемой оптимизации, используя функцию 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)

r optimization log-likelihood
person E. Pesantez    schedule 12.06.2020    source источник
comment
Не могли бы вы немного уточнить свой вопрос. Мне это не совсем ясно. У вас есть проблемы с поиском алгоритма, который создает хорошие отправные точки (т.е. хорошо распределяется по пространству параметров)? Или вы надеялись, что функция optim справится с этим сама? В этом случае я не думаю, что это возможно. Вам придется написать свой собственный цикл: 1. выбрать начальную точку, 2. оптимизировать, 3. повторить попытку, если решение неприемлемо.   -  person Jan    schedule 12.06.2020
comment
Эй, я не хочу создавать алгоритм, который создает хорошие отправные точки. Допустим, я придумываю возможные начальные значения, которые я хотел бы попробовать, скажем, их сетку. Я хочу, чтобы оптим перепробовал все комбинации сетки.   -  person E. Pesantez    schedule 12.06.2020

Ответы (1)


arrow_upward
0
arrow_downward

Судя по вашему комментарию, вы хотите запустить поиск по сетке. Я буду использовать простой пример из ?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
person Enrico Schumann    schedule 12.06.2020