I'm writing a Lincoln-Petersen MLE function in R. I have two versions of my distribution function, one that uses the gamma function, and one that uses Rs internal dmultinom (which uses lgamma). Both return identical results when I plug in sample values, however the dmultinom version fails to provide acceptable MLE estimates. I'm curious as to what is causing this.
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)))
In the second function, the gradient for the N parameter appears to be 0 and the MLE never moves away from the starting parameter. Why would the gradient for the second version be different than the first?
dmultinom(c(m_2,n_1-m_2,n_2-m_2,N-r),prob=......) is the same as dmultinom(c(m_2,n_1-m_2,n_2-m_2,as.integer(N-r)),prob=......). So a small variation of N does not change this value, because as.integer(N-r) does not change:
> 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
So optim "thinks" that the first parameter, N, has no impact on the value of the function.