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" <HDoran at air.org>
To: <r-help at 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 at 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