[R] optim works on command-line but not inside a function

[Damokun]

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.