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.
Jonathan P Daily
2010-Nov-03 14:26 UTC
[R] optim works on command-line but not inside a function
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 and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]]
Diederik Roijers
2010-Nov-03 18:09 UTC
[R] FW: optim works on command-line but not inside a function
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]]
Berend Hasselman
2010-Nov-03 20:03 UTC
[R] optim works on command-line but not inside a function
Damokun wrote:> > 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" >In your command line you have set t=st and X=sx. However in the alternative you do: IRT.estimate.abFromThetaX(sx, st, c(0,0)) Therefore you are assigning sx to t and st to X in the IRT.estimate.abFromThetaX function, which is reversed from your command line call. You should switch sx and st in the function call: IRT.estimate.abFromThetaX(st, sx, c(0,0)) and then all will be well. best Berend If yoou -- View this message in context: http://r.789695.n4.nabble.com/optim-works-on-command-line-but-not-inside-a-function-tp3025414p3026099.html Sent from the R help mailing list archive at Nabble.com.