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
day24
svn rev37909
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.
