multiple parameter optimization with optim()

5 messages
Open this post in threaded view
|

multiple parameter optimization with optim()

 I am trying to generalize a working piece of code for a single parameter to a multiple parameter problem. Reproducible code is below. The parameters to be estimated are a, b, and c. The estimation problem is such that there is one set of a, b, c parameters for each column of the data. Hence, in this sample data with 20 columns, there are 20 a params, 20 b-params, and 20 c-params. Because I am estimating so many parameters, I am not certain that I have indicated to the function properly the right number of params to estimate and also if I have generated starting values in a sufficient way. Thanks for any help. Harold dat <- replicate(20, sample(c(0,1), 2000, replace = T)) library(stat mod) qq <- gauss.quad.prob(Q, dist = 'normal', mu = 0, sigma = 1) nds <- qq\$nodes wts <- qq\$weights fn <- function(params){ a <- params[1:ncol(dat)] b <- params[1:ncol(dat)] c <- params[1:ncol(dat)] L <- sapply(1:ncol(dat), function(i) dbinom(dat[,i], 1, c + ((1 - c)/(1 + exp(-1.7 * a * (nds[i] - b)))) * wts[i])) r1 <- prod(colSums(L * wts)) -log(r1) } startVal <- rep(.5, ncol(dat)) opt <- optim(startVal, fn)         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: multiple parameter optimization with optim()

 Some observations -- no solution here though: 1) the code is not executable. I tried. Maybe that makes it reproducible!     Typos such as "stat mod", undefined Q etc. 2) My experience is that any setup with a ?apply approach that doesn't then check to see that the structure of the data is correct has a high probability of failure due to mismatch with the optimizer requirements. It's worth being VERY pedestrian in setting up optimization functions and checking obsessively that you get what you expect and that there are no regions you are likely to wander into with divide by 0, log(negative), etc. 3) optim() is a BAD choice here. I wrote the source for three of the codes, and the one most appropriate for many parameters (CG) I have been deprecating for about 30 years. Use Rcgmin or something else instead. 4) If possible, analytic gradients are needed for CG like codes. You probably need to dig out some source code for dbinom() to do this, but your function is not particularly complicated, and doesn't have "if" statements etc. However, you could test a case using the numDeriv gradient that is an option for Rcgmin, but it will be painfully slow. For a one-off computation, that may still be acceptable. JN On 15-02-18 06:00 AM, [hidden email] wrote: > Message: 37 > Date: Tue, 17 Feb 2015 23:03:24 +0000 > From: "Doran, Harold" <[hidden email]> > To: "[hidden email]" <[hidden email]> > Subject: [R] multiple parameter optimization with optim() > Message-ID: > Content-Type: text/plain; charset="UTF-8" > > I am trying to generalize a working piece of code for a single parameter to a multiple parameter problem. Reproducible code is below. The parameters to be estimated are a, b, and c. The estimation problem is such that there is one set of a, b, c parameters for each column of the data. Hence, in this sample data with 20 columns, there are 20 a params, 20 b-params, and 20 c-params. > > Because I am estimating so many parameters, I am not certain that I have indicated to the function properly the right number of params to estimate and also if I have generated starting values in a sufficient way. > > Thanks for any help. > Harold > > dat <- replicate(20, sample(c(0,1), 2000, replace = T)) > library(stat mod) > qq <- gauss.quad.prob(Q, dist = 'normal', mu = 0, sigma = 1) > nds <- qq\$nodes > wts <- qq\$weights > fn <- function(params){ > a <- params[1:ncol(dat)] > b <- params[1:ncol(dat)] > c <- params[1:ncol(dat)] > L <- sapply(1:ncol(dat), function(i) dbinom(dat[,i], 1, c + ((1 - c)/(1 + exp(-1.7 * a * (nds[i] - b)))) * wts[i])) > r1 <- prod(colSums(L * wts)) > -log(r1) > } > startVal <- rep(.5, ncol(dat)) > opt <- optim(startVal, fn) > > [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|