[R] Optim question

[Dimitris Rizopoulos]

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