[R] problems with loop

From: Simon Ruegg <s.ruegg_at_access.unizh.ch>
Date: Sat 26 Aug 2006 - 22:47:16 EST


Dear all,  

I am trying to evaluate the optimisation behaviour of a function. Originally I have optimised a model with real data and got a set of parameters. Now I am creating simulated data sets based on these estimates. With these simulations I am estimating the parameters again to see how variable the estimation is. To this end I have written a loop which should generate a new simulated data set each time. However, the optimisation algorithm, which works fine if only one data set is used, does not recognise the simulated data in the loop. Can anyone tell me where the error is? The code is below.  

Thanks for your help  

Simon    

# The loop in which nothing works anymore. The NLL in the optim function appears not to recognise the "new" data set. The functions used by this loop are given below.  

sim.estim=function(r)

{

      res=matrix(nrow=r,ncol=5)

      for (s in 1:r)

      {

            new=new.set()  

Min=optim(c(0.5,0.5,0.01,0.1),NLL,method="L-BFGS-B",lower=c(0,0,0,0))

            res[s,1]=Min$par[1]

            res[s,2]=Min$par[2]

            res[s,3]=Min$par[3]

            res[s,4]=Min$par[4]

            res[s,5]=Min$value[1]

      }

return(res)

}
     

# function to create a new data set. Works fine on its own  

new.set=function()

{

      tmp=matrix(nrow=510,ncol=2)

      v=numeric()

      sim=Model_1(1.2761707, 0.1953354, 2.7351930, 0.1032929)  

      tmp[,1]=sample(sim[,1],510,replace=T)

      v=runif(510,0,1)  

      for (k in 1 : 510)

      {

            for (n in 1 : length(sim[,1]))

            {

                  if (tmp[k,1]==sim[n,1] && v[k]<=sim[n,2]){tmp[k,2]=1}

                  if(tmp[k,1]==sim[n,1] && v[k]>sim[n,2]) {tmp[k,2]=0}

            }

      }

return(tmp)

}
   

# The model to be fitted  

Model_1=function(beta0_mod,mu_mod,k_mod,I_mod)

{

      parms = c(beta0=beta0_mod, mu=mu_mod, k=k_mod)

      my.atol = 1e-6  

times1=c(0,0.002884615,0.003846154,0.005769231,0.009615385,0.01,0.019230769)

      times2=c(0.028846154,0.038461538,0.057692308, 0.076923077,0.083333333,0.096153846)  

times3=c(0.115384615,0.125,0.134615385,0.166666667,0.25,0.416666667,0.5)

      times4=c(0.58333333,0.666666667,0.75,0.83333333,0.916666667,1)

      times5=c( 1.166666667,1.25,1.5,1.6,1.75,2,2.5)

      times6=seq(3,20)

      times = c(times1,times2,times3,times4,times5,times6)

      ODE_sys = function(t, x, w) #w=c("beta0","mu","k")

             {

                  I=x[1]

                                   

            dI=-w["mu"]*I+w["beta0"]*exp(-w["k"]*t)*(1-I)

            

            list(c(dI),c(S=1-I))

             }

 

      output = lsoda(c(I_mod),times,ODE_sys, parms, rtol=1e-6, atol=
my.atol)

      return(output)

}
 

# The negative log likelihood of this model given a data set "new". This works fine if it is used in the optim() routine with only one data set.  

NLL=function(coef)

{

      beta0=coef[1]

      mu=coef[2]

      k=coef[3]

      I=coef[4]             

      data_sim=Model_1(beta0,mu,k,I)  

      LLsum=0 #log likelihood sum       

      for (j in 1:length(data_sim[,1]))

            {

            if (data_sim[j,2]<0 || data_sim[j,3]<0)

            {return(1500)}  

            for (i in 1:length(new[,1]))

            {

                  if (new[i,1]==data_sim[j,1] && is.na(new[i,1])==F &&
is.na(data_sim[j,1])==F)
                  # if age of real data = age of simulated model and not
equal to NA
                  {

                        if (is.na(new[i,2])==F && is.na(data_sim[j,2])==F &&
is.na(data_sim[j,3])==F)
                        # if state of real data and state of simulated model
are not equal to NA
                        {

                             tmp1=new[i,2]*log(data_sim[j,2])

                             tmp2=(1-new[i,2])*log(data_sim[j,3])

 

                             if (tmp1=="-Inf" || tmp1=="NaN"){tmp1=-700}

                             if (tmp2=="-Inf" || tmp2=="NaN"){tmp2=-700}

                             LLsum=LLsum+tmp1+tmp2

                        }

                  }

                  

            }

      }

 

return(-LLsum)

}
         


Simon Ruegg

Dr.med.vet., PhD student

Institute for Parasitology

Winterthurstr. 266a

8057 Zurich

Switzerland  

phone: +41 44 635 85 93

fax: +41 44 635 89 07

e-mail: s.ruegg@access.unizh.ch  

        [[alternative HTML version deleted]]



R-help@stat.math.ethz.ch 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. Received on Sat Aug 26 22:51:00 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 27 Aug 2006 - 04:22:44 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.