Well, the function should not be able to be infinite as IRT.llZetaLambdaCorrNan
is a sum of products of either one or zero and log(1+exp(x)) or
log(1+exp(-x)) (these logs are always bigger or equal to log(1)=0) Further
more, I bounded x to be finite to fix my problem (as I expected that it
might try x->Inf. But this did not help.

And it is a mystery to me why it would work on the command line, and not as
part of a function (it is just one call, and exactly the same one too.) (I
tried this in order to find proper intervals and start values where the
error would not arise. But to my surprise it just gave normal values when I
used the same settings as in the function.)

Thanks,
Diederik


On 3 November 2010 15:41, Roijers, D.M. <D.M.Roijers@students.uu.nl> wrote:

>
> -------------------------------------------
> *From:* Jonathan P Daily[SMTP:JDAILY@USGS.GOV <SMTP%3AJDAILY@USGS.GOV>]
> *Sent:* Wednesday, November 03, 2010 3:26:09 PM
> *To:* Damokun
> *Cc:* r-help@r-project.org; r-help-bounces@r-project.org
> *Subject:* Re: [R] optim works on command-line but not inside a function
> *Auto forwarded by a Rule*
>
>
> 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 <dmroijer@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-bounces@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<http://www.r-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>
>
>

	[[alternative HTML version deleted]]

