Почему оптимизация дмультином MLE ложно сходится в R? - PullRequest
Новая
       15

Почему оптимизация дмультином MLE ложно сходится в R?

1 голос
/ 19 июня 2019

Я пишу MLE-функцию Линкольна-Петерсена в R. У меня есть две версии моей функции распределения, одна, которая использует гамма-функцию, и другая, которая использует Rs с внутренним dmultinom (который использует lgamma). Оба возвращают идентичные результаты, когда я вставляю значения выборки, однако версия dmultinom не может обеспечить приемлемые оценки MLE. Мне любопытно, что вызывает это. 

llik_lincoln_multinom1 <- function(par,n_1,n_2,m_2) {
  N   <- par[1]
  p_1 <- par[2]
  p_2 <- par[3]
  q_1 <- 1-p_1
  q_2 <- 1-p_2
  r   <- n_1 + n_2 - m_2
  l   <- (gamma(sum(c(m_2,n_1-m_2,n_2-m_2,N-r))+1)/prod(gamma(c(m_2,n_1-m_2,n_2-m_2,N-r)+1)))*(((p_1*p_2)^m_2)*((p_1*q_2)^(n_1-m_2))*((q_1*p_2)^(n_2-m_2))*((q_1*q_2)^(N-r)))
  return(-log(l))
}
optim(par=c(20,0.1,0.1),fn=llik_lincoln_multinom1,
      n_1=10,n_2=10,m_2=2,
      lower=c(18,0,0),
      upper=c(Inf,1,1),
      method="L-BFGS-B",
      control=list(parscale=c(100,1,1)))

llik_lincoln_multinom2 <- function(par,n_1,n_2,m_2) {
  N   <- par[1]
  p_1 <- par[2]
  p_2 <- par[3]
  q_1 <- 1-p_1
  q_2 <- 1-p_2
  r   <- n_1 + n_2 - m_2
  l   <- dmultinom(c(m_2,n_1-m_2,n_2-m_2,N-r),prob=c(p_1*p_2,p_1*q_2,q_1*p_2,q_1*q_2))
  return(-log(l))
}
optim(par=c(20,0.1,0.1),fn=llik_lincoln_multinom2,
      n_1=10,n_2=10,m_2=2,
      lower=c(18,0,0),
      upper=c(Inf,1,1),
      method="L-BFGS-B",
      control=list(parscale=c(100,1,1)))

Во второй функции градиент для параметра N представляется равным 0, и MLE никогда не удаляется от начального параметра. Почему градиент для второй версии отличается от первого?

1 Ответ

0 голосов
/ 20 июня 2019

dmultinom(c(m_2,n_1-m_2,n_2-m_2,N-r),prob=......) совпадает с dmultinom(c(m_2,n_1-m_2,n_2-m_2,as.integer(N-r)),prob=......).Таким образом, небольшое изменение N не меняет это значение, потому что as.integer(N-r) не изменяется: 

> llik_lincoln_multinom2(c(20, 0.5, 0.5), 10, 10, 2)
[1] 7.985771
> llik_lincoln_multinom2(c(20.1, 0.5, 0.5), 10, 10, 2)
[1] 7.985771
> llik_lincoln_multinom2(c(20.2, 0.5, 0.5), 10, 10, 2)
[1] 7.985771

Так что optim «думает», что первый параметр, N, не оказывает влиянияна значение функции.

...