[R] nlminb and optim
Ravi Varadhan
rvaradhan at jhmi.edu
Thu Sep 30 00:26:58 CEST 2010
Can you send your code and data as separate files so we can get it into R easily?
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvaradhan at jhmi.edu
----- Original Message -----
From: "Doran, Harold" <HDoran at air.org>
Date: Wednesday, September 29, 2010 10:36 am
Subject: [R] nlminb and optim
To: R-help <R-help at r-project.org>
> 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 at r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list