[R] question on "optim"

[Ravi Varadhan]

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.