I am using both nlminb and optim to get MLEs from a likelihood function I have
developed. AFAIK, the model I has not been previously used in this way and so I
am struggling a bit to unit test my code since I don't have another data set to
compare this kind of estimation to.
The likelihood I have is (in tex below)
\begin{equation}
\label{eqn:marginal}
L(\beta) = \prod_{s=1}^N \int \prod_{i=1}^K\frac{e^{x_{is}(\theta_s-\beta_i)}}
{x_{is}!e^{e^(\theta_s-\beta_i)}} f(\theta)d(\theta)
\end{equation}
Where I view $\theta$ as a nuisance parameter and so I integrate it out of the
likelihood. The goal is to get parameter estimates for $\beta$. The integral
cannot be easily evaluated so I approximate it as:
\begin{equation}
\label{eqn:marginal2}
L(\beta) \approx \prod_{s=1}^N \sum_{q=1}^Q
\prod_{i=1}^K\frac{e^{x_{is}(\theta_{q}-\beta_i)}}
{x_{is}!e^{e^(\theta_{q}-\beta_i)}} w_q
\end{equation}
\noindent where $\theta_{q}$ is the node at quadrature point $q = 1, \ldots, Q$
and $w_q$ is the weight at quadrature point $q$.
For now, I am assuming $f(\theta)$ is Uniform but this may change and that is
not a major issue. Now, I have written a function using both nlminb and optim.
I have copied my function below where the arguments are the data set, Q for the
number of quadrature points, and an option for providing starting values. I am
running into some flow control problems and so I hack it a bit and fix a lower
bound as you may see in the function.
Here is the issue at hand. First, nlminb and optim give different parameter
estimates. Optim tells me it converged but nlminb tells me I get relative
convergence. If I use different number of quadrature points in optim, I get
very different parameter estimates and this should not happen. But, varying the
number of quadrature points in nlminb yields the same parameter estimates, but
I am not sure the model converged.
I have also provided the data set I am working with below as well as
sessionInfo(). There may be many issues going on here and am grateful for any
reactions you all may have. I *think* my likelihood is written properly and I
*think* I am handling flow control issues in a relatively decent way. I may be
misinterpreting optim or nlminb reports.
In any event, should anyone take the time to review this and offers pointers I
would be grateful.
Harold
fit <- function(data, Q, startVal = NULL, ...){
if(!is.null(startVal) && length(startVal) != ncol(data) ){
stop("Length of argument startVal not equal to
the number of parameters estimated")
}
if(!is.null(startVal)){
startVal <- startVal
} else {
startVal <- rep(0, ncol(data))
}
#qq <- gauss.quad.prob(Q, dist = 'normal', mu = 0, sigma=1)
qq <- gauss.quad.prob(Q, dist = 'uniform', alpha = -5, beta=5)
rr1 <- matrix(0, nrow = Q, ncol = nrow(data))
data <- as.matrix(data)
L <- nrow(data)
C <- ncol(data)
fn <- function(b){
b <- b[1:C]
for(j in 1:Q){
for(i in 1:L){
rr1[j,i] <-
prod((exp(data[i,]*(qq$nodes[j]-b))) / (factorial(data[i,]) *
exp(exp(qq$nodes[j]-b)))) *
qq$weights[j]
}
}
rr1[rr1==0] <- 1e-10
-sum(log(colSums(rr1)))
}
#opt <- optim(startVal, fn, method = "BFGS", hessian = FALSE)
opt <- nlminb(startVal, fn)
#list("coefficients" = opt$par, "LogLik" = -opt$value,
"Std.Error" = sqrt(diag(solve(opt$hessian))))
#list("coefficients" = opt$par, "LogLik" = -opt$value)
opt
}
fit(data, Q = 10)
> data
cindyR cjR ohsR sdhpR
[1,] 24 14 6 9
[2,] 19 14 6 10
[3,] 19 14 6 10
[4,] 17 18 6 8
[5,] 19 17 5 8
[6,] 17 17 6 9
[7,] 13 15 5 11
[8,] 19 13 6 8
[9,] 21 10 4 11
[10,] 16 15 5 9
[11,] 14 16 5 10
[12,] 14 13 5 10
[13,] 17 15 5 8
[14,] 16 18 4 8
[15,] 16 14 5 8
[16,] 18 13 5 8
[17,] 15 14 6 7
[18,] 16 17 5 7
[19,] 18 12 5 8
[20,] 18 9 4 9
[21,] 17 9 4 10
[22,] 18 12 4 8
[23,] 15 15 5 7
[24,] 16 9 5 8
[25,] 18 11 5 7
[26,] 18 10 4 9
[27,] 15 15 4 7
[28,] 16 11 4 9
[29,] 11 16 5 8
[30,] 11 16 5 8
[31,] 16 13 4 8
[32,] 13 15 5 7
[33,] 17 17 2 9
[34,] 18 8 3 9
[35,] 13 14 4 8
[36,] 16 11 4 8
[37,] 16 11 4 8
[38,] 15 8 5 9
[39,] 14 15 4 8
[40,] 15 14 4 7
[41,] 14 11 4 8
[42,] 13 12 4 9
[43,] 15 13 3 8
[44,] 12 17 3 8
[45,] 20 0 4 8
[46,] 14 10 4 9
[47,] 11 15 5 7
[48,] 14 12 4 8
[49,] 11 14 4 8
[50,] 13 11 5 7
[51,] 12 13 5 6
[52,] 16 7 5 7
[53,] 10 13 4 8
[54,] 15 12 3 7
[55,] 16 7 4 8
[56,] 16 11 3 7
[57,] 12 15 4 7
[58,] 19 2 4 8
[59,] 13 13 4 7
[60,] 15 12 3 7
[61,] 14 15 3 7
[62,] 19 10 2 8
[63,] 13 7 5 7
[64,] 16 7 4 7
[65,] 16 2 5 7
[66,] 12 17 3 7
[67,] 13 6 4 9
[68,] 15 7 4 7
[69,] 13 11 3 8
[70,] 15 7 4 7
[71,] 12 14 4 6
[72,] 14 6 4 8
[73,] 15 6 4 7
[74,] 13 11 4 7
[75,] 15 8 4 7
[76,] 13 1 6 7
[77,] 12 13 3 7
[78,] 13 6 4 8
[79,] 16 10 2 7
[80,] 11 20 3 6
[81,] 12 9 4 7
[82,] 16 6 3 8
[83,] 18 6 2 8
[84,] 12 12 4 6
[85,] 14 10 2 7
[86,] 15 1 4 8
[87,] 12 9 3 8
[88,] 12 10 4 6
[89,] 14 6 3 7
[90,] 15 0 5 7
[91,] 17 3 1 9
[92,] 14 7 2 7
[93,] 13 10 2 8
[94,] 13 5 4 7
[95,] 18 2 2 7
[96,] 15 5 3 7
[97,] 11 3 5 7
[98,] 11 6 4 6
[99,] 16 1 3 8
[100,] 13 7 3 6
[101,] 13 4 4 6
[102,] 15 9 2 6
[103,] 13 4 3 7
[104,] 14 4 3 7
[105,] 17 0 2 8
[106,] 13 3 2 8
[107,] 9 7 4 7
[108,] 11 4 4 6
[109,] 7 12 3 6
[110,] 9 7 4 6
[111,] 14 2 1 9
[112,] 13 15 0 6
[113,] 9 5 3 8
[114,] 11 5 4 6
[115,] 13 2 2 8
[116,] 9 4 4 6
[117,] 14 2 2 6
[118,] 13 8 0 8
[119,] 11 7 3 5
[120,] 16 1 1 7
[121,] 13 1 3 6
[122,] 11 4 2 7
[123,] 11 1 2 8
[124,] 13 2 2 6
[125,] 11 5 2 6
[126,] 14 0 2 7
[127,] 11 7 2 6
[128,] 8 8 3 6
[129,] 10 4 3 6
[130,] 13 0 2 6
[131,] 10 2 2 7
[132,] 10 1 3 6
[133,] 11 2 2 6
[134,] 11 0 2 6
[135,] 10 1 3 6
[136,] 11 1 2 6
[137,] 14 1 2 5
[138,] 10 1 1 7
[139,] 12 0 2 6
[140,] 10 0 2 6
[141,] 11 0 0 6
> sessionInfo()
R version 2.10.1 (2009-12-14)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] MiscPsycho_1.6 statmod_1.4.6 lattice_0.17-26 gdata_2.8.0
loaded via a namespace (and not attached):
[1] grid_2.10.1 gtools_2.6.2 tools_2.10.1
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
