Ошибка вsolve.default(oout$hessian): подпрограмма Lapack dgesv: система точно в единственном числе: U[1,1] = 0

Я использую метод максимального правдоподобия для оценки набора параметров. Теперь я собираюсь использовать функцию mle из пакета stats4 в R, чтобы сделать вероятность профиля для одного из параметров. Для этого мне нужно исправить один из параметров при вызове функции mle. Вот код:

fr <- function(x1, x2, x3) {   
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2 + x3
}

out <- mle(fr,start = list(x1=1, x2=2, x3=3), method="Nelder-Mead",
           control=list(trace=4), fixed = list(x2=1))

и я получаю эту ошибку:

Ошибка вsolve.default(oout$hessian): подпрограмма Lapack dgesv: система точно в единственном числе: U[1,1] = 0

Если я не использую вариант fixed, то этой ошибки у меня нет, но результат не профильный правдоподобие. Не могли бы вы сообщить мне, как я могу решить эту проблему?


r optimization statistics mathematical-optimization log-likelihood
person Maryam    schedule 02.05.2021    source источник
comment
Нам действительно нужен минимальный воспроизводимый пример, чтобы справиться с этим...   -  person Ben Bolker    schedule 03.05.2021
comment
@BenBolker Я отредактировал свой пост и привел пример   -  person Maryam    schedule 03.05.2021

Ответы (1)


arrow_upward
1
arrow_downward

tl;dr Я не уверен, что ваша целевая функция имеет смысл, я предполагаю, что у вас опечатка. (Более того, если ваша целевая функция работает с mle, вам не нужно явно задавать fixed: метод profile автоматически вычислит для вас профили правдоподобия...)

Давайте начнем с полной модели и воспользуемся optim(), а не stats4::mle() (я знаю, что вы хотите вернуться к mle, чтобы можно было выполнить профилирование правдоподобия, но отлаживать optim() проблемы немного проще, так как на один слой кода меньше перекопать.)

Поскольку optim() нужна целевая функция, которая принимает вектор, а не список аргументов, напишите оболочку (мы также могли бы использовать do.call(fr, as.list(p))):

fr0 <- function(p) {
  fr(x1=p[1], x2=p[2], x3=p[3])
}
opt1 <- optim(fn=fr0, par=c(1,2,3), method="Nelder-Mead")

Полученные результаты:

$par
[1]    28.51486   812.09978 -7095.39630

$value
[1] -6238.881

$counts
function gradient 
     502       NA 

$convergence
[1] 1

Обратите внимание, что значение x[3] сильно отрицательное, как и значение целевой функции, а код сходимости не равен нулю: в частности, это означает (из ?optim):

«1» указывает, что предел итераций «maxit» был достигнут.

Если мы установим control=list(maxit=2000) и попробуем еще раз x3, целевая функция станет еще меньше, а код сходимости по-прежнему равен 1!

Затем мы более внимательно смотрим на целевую функцию и замечаем, что она достигает -Inf при x3 → Inf, так что мы никогда не получим ответ. (Предположительно, в какой-то момент мы доберемся до задачи с плавающей запятой, но 10 миллионов итераций дают нам только -1e17...)

Если я изменю x3 на x3^2 в вашей функции, все будет работать нормально ... может быть, это то, что вы хотели ... ???

person Ben Bolker    schedule 07.05.2021