Hi Harold, you're probably looking for something like:
rasch.max2 <- function(x, betas){ opt <- function(theta){ -sum(dbinom(x, 1, plogis(theta - betas), log = TRUE)) } out <- optim(log(sum(x)/(length(x)/sum(x))), opt, method = "BFGS", hessian = TRUE) cat('theta is about', round(out$par, 2), ', se', 1/sqrt(out$hes),'\n') } rasch.max(c(1, 1, 0, 0), c(-1, .5, 0, 1)) rasch.max2(c(1, 1, 0, 0), c(-1, .5, 0, 1)) rasch.max(c(1, 0, 1, 1), c(-1, .5, 0, 1)) rasch.max2(c(1, 0, 1, 1), c(-1, .5, 0, 1)) I hope it helps. Best, Dimitris ---- Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm ----- Original Message ----- From: "Doran, Harold" <[EMAIL PROTECTED]> To: <r-help@stat.math.ethz.ch> Sent: Thursday, August 24, 2006 2:54 PM Subject: [R] Optim question > This is a very basic question, but I am a bit confused with optim. I > want to get the MLEs using optim which could replace the > newton-raphson > code I have below which also gives the MLEs. The function takes as > input > a vector x denoting whether a respondent answered an item correctly > (x=1) or not (x=0). It also takes as input a vector b_vector, and > these > are parameters of test items (Rasch estimates in this case) > > For example, here is how my current function operates. > >> rasch.max(c(1,1,0,0), c(-1,.5,0,1)) > theta is about 0.14 , se 1.063972 > > I'm not quite sure how to accomplish the same thing using optim. Can > anyone offer a suggestion? > > rasch.max <- function(x, b_vector){ > p <- numeric(length(b_vector)) > theta <- log(sum(x)/(length(x)/sum(x))) # This is a starting value > for theta > rasch <- function(theta,b) 1/ (1 + exp(b-theta)) > old <- 0 > updated <- 5 > while(abs(old-updated) > .001){ > old <- updated > for(k in seq(along=b_vector)) p[k] <- rasch(theta,b_vector[k]) > first_deriv <- sum(x) - sum(p) > second_deriv <- sum((1-p)*-p) > change <- (first_deriv/second_deriv) > theta <- theta - change # This is the updated theta > updated <- change > } > cat('theta is about', round(theta,2), ', se', 1/sqrt(-second_deriv), > '\n') > } > > Harold > >> version > _ > platform i386-pc-mingw32 > arch i386 > os mingw32 > system i386, mingw32 > status > major 2 > minor 3.0 > year 2006 > month 04 > day 24 > svn rev 37909 > language R > version.string Version 2.3.0 (2006-04-24) > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@stat.math.ethz.ch 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. > Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm ______________________________________________ R-help@stat.math.ethz.ch 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.