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
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
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 <heyskywalker at yahoo.com>
To: R <r-help at 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 at 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.
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: rvaradhan at jhmi.edu ----- Original Message ----- From: Hey Sky <heyskywalker at yahoo.com> Date: Tuesday, September 7, 2010 10:40 am Subject: [R] question on "optim" To: R <r-help at 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 at r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.
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 <heyskywalker at yahoo.com> > To: R <r-help at r-project.org> > Subject: [R] question on "optim" > Message-ID: <741050.53212.qm at 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 >
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.
Maybe Matching Threads
- optim with BFGS--what may lead to this, a strange thing happened
- replace loops with matrix
- Marginal Effect larger than 1 for a binary variable (summary.Design after lrm)
- How to remove the quote "" in the data frame?
- Marginal Effects of continuous variable in lrm model (Design package)