Re: [R] question about optim
thank you very much for your response!
But I have now another problem: optim( ) gives me the exactly same initial
values for the 4 parameters. I have tried many initial values, and it always
didnt change them.
opt$convergence = 1. I increased the number of iterations, but the problem is
still there. I also made sure that the function doesnt return NAN or INF
values...
Does anyone have any idea where the problem could come from ?
Many thanks.
To: r-h...@stat.math.ethz.ch
From: bbol...@gmail.com
Date: Tue, 12 Apr 2011 17:24:05 +
Subject: Re: [R] question about optim
chirine wolley wolley.chirine at hotmail.com writes:
Dear R-users, I would like to use optim( ) to minimize a function
which depends on 4 parameters: 2 vectors, a scalar, and a matrix.
And I have a hard to define the parameters at the beginning of the
function, and then to call optim. Indeed, all the examples I have
seen dont treat cases where parameters are not all real. Here is my
code, it doesnt work but its just to show you where is exactly my
problem:
suppose the dimensions of X are (Xc,Xr)
and the dimensions of x3 (I can't easily figure out what they should
be) are (m,n)
g=function(x,matrice)
{
x1 = x[1:Xc] # 1st parameter x1 is a vector
x2 = x[Xc+1] # 2nd parameter x2 is real
x3 - matrix(x[(Xc+2):(Xc+1+m*n)],nrow=m,ncol=n)
## be careful about column- vs row-encoding, possibly use byrow=TRUE
x4 - x[(Xc+1+m*n+1):length(x)] #4th paramater is a vector
res1=rep(0,nrow(X))
res2=matrix(0,nrow=nrow(X),ncol=ncol(Y))
for (i in 1:nrow(X))
{
res1[i]=log(1/(1+exp(-t(x1)%*%X[i,]-x2)))
for (t in 1:ncol(Y))
{
res2[i,t]=log(((1+exp(-t(x3[,t])%*%X[i,]-x4[t]))/(sqrt(2*pi)))*
(exp(-0.5*((1+exp(-t(x3[,t])%*%X[i,]-x4[t]))*(1-Y[i,t]))^(2
}
}
sum(res1)+sum(res2) # the function to minimize
}
to call optim(), concatenate reasonable initial guesses for your
parameters. (I have no idea what 'c' is in your code below ...)
opt=optim(c(alpha[,c+1],beta[c+1],w,gamma),g)
### and how can we call optim ???
__
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.
[R] question about optim
Dear R-users,
I would like to use optim( ) to minimize a function which depends on 4
parameters: 2 vectors, a scalar, and a matrix.
And I have a hard to define the parameters at the beginning of the function,
and then to call optim. Indeed, all the examples I have seen dont treat cases
where parameters are not all real.
Here is my code, it doesnt work but its just to show you where is exactly my
problem:
g=function(x,matrice)
{
x1 = x[1:ncol(X)] # 1st parameter x1 is a vector
x2 = x[ncol(X)+1]# 2nd parameter x2 is real
x3 = x(0,166,5) # 3rd parameter is a matrix How should one write it ???
x4 = x[ncol(X)+1:ncol(X)+1+ncol(Y)] #4th paramater is a vector
res1=rep(0,nrow(X))
res2=matrix(0,nrow=nrow(X),ncol=ncol(Y))
for (i in 1:nrow(X))
{
res1[i]=log(1/(1+exp(-t(x1)%*%X[i,]-x2)))
for (t in 1:ncol(Y))
{
res2[i,t]=log(((1+exp(-t(x3[,t])%*%X[i,]-x4[t]))/(sqrt(2*pi)))*
(exp(-0.5*((1+exp(-t(x3[,t])%*%X[i,]-x4[t]))*(1-Y[i,t]))^(2
}
}
sum(res1)+sum(res2)# the function to minimize
}
opt=optim(c(alpha[,c+1],beta[c+1],w,gamma),g) ### and how can we call optim ???
Thank you in advance!
[[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.
Re: [R] question about optim
chirine wolley wolley.chirine at hotmail.com writes:
Dear R-users, I would like to use optim( ) to minimize a function
which depends on 4 parameters: 2 vectors, a scalar, and a matrix.
And I have a hard to define the parameters at the beginning of the
function, and then to call optim. Indeed, all the examples I have
seen dont treat cases where parameters are not all real. Here is my
code, it doesnt work but its just to show you where is exactly my
problem:
suppose the dimensions of X are (Xc,Xr)
and the dimensions of x3 (I can't easily figure out what they should
be) are (m,n)
g=function(x,matrice)
{
x1 = x[1:Xc] # 1st parameter x1 is a vector
x2 = x[Xc+1]# 2nd parameter x2 is real
x3 - matrix(x[(Xc+2):(Xc+1+m*n)],nrow=m,ncol=n)
## be careful about column- vs row-encoding, possibly use byrow=TRUE
x4 - x[(Xc+1+m*n+1):length(x)] #4th paramater is a vector
res1=rep(0,nrow(X))
res2=matrix(0,nrow=nrow(X),ncol=ncol(Y))
for (i in 1:nrow(X))
{
res1[i]=log(1/(1+exp(-t(x1)%*%X[i,]-x2)))
for (t in 1:ncol(Y))
{
res2[i,t]=log(((1+exp(-t(x3[,t])%*%X[i,]-x4[t]))/(sqrt(2*pi)))*
(exp(-0.5*((1+exp(-t(x3[,t])%*%X[i,]-x4[t]))*(1-Y[i,t]))^(2
}
}
sum(res1)+sum(res2)# the function to minimize
}
to call optim(), concatenate reasonable initial guesses for your
parameters. (I have no idea what 'c' is in your code below ...)
opt=optim(c(alpha[,c+1],beta[c+1],w,gamma),g)
### and how can we call optim ???
__
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.
Re: [R] question on optim() fn.
As I said before, you get NaNs whenever the iterates are such that they are mathematically infeasible. This could very much depend on the starting value for the optimizer, as you have observed, because the trajectories of the optimizer can be widely different when started from diffrent parameter values. KKT stands for Karash-Kuhn-Tucker conditions. They provide the standard, necessary conditions for a parameter vector to be a local optimum in smooth problems. The first order conditions relates to the stationarity of the objective function and can be tested using its gradient (or the augmented Lagrangian in the case of constrained optimization), and the second order condition pertains to the curvature of the objective function (or the augmented Lagrangian). In convex optimization, KKT conditions are also sufficient. optim() does not test for the KKT conditions. You can use optimx() to do this, although the results of the tests of first and second order conditions are sensitive to scaling of parameters and objective function. Ravi. Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu - Original Message - From: sharkoil whitewhale2...@gmail.com Date: Thursday, October 21, 2010 1:09 am Subject: Re: [R] question on optim() fn. To: r-help@r-project.org Hi Ravi, Thanks a lot for your informations. I think you are right. Without the constraints of the parameters, it's easy to get the NaN Produced warning message. The pdfs I worked on are Burr, Lognormal and Inverse Gaussian distributions. Lognormal distribution is ok, but I always have warning message to get mle for Burr and Inverse Gaussian. Also I found out that different initial values make difference. Some will get the valid MLE without warning message, some will have warning message. Don't why! Probably I have to try optimx() function. By the way, what is KKT test? Does optim() function have this output I did get the MLE with convergence being 0. With your tips, I will assume that I have got the valid MLE even though I had the warning message. Also I tried to use proc nlmixed procedure in sas to estimate the MlE. It seems the optim() function and proc nlmixed give me the similar MLE. Regards -- View this message in context: Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code. __ 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.
Re: [R] question on optim() fn.
Hi, You are doing the right thing by being concerned about the warning messages. The warning message of `NaN' usually arises due to illegal arithmetic operations (in the real number field) such as logarithm or square root of a negative number. This could happen during the intermediate stages of an iterative procedure. If you do not like these warnings, or are concerned about the solutions, you should try and set bounds on the parameters based on your knowledge about the statistical model. Then you would use a box-constrained optimization algorithm. If you have more complicated constraints, you may need more specialized optimization algorithms. It is hard to tell whether your solution is valid or not without seeing your model and the output from optim. If your convergence indicator is 0, this usually (but not necessarily) means that you have a valid solution, although I would still check the first and second order KKT optimality conditions. You can also use `optimx' package, which provides a test of KKT conditions, and also allows the use of other optimizers in R. Ravi. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of sharkoil Sent: Tuesday, October 19, 2010 11:55 PM To: r-help@r-project.org Subject: [R] question on optim() fn. i am trying to use optim() fn in R to estimate mle of my pdf. even though I am able to get the results, but there are always warning message, which says that NaN produced. I am not very sure if I should ignore these warning message since I have got solution without error message. I am wondering if these message is an indicator of bad mle. How I can debug to get rid of these warning message? thanks a lot. -- View this message in context: http://r.789695.n4.nabble.com/question-on-optim-fn-tp3003217p3003217.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. __ 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.
Re: [R] question on optim() fn.
Hi Ravi, Thanks a lot for your informations. I think you are right. Without the constraints of the parameters, it's easy to get the NaN Produced warning message. The pdfs I worked on are Burr, Lognormal and Inverse Gaussian distributions. Lognormal distribution is ok, but I always have warning message to get mle for Burr and Inverse Gaussian. Also I found out that different initial values make difference. Some will get the valid MLE without warning message, some will have warning message. Don't why! Probably I have to try optimx() function. By the way, what is KKT test? Does optim() function have this output I did get the MLE with convergence being 0. With your tips, I will assume that I have got the valid MLE even though I had the warning message. Also I tried to use proc nlmixed procedure in sas to estimate the MlE. It seems the optim() function and proc nlmixed give me the similar MLE. Regards -- View this message in context: http://r.789695.n4.nabble.com/question-on-optim-fn-tp3003217p3005004.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.
[R] question on optim() fn.
i am trying to use optim() fn in R to estimate mle of my pdf. even though I am able to get the results, but there are always warning message, which says that NaN produced. I am not very sure if I should ignore these warning message since I have got solution without error message. I am wondering if these message is an indicator of bad mle. How I can debug to get rid of these warning message? thanks a lot. -- View this message in context: http://r.789695.n4.nabble.com/question-on-optim-fn-tp3003217p3003217.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.
Re: [R] question on optim
thanks. Ravi and Nash. I will read the new package and may use it after I am familiar with it. I may bother both of you when I have questions.thanks for that in advance. Nan from Montreal Hi Nan, You can take a look at the optimx package on CRAN. John Nash and I wrote this package to help lay and sophisticated users alike. This package unifies various optimization algorithms in R for smooth, box-constrained optimization. It has features for checking objective function, gradient (and hessian) specifications. It checks for potential problems due to poor scaling; checks feasibility of starting values. It provides diagnostics (KKT conditions) on whether or not a local optimum has been located. It also allows the user to run various optimization algorithms in one simple call, which is essentially identical to optim call. This feature can be especially useful for developers to benchmark different algorithms and choose the best one for their class of problems. http://cran.r-project.org/web/packages/optimx/index.html Ravi. __ 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.
Re: [R] question on optim
Hi Nan, You can take a look at the optimx package on CRAN. John Nash and I wrote this package to help lay and sophisticated users alike. This package unifies various optimization algorithms in R for smooth, box-constrained optimization. It has features for checking objective function, gradient (and hessian) specifications. It checks for potential problems due to poor scaling; checks feasibility of starting values. It provides diagnostics (KKT conditions) on whether or not a local optimum has been located. It also allows the user to run various optimization algorithms in one simple call, which is essentially identical to optim call. This feature can be especially useful for developers to benchmark different algorithms and choose the best one for their class of problems. http://cran.r-project.org/web/packages/optimx/index.html Ravi. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Hey Sky Sent: Tuesday, September 07, 2010 2:48 PM To: Ben Bolker; r-h...@stat.math.ethz.ch Subject: Re: [R] question on optim thanks. Ben after read your email, I realized the initial value of w[5]=0 is a stupid mistake. and I have changed it. but I am sorry I cannot reproduce the result, convergence, as you get. the error message is non-finite finite difference value [12]. any suggestion about it? and could you plz recommend some R books on optimization, such as tips for setup gradient and others, or common mistakes? thanks Nan - Original Message From: Ben Bolker bbol...@gmail.com To: r-h...@stat.math.ethz.ch Sent: Tue, September 7, 2010 11:15:43 AM Subject: Re: [R] question on quot;optimquot; Hey Sky heyskywalker at yahoo.com writes: I do not know how to describe my question. I am a new user for R and write the following code for a dynamic labor economics model and use OPTIM to get optimizations and parameter values. the following code does not work due to the equation: wden[,i]-dnorm((1-regw[,i])/w[5])/w[5] where w[5] is one of the parameters (together with vector a, b and other elements in vector w) need to be estimated. if I delete the w[5] from the upper equation. that is: wden[,i]-dnorm(1-regw[,i]) optim will give me the estimated parameters. Thank you for the reproducible example! The problem is that you are setting the initial value of w[5] to zero, and then trying to divide by it ... I find that guess-rep(0,times=npar) guess[16] - 1 system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=TRUE, control=list(trace=TRUE))) seems to work OK (I have no idea if the answers are sensible are not ...) If you're going to be doing a lot of this it might be wise to see if you can specify the gradient of your objective function for R -- it will speed up and stabilize the fitting considerably. By the way, you should be careful with this function: if we try this with Nelder-Mead instead, it appears to head to a set of parameters that lead to some sort of singularity in the objective function: system.time(r2-optim(guess,myfunc1,data=mydata, method=Nelder-Mead,hessian=FALSE, control=list(trace=TRUE,maxit=5000))) ## still thinks it hasn't converged, but objective function is ## much smaller ## plot 'slice' through objective space where 0 corresponds to ## fit-1 parameters and 1 corresponds to fit-2 parameters; ## adapted from emdbook::calcslice range - seq(-0.1,1.1,length=400) slicep - seq(range[1], range[2], length = 400) slicepars - t(sapply(slicep, function(x) (1 - x) * r1$par + x * r2$par)) v - apply(slicepars, 1, myfunc1) plot(range,v,type=l) Ideally, you should be able to look at the parameters of fit #2 and figure out (a) what the result means in terms of labor economics and (b) how to keep the objective function from going there, or at least identifying when it does. Ben Bolker __ 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. __ 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. __ 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.
[R] question on optim
Ben and Ravi have already pointed out that you have a problem with a non-computable (e.g., divide by zero or similar) objective function. This is common. optim() has mostly unconstrained optimizers. L-BFGS-B does handle box or bounds constraints. There is also Rvmmin, which is supposedly the same algorithm as optim's BFGS (I'm the author, but not the implementor of optim). However, Rvmmin has bounds constraints, so you could constrain the parameter going to zero to be .1 or something like that. You don't want to use a bound that would cause the computational failure, because the numerical gradient routine takes a step, essentially walking off the cliff. If you have analytic gradients, you can do much better usually. Note that bounds in Nelder-Mead are fairly awkward to use, and our codes are not very good for that. You might also try bobyqa from the minqa package, but be warned that it is an F1 race car relative to the Chevy of Rvmmin. Very good when set up right, but sometimes tricky to set up. John Nash Message: 33 Date: Tue, 7 Sep 2010 07:38:55 -0700 (PDT) From: Hey Sky heyskywal...@yahoo.com To: R r-help@r-project.org Subject: [R] question on optim Message-ID: 741050.53212...@web113920.mail.gq1.yahoo.com Content-Type: text/plain; charset=iso-8859-1 Hey, R users I do not know how to describe my question. I am a new user for R and write the following?code for a dynamic labor economics?model and use OPTIM to get optimizations and parameter values. the following code does not work due to the?equation: ?? wden[,i]-dnorm((1-regw[,i])/w[5])/w[5] where w[5]?is one of the parameters (together with vector a, b and other elements in vector w)?need to be estimated. if I delete the w[5] from the upper equation. that is: ?wden[,i]-dnorm(1-regw[,i]) optim will give me the estimated parameters. but the paramter w[5] is a key one for me. now my questions is: what reason lead to the first?equation does not work and the way to correct it to make it work? ?I wish I made the queston clear and thanks for any suggestion. Nan from Montreal __ 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.
Re: [R] question on optim
Dear R-list members, I am using R 2.11.1 on Windows XP. When I try to install package optimx through the GUI menu Packages / Install packages, this package does not appear in the list that opens up (chosen from the Austria CRAN site). The package is listed on Austria's CRAN web page, but today (8 September 2010) it does not show in the list obtained through the menu. Thank you. Paulo Barata (Rio de Janeiro - Brazil) -- -- On 8/9/2010 11:01, Ravi Varadhan wrote: Hi Nan, You can take a look at the optimx package on CRAN. John Nash and I wrote this package to help lay and sophisticated users alike. This package unifies various optimization algorithms in R for smooth, box-constrained optimization. It has features for checking objective function, gradient (and hessian) specifications. It checks for potential problems due to poor scaling; checks feasibility of starting values. It provides diagnostics (KKT conditions) on whether or not a local optimum has been located. It also allows the user to run various optimization algorithms in one simple call, which is essentially identical to optim call. This feature can be especially useful for developers to benchmark different algorithms and choose the best one for their class of problems. http://cran.r-project.org/web/packages/optimx/index.html Ravi. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Hey Sky Sent: Tuesday, September 07, 2010 2:48 PM To: Ben Bolker; r-h...@stat.math.ethz.ch Subject: Re: [R] question on optim thanks. Ben after read your email, I realized the initial value of w[5]=0 is a stupid mistake. and I have changed it. but I am sorry I cannot reproduce the result, convergence, as you get. the error message isnon-finite finite difference value [12]. any suggestion about it? and could you plz recommend some R books on optimization, such as tips for setup gradient and others, or common mistakes? thanks Nan - Original Message From: Ben Bolkerbbol...@gmail.com To: r-h...@stat.math.ethz.ch Sent: Tue, September 7, 2010 11:15:43 AM Subject: Re: [R] question onquot;optimquot; Hey Skyheyskywalkerat yahoo.com writes: I do not know how to describe my question. I am a new user for R and write the following code for a dynamic labor economics model and use OPTIM to get optimizations and parameter values. the following code does not work due to the equation: wden[,i]-dnorm((1-regw[,i])/w[5])/w[5] where w[5] is one of the parameters (together with vector a, b and other elements in vector w) need to be estimated. if I delete the w[5] from the upper equation. that is: wden[,i]-dnorm(1-regw[,i]) optim will give me the estimated parameters. Thank you for the reproducible example! The problem is that you are setting the initial value of w[5] to zero, and then trying to divide by it ... I find that guess-rep(0,times=npar) guess[16]- 1 system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=TRUE, control=list(trace=TRUE))) seems to work OK (I have no idea if the answers are sensible are not ...) If you're going to be doing a lot of this it might be wise to see if you can specify the gradient of your objective function for R -- it will speed up and stabilize the fitting considerably. By the way, you should be careful with this function: if we try this with Nelder-Mead instead, it appears to head to a set of parameters that lead to some sort of singularity in the objective function: system.time(r2-optim(guess,myfunc1,data=mydata, method=Nelder-Mead,hessian=FALSE, control=list(trace=TRUE,maxit=5000))) ## still thinks it hasn't converged, but objective function is ## much smaller ## plot 'slice' through objective space where 0 corresponds to ## fit-1 parameters and 1 corresponds to fit-2 parameters; ## adapted from emdbook::calcslice range- seq(-0.1,1.1,length=400) slicep- seq(range[1], range[2], length = 400) slicepars- t(sapply(slicep, function(x) (1 - x) * r1$par + x * r2$par)) v- apply(slicepars, 1, myfunc1) plot(range,v,type=l) Ideally, you should be able to look at the parameters of fit #2 and figure out (a) what the result means in terms of labor economics and (b) how to keep the objective function from going there, or at least identifying when it does. Ben Bolker __ 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. __ 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
Re: [R] question on optim
The windows build is not on CRAN at the moment. It might be there in a few days. In the meanwhile you can get the windows binary from the optimizer project on r-forge: https://r-forge.r-project.org/R/?group_id=395 Ravi. -Original Message- From: Paulo Barata [mailto:pbar...@infolink.com.br] Sent: Wednesday, September 08, 2010 4:26 PM To: Ravi Varadhan Cc: 'Hey Sky'; 'Ben Bolker'; r-h...@stat.math.ethz.ch Subject: Re: [R] question on optim Dear R-list members, I am using R 2.11.1 on Windows XP. When I try to install package optimx through the GUI menu Packages / Install packages, this package does not appear in the list that opens up (chosen from the Austria CRAN site). The package is listed on Austria's CRAN web page, but today (8 September 2010) it does not show in the list obtained through the menu. Thank you. Paulo Barata (Rio de Janeiro - Brazil) -- -- On 8/9/2010 11:01, Ravi Varadhan wrote: Hi Nan, You can take a look at the optimx package on CRAN. John Nash and I wrote this package to help lay and sophisticated users alike. This package unifies various optimization algorithms in R for smooth, box-constrained optimization. It has features for checking objective function, gradient (and hessian) specifications. It checks for potential problems due to poor scaling; checks feasibility of starting values. It provides diagnostics (KKT conditions) on whether or not a local optimum has been located. It also allows the user to run various optimization algorithms in one simple call, which is essentially identical to optim call. This feature can be especially useful for developers to benchmark different algorithms and choose the best one for their class of problems. http://cran.r-project.org/web/packages/optimx/index.html Ravi. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Hey Sky Sent: Tuesday, September 07, 2010 2:48 PM To: Ben Bolker; r-h...@stat.math.ethz.ch Subject: Re: [R] question on optim thanks. Ben after read your email, I realized the initial value of w[5]=0 is a stupid mistake. and I have changed it. but I am sorry I cannot reproduce the result, convergence, as you get. the error message isnon-finite finite difference value [12]. any suggestion about it? and could you plz recommend some R books on optimization, such as tips for setup gradient and others, or common mistakes? thanks Nan - Original Message From: Ben Bolkerbbol...@gmail.com To: r-h...@stat.math.ethz.ch Sent: Tue, September 7, 2010 11:15:43 AM Subject: Re: [R] question onquot;optimquot; Hey Skyheyskywalkerat yahoo.com writes: I do not know how to describe my question. I am a new user for R and write the following code for a dynamic labor economics model and use OPTIM to get optimizations and parameter values. the following code does not work due to the equation: wden[,i]-dnorm((1-regw[,i])/w[5])/w[5] where w[5] is one of the parameters (together with vector a, b and other elements in vector w) need to be estimated. if I delete the w[5] from the upper equation. that is: wden[,i]-dnorm(1-regw[,i]) optim will give me the estimated parameters. Thank you for the reproducible example! The problem is that you are setting the initial value of w[5] to zero, and then trying to divide by it ... I find that guess-rep(0,times=npar) guess[16]- 1 system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=TRUE, control=list(trace=TRUE))) seems to work OK (I have no idea if the answers are sensible are not ...) If you're going to be doing a lot of this it might be wise to see if you can specify the gradient of your objective function for R -- it will speed up and stabilize the fitting considerably. By the way, you should be careful with this function: if we try this with Nelder-Mead instead, it appears to head to a set of parameters that lead to some sort of singularity in the objective function: system.time(r2-optim(guess,myfunc1,data=mydata, method=Nelder-Mead,hessian=FALSE, control=list(trace=TRUE,maxit=5000))) ## still thinks it hasn't converged, but objective function is ## much smaller ## plot 'slice' through objective space where 0 corresponds to ## fit-1 parameters and 1 corresponds to fit-2 parameters; ## adapted from emdbook::calcslice range- seq(-0.1,1.1,length=400) slicep- seq(range[1], range[2], length = 400) slicepars- t(sapply(slicep, function(x) (1 - x) * r1$par + x * r2$par)) v- apply(slicepars, 1, myfunc1) plot(range,v,type=l) Ideally, you should be able to look at the parameters of fit #2 and figure out (a) what the result means in terms of labor economics and (b) how to keep the objective function from going
[R] question on optim
Hey, R users
I do not know how to describe my question. I am a new user for R and write the
following code for a dynamic labor economics model and use OPTIM to get
optimizations and parameter values. the following code does not work due to
the equation:
wden[,i]-dnorm((1-regw[,i])/w[5])/w[5]
where w[5] is one of the parameters (together with vector a, b and other
elements in vector w) need to be estimated. if I delete the w[5] from the upper
equation. that is:
wden[,i]-dnorm(1-regw[,i])
optim will give me the estimated parameters. but the paramter w[5] is a key one
for me.
now my questions is: what reason lead to the first equation does not work and
the way to correct it to make it work?
I wish I made the queston clear and thanks for any suggestion.
Nan
from Montreal
#the function
myfunc1-function(par,data) {
# define the parameter matrix used in following part
vbar2-matrix(0,n,nt)
vbar3-matrix(0,n,nt)
v8 -matrix(0,n,nt)
regw-matrix(0,n,nt)
wden-matrix(0,n,nt)
lia-matrix(0,n,nt)
ccl-matrix(1,n,ns)
eta-c(0,0)
# setup the parts for loglikelihood
q1-exp(par[1])
pr1-q1/(1+q1)
pr2-1-pr1
eta[2]-par[2]
a-par[3:6]
b-par[7:11]
w-par[12:npar]
for(m in 1:ns){
for(i in 1:nt){
regw[,i]-w[1]+ w[2]*eta[m]+exp(w[3]+w[4]*eta[m])*actr[,i]+w[5]*acwrk[,i]
vbar2[,i]=a[1]+ eta[m]+regw[,i]*a[2]+acwrk[,i]*a[3]+actr[,i]*a[4]
vbar3[,i]=b[1]+b[2]*eta[m]+regw[,i]*b[3]+acwrk[,i]*b[4]+actr[,i]*b[5]
v8[,i]=1+exp(vbar2[,i])+exp(vbar3[,i])
for(j in 1:n){
if (home[j,i]==1) lia[j,i]=1/v8[j,i]
if (wrk[j,i]==1) lia[j,i]=exp(vbar2[j,i])/v8[j,i]
if (tr[j,i]==1) lia[j,i]=exp(vbar3[j,i])/v8[j,i]
}
wden[,i]-dnorm((1-regw[,i])/w[5])/w[5]
ccl[,m]-lia[,i]*ccl[,m]*wden[,i]
}
}
func-pr1*ccl[,1]+pr2*ccl[,2]
f-sum(log(func))
return(-f)
}
#***
# main program ** gen random data and get the optimization **
nt-16 # number of periods
ns-2 # number of person type
n-50 # number of people
npar-17 # number of parameters
tr-matrix(0,n,nt)
wrk-matrix(0,n,nt)
home-matrix(0,n,nt)
actr-matrix(0,n,nt)
acwrk-matrix(0,n,nt)
for(i in 1:nt){
tr[,i]-round(runif(n))
home[,i]-round(runif(n))
}
for(i in 1:nt){
for(k in 1:n){
if(tr[k,i]==1 home[k,i]==1) home[k,i]=0
wrk[k,i]- 1-tr[k,i]-home[k,i]
}
}
actr[,1]-tr[,1]
acwrk[,1]-wrk[,1]
for(j in 2:nt){
actr[,j]-actr[,j-1]+tr[,j]
acwrk[,j]-acwrk[,j-1]+wrk[,j]
}
mydata-cbind(tr,wrk,home,actr,acwrk)
guess-rep(0,times=npar)
system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=T))
r1$par
__
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.
Re: [R] question on optim
sorry. there is a type in the following code. there is no w[5]*acwrk[,i] in the
regw equation. the right one should be as following:
regw[,i]-w[1]+ w[2]*eta[m]+exp(w[3]+w[4]*eta[m])*actr[,i]
- Original Message
From: Hey Sky heyskywal...@yahoo.com
To: R r-help@r-project.org
Sent: Tue, September 7, 2010 10:38:55 AM
Subject: [R] question on optim
Hey, R users
I do not know how to describe my question. I am a new user for R and write the
following code for a dynamic labor economics model and use OPTIM to get
optimizations and parameter values. the following code does not work due to
the equation:
wden[,i]-dnorm((1-regw[,i])/w[5])/w[5]
where w[5] is one of the parameters (together with vector a, b and other
elements in vector w) need to be estimated. if I delete the w[5] from the upper
equation. that is:
wden[,i]-dnorm(1-regw[,i])
optim will give me the estimated parameters. but the paramter w[5] is a key one
for me.
now my questions is: what reason lead to the first equation does not work and
the way to correct it to make it work?
I wish I made the queston clear and thanks for any suggestion.
Nan
from Montreal
#the function
myfunc1-function(par,data) {
# define the parameter matrix used in following part
vbar2-matrix(0,n,nt)
vbar3-matrix(0,n,nt)
v8 -matrix(0,n,nt)
regw-matrix(0,n,nt)
wden-matrix(0,n,nt)
lia-matrix(0,n,nt)
ccl-matrix(1,n,ns)
eta-c(0,0)
# setup the parts for loglikelihood
q1-exp(par[1])
pr1-q1/(1+q1)
pr2-1-pr1
eta[2]-par[2]
a-par[3:6]
b-par[7:11]
w-par[12:npar]
for(m in 1:ns){
for(i in 1:nt){
regw[,i]-w[1]+ w[2]*eta[m]+exp(w[3]+w[4]*eta[m])*actr[,i]+w[5]*acwrk[,i]
vbar2[,i]=a[1]+ eta[m]+regw[,i]*a[2]+acwrk[,i]*a[3]+actr[,i]*a[4]
vbar3[,i]=b[1]+b[2]*eta[m]+regw[,i]*b[3]+acwrk[,i]*b[4]+actr[,i]*b[5]
v8[,i]=1+exp(vbar2[,i])+exp(vbar3[,i])
for(j in 1:n){
if (home[j,i]==1) lia[j,i]=1/v8[j,i]
if (wrk[j,i]==1) lia[j,i]=exp(vbar2[j,i])/v8[j,i]
if (tr[j,i]==1) lia[j,i]=exp(vbar3[j,i])/v8[j,i]
}
wden[,i]-dnorm((1-regw[,i])/w[5])/w[5]
ccl[,m]-lia[,i]*ccl[,m]*wden[,i]
}
}
func-pr1*ccl[,1]+pr2*ccl[,2]
f-sum(log(func))
return(-f)
}
#***
# main program ** gen random data and get the optimization **
nt-16 # number of periods
ns-2 # number of person type
n-50 # number of people
npar-17 # number of parameters
tr-matrix(0,n,nt)
wrk-matrix(0,n,nt)
home-matrix(0,n,nt)
actr-matrix(0,n,nt)
acwrk-matrix(0,n,nt)
for(i in 1:nt){
tr[,i]-round(runif(n))
home[,i]-round(runif(n))
}
for(i in 1:nt){
for(k in 1:n){
if(tr[k,i]==1 home[k,i]==1) home[k,i]=0
wrk[k,i]- 1-tr[k,i]-home[k,i]
}
}
actr[,1]-tr[,1]
acwrk[,1]-wrk[,1]
for(j in 2:nt){
actr[,j]-actr[,j-1]+tr[,j]
acwrk[,j]-acwrk[,j-1]+wrk[,j]
}
mydata-cbind(tr,wrk,home,actr,acwrk)
guess-rep(0,times=npar)
system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=T))
r1$par
__
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.
__
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.
Re: [R] question on optim
Hi,
I do not see how `data' is used in your objective function.
The objective function is not even evaluable at the initial guess.
myfunc1(guess, mydata)
[1] NaN
I also think that some of the parameters may have to be constrained, for
example, par[1] 0. At a minimum, make sure that the obj fn is correctly
specified before we can start tackling other issues.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Hey Sky heyskywal...@yahoo.com
Date: Tuesday, September 7, 2010 10:40 am
Subject: [R] question on optim
To: R r-help@r-project.org
Hey, R users
I do not know how to describe my question. I am a new user for R and
write the
following code for a dynamic labor economics model and use OPTIM to
get
optimizations and parameter values. the following code does not work
due to
the equation:
wden[,i]-dnorm((1-regw[,i])/w[5])/w[5]
where w[5] is one of the parameters (together with vector a, b and
other
elements in vector w) need to be estimated. if I delete the w[5] from
the upper
equation. that is:
wden[,i]-dnorm(1-regw[,i])
optim will give me the estimated parameters. but the paramter w[5] is
a key one
for me.
now my questions is: what reason lead to the first equation does not
work and
the way to correct it to make it work?
I wish I made the queston clear and thanks for any suggestion.
Nan
from Montreal
#the function
myfunc1-function(par,data) {
# define the parameter matrix used in following part
vbar2-matrix(0,n,nt)
vbar3-matrix(0,n,nt)
v8 -matrix(0,n,nt)
regw-matrix(0,n,nt)
wden-matrix(0,n,nt)
lia-matrix(0,n,nt)
ccl-matrix(1,n,ns)
eta-c(0,0)
# setup the parts for loglikelihood
q1-exp(par[1])
pr1-q1/(1+q1)
pr2-1-pr1
eta[2]-par[2]
a-par[3:6]
b-par[7:11]
w-par[12:npar]
for(m in 1:ns){
for(i in 1:nt){
regw[,i]-w[1]+ w[2]*eta[m]+exp(w[3]+w[4]*eta[m])*actr[,i]+w[5]*acwrk[,i]
vbar2[,i]=a[1]+ eta[m]+regw[,i]*a[2]+acwrk[,i]*a[3]+actr[,i]*a[4]
vbar3[,i]=b[1]+b[2]*eta[m]+regw[,i]*b[3]+acwrk[,i]*b[4]+actr[,i]*b[5]
v8[,i]=1+exp(vbar2[,i])+exp(vbar3[,i])
for(j in 1:n){
if (home[j,i]==1) lia[j,i]=1/v8[j,i]
if (wrk[j,i]==1) lia[j,i]=exp(vbar2[j,i])/v8[j,i]
if (tr[j,i]==1) lia[j,i]=exp(vbar3[j,i])/v8[j,i]
}
wden[,i]-dnorm((1-regw[,i])/w[5])/w[5]
ccl[,m]-lia[,i]*ccl[,m]*wden[,i]
}
}
func-pr1*ccl[,1]+pr2*ccl[,2]
f-sum(log(func))
return(-f)
}
#***
# main program ** gen random data and get the optimization **
nt-16 # number of periods
ns-2 # number of person type
n-50 # number of people
npar-17 # number of parameters
tr-matrix(0,n,nt)
wrk-matrix(0,n,nt)
home-matrix(0,n,nt)
actr-matrix(0,n,nt)
acwrk-matrix(0,n,nt)
for(i in 1:nt){
tr[,i]-round(runif(n))
home[,i]-round(runif(n))
}
for(i in 1:nt){
for(k in 1:n){
if(tr[k,i]==1 home[k,i]==1) home[k,i]=0
wrk[k,i]- 1-tr[k,i]-home[k,i]
}
}
actr[,1]-tr[,1]
acwrk[,1]-wrk[,1]
for(j in 2:nt){
actr[,j]-actr[,j-1]+tr[,j]
acwrk[,j]-acwrk[,j-1]+wrk[,j]
}
mydata-cbind(tr,wrk,home,actr,acwrk)
guess-rep(0,times=npar)
system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=T))
r1$par
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
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.
Re: [R] question on optim
thanks. Ben after read your email, I realized the initial value of w[5]=0 is a stupid mistake. and I have changed it. but I am sorry I cannot reproduce the result, convergence, as you get. the error message is non-finite finite difference value [12]. any suggestion about it? and could you plz recommend some R books on optimization, such as tips for setup gradient and others, or common mistakes? thanks Nan - Original Message From: Ben Bolker bbol...@gmail.com To: r-h...@stat.math.ethz.ch Sent: Tue, September 7, 2010 11:15:43 AM Subject: Re: [R] question on quot;optimquot; Hey Sky heyskywalker at yahoo.com writes: I do not know how to describe my question. I am a new user for R and write the following code for a dynamic labor economics model and use OPTIM to get optimizations and parameter values. the following code does not work due to the equation: wden[,i]-dnorm((1-regw[,i])/w[5])/w[5] where w[5] is one of the parameters (together with vector a, b and other elements in vector w) need to be estimated. if I delete the w[5] from the upper equation. that is: wden[,i]-dnorm(1-regw[,i]) optim will give me the estimated parameters. Thank you for the reproducible example! The problem is that you are setting the initial value of w[5] to zero, and then trying to divide by it ... I find that guess-rep(0,times=npar) guess[16] - 1 system.time(r1-optim(guess,myfunc1,data=mydata, method=BFGS,hessian=TRUE, control=list(trace=TRUE))) seems to work OK (I have no idea if the answers are sensible are not ...) If you're going to be doing a lot of this it might be wise to see if you can specify the gradient of your objective function for R -- it will speed up and stabilize the fitting considerably. By the way, you should be careful with this function: if we try this with Nelder-Mead instead, it appears to head to a set of parameters that lead to some sort of singularity in the objective function: system.time(r2-optim(guess,myfunc1,data=mydata, method=Nelder-Mead,hessian=FALSE, control=list(trace=TRUE,maxit=5000))) ## still thinks it hasn't converged, but objective function is ## much smaller ## plot 'slice' through objective space where 0 corresponds to ## fit-1 parameters and 1 corresponds to fit-2 parameters; ## adapted from emdbook::calcslice range - seq(-0.1,1.1,length=400) slicep - seq(range[1], range[2], length = 400) slicepars - t(sapply(slicep, function(x) (1 - x) * r1$par + x * r2$par)) v - apply(slicepars, 1, myfunc1) plot(range,v,type=l) Ideally, you should be able to look at the parameters of fit #2 and figure out (a) what the result means in terms of labor economics and (b) how to keep the objective function from going there, or at least identifying when it does. Ben Bolker __ 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. __ 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.
Re: [R] question on optim
On 10-09-07 02:48 PM, Hey Sky wrote: thanks. Ben after read your email, I realized the initial value of w[5]=0 is a stupid mistake. and I have changed it. but I am sorry I cannot reproduce the result, convergence, as you get. the error message is non-finite finite difference value [12]. any suggestion about it? When I 'fix' the objective function as you specified in your second message, I get into trouble too. and could you plz recommend some R books on optimization, such as tips for setup gradient and others, or common mistakes? thanks Nan I'm afraid I don't know of a great reference, although pp. 340-341 of http://www.math.mcmaster.ca/~bolker/emdbook/book.pdf do give some basic trouble-shooting suggestions (nothing about gradients, though, but the example in ?optim does use gradients). I would say that my best _general_ advice is to understand what all the parameters of your model mean: what are reasonable starting values and reasonable ranges? Use control(parscale) to tell R the approximate expected order of magnitude of each parameter. You may be able to keep the model from getting into trouble by using method=L-BFGS-B and bounding the parameters within a reasonable range. For more general troubleshooting: try different optimization methods (in particular, Nelder-Mead doesn't need to compute finite differences); trap non-finite (NA, NaN, Inf, -Inf) that occur in the function and report them and/or stop the function (see ?browser) and/or replace them with large finite values; use control(trace=TRUE) or add cat() statements to your function to see where the optimizer is trying to go. good luck, Ben Bolker [[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.
