As the error message says, the values of your function must be finite in
order to run the algorithm.
Some part of your loop is passing arguments (inits maybe... you only tried
(0,0) in the CLI example) that cause IRT.llZetaLambdaCorrNan to be
infinite.
--------------------------------------
Jonathan P. Daily
Technician - USGS Leetown Science Center
11649 Leetown Road
Kearneysville WV, 25430
(304) 724-4480
"Is the room still a room when its empty? Does the room,
the thing itself have purpose? Or do we, what's the word... imbue it."
- Jubal Early, Firefly
From:
Damokun <dmroi...@students.cs.uu.nl>
To:
r-help@r-project.org
Date:
11/03/2010 10:19 AM
Subject:
[R] optim works on command-line but not inside a function
Sent by:
r-help-boun...@r-project.org
Dear all,
I am trying to optimize a logistic function using optim, inside the
following functions:
#Estimating a and b from thetas and outcomes by ML
IRT.estimate.abFromThetaX <- function(t, X, inits, lw=c(-Inf,-Inf),
up=rep(Inf,2)){
optRes <- optim(inits, method="L-BFGS-B", fn=IRT.llZetaLambdaCorrNan,
gr=IRT.gradZL,
lower=lw, upper=up, t=t, X=X)
c(optRes$par[2], -(optRes$par[1]/optRes$par[2]) )
}
#Estimating a and b from thetas and outcomes by ML, avoiding 0*log(0)
IRT.estimate.abFromThetaX2 <- function(tar, Xes, inits, lw=c(-Inf,-Inf),
up=rep(Inf,2)){
optRes <- optim(inits, method="L-BFGS-B", fn=IRT.llZetaLambdaCorrNan,
gr=IRT.gradZL,
lower=lw, upper=up, t=tar, X=Xes)
c(optRes$par[2], -(optRes$par[1]/optRes$par[2]) )
}
The problem is that this does not work:
> IRT.estimate.abFromThetaX(sx, st, c(0,0))
Error in optim(inits, method = "L-BFGS-B", fn = IRT.llZetaLambdaCorrNan, :
L-BFGS-B needs finite values of 'fn'
But If I try the same optim call on the command line, with the same data,
it
works fine:
> optRes <- optim(c(0,0), method="L-BFGS-B", fn=IRT.llZetaLambdaCorrNan,
+ gr=IRT.gradZL,
+ lower=c(-Inf, -Inf), upper=c(Inf, Inf), t=st, X=sx)
> optRes
$par
[1] -0.6975157 0.7944972
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
Does anyone have an idea what this could be, and what I could try to avoid
this error? I tried bounding the parameters, with lower=c(-10, -10) and
upper=... but that made no difference.
Thanks,
Diederik Roijers
Utrecht University MSc student.
------
PS: the other functions I am using are:
#IRT.p is the function that represents the probability
#of a positive outcome of an item with difficulty b,
#discriminativity a, in combination with a student with
#competence theta.
IRT.p <- function(theta, a, b){
epow <- exp(-a*(theta-b))
result <- 1/(1+epow)
result
}
# = IRT.p^-1 ; for usage in the loglikelihood
IRT.oneOverP <- function(theta, a, b){
epow <- exp(-a*(theta-b))
result <- (1+epow)
result
}
# = (1-IRT.p)^-1 ; for usage in the loglikelihood
IRT.oneOverPneg <- function(theta, a, b){
epow <- exp(a*(theta-b))
result <- (1+epow)
result
}
#simulation-based sample generation of thetas and outcomes
#based on a given a and b. (See IRT.p) The sample-size is n
IRT.generateSample <- function(a, b, n){
x<-rnorm(n, mean=b, sd=b/2)
t<-IRT.p(x,a,b)
ch<-runif(length(t))
t[t>=ch]=1
t[t<ch]=0
cbind(x,t)
}
#This loglikelihood function is based on the a and be parameters,
#and requires thetas as input in X, and outcomes in t
#prone to give NaN errors due to 0*log(0)
IRT.logLikelihood2 <- function(params, t, X){
pos<- sum(t * log(IRT.p(X,params[1],params[2])))
neg<- sum( (1-t) * log( (1-IRT.p(X,params[1],params[2])) ) )
-pos-neg
}
#Avoiding NaN problems due to 0*log(0)
#otherwise equivalent to IRT.logLikelihood2
IRT.logLikelihood2CorrNan <- function(params, t, X){
pos<- sum(t * log(IRT.oneOverP(X,params[1],params[2])))
neg<- sum((1-t) * log(IRT.oneOverPneg(X,params[1],params[2])))
-pos-neg
}
#IRT.p can also be espressed in terms of z and l
#where z=-ab and l=a <- makes it a standard logit function
IRT.pZL <- function(theta, z, l){
epow <- exp(-(z+l*theta))
result <- 1/(1+epow)
result
}
#as IRT.oneOverP but now for IRT.pZL
IRT.pZLepos <- function(theta, z, l){
epow <- exp(-(z+l*theta))
result <- (1+epow)
result
}
#as IRT.oneOverPneg but now for IRT.pZL
IRT.pZLeneg <- function(theta, z, l){
epow <- exp(z+l*theta)
result <- (1+epow)
result
}
#The loglikelihood of IRT, but now expressed in terms of z and l
IRT.llZetaLambda <- function(params, t, X){
pos<- sum(t * log(IRT.pZL( X,params[1],params[2]) ))
neg<- sum( (1-t) * log( (1-IRT.pZL(X,params[1],params[2] )) ) )
-pos-neg
}
#Same as IRT.logLikelihood2CorrNan but for IRT.llZetaLambda
IRT.llZetaLambdaCorrNan <- function(params, t, X){
pos <- sum(t * log(IRT.pZLepos( X,params[1],params[2]) ))
neg <- sum((1-t) * log(IRT.pZLeneg(X,params[1],params[2]) ))
pos+neg
}
#Gradient of IRT.llZetaLambda
IRT.gradZL <- function(params, t, X){
res<-numeric(length(params))
res[1] <- sum(t-IRT.pZL( X,params[1],params[2] ))
res[2] <- sum(X*(t-IRT.pZL( X,params[1],params[2] )))
-res
}
#And to create the sample:
s <- IRT.generateSample(0.8, 1, 50)
sx <- s[,1]
st <- s[,2]
IRT.estimate.abFromThetaX(sx, st, c(0,0))
--
View this message in context:
http://r.789695.n4.nabble.com/optim-works-on-command-line-but-not-inside-a-function-tp3025414p3025414.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org 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.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org 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.
