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]]