[R] question on "optim"

[Hey Sky]

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.