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.
