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.
