为什么 dmultinom MLE 的优化在 R 中错误地收敛？

Question

我正在用 R 编写 Lincoln-Petersen MLE 函数。我有两个版本的分布函数，一个使用 gamma 函数，另一个使用 Rs 内部 dmultinom（使用 lgamma）。当我插入样本值时，两个 return 相同的结果，但是 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 永远不会离开起始参数。为什么第二个版本的渐变与第一个不同？

Answer 1

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 "thinks" 第一个参数 N 对函数的值没有影响。

为什么 dmultinom MLE 的优化在 R 中错误地收敛？

Why is optimization of dmultinom MLE falsely converging in R?

optimization

r

probability