[R] loop in optim

[Joshua Wiley]

On Tue, Jul 5, 2011 at 7:07 PM, EdBo <n.bowora at gmail.com> wrote:
> Hi Josh
>
> I have run the code and the structure of the output is what I wanted.
> However, the code is giving an identical result for all runs.

Right because the object "a" stays the same for all runs.

>
> I have attached the code I ran below as well as the output. I have just
> changed number of runs to match with the size of the data.
>
> a=read.table("D:/hope.txt",header=T)
> attach(a)
> a
>  #likilihood function
> llik = function(x)
>    {
>     al_j=x[1]; au_j=x[2]; sigma_j=x[3];  b_j=x[4]
>     sum(na.rm=T,
>         ifelse(a$R_j< 0, -log(1/(2*pi*(sigma_j^2)))-
>                            (1/(2*(sigma_j^2))*(a$R_j+al_j-b_j*a$R_m))^2,
>
>          ifelse(a$R_j>0 , -log(1/(2*pi*(sigma_j^2)))-
>                            (1/(2*(sigma_j^2))*(a$R_j+au_j-b_j*a$R_m))^2,
>
>                                        -log(ifelse (( pnorm (au_j, mean=b_j * a$R_m,
>                                                sd= sqrt(sigma_j^2))-
>                          pnorm(al_j, mean=b_j * a$R_m, sd=sqrt (sigma_j^2)
> )) > 0,
>                                (pnorm (au_j,mean=b_j * a$R_m, sd= sqrt(sigma_j^2))-
>                            pnorm(al_j, mean=b_j * a$R_m, sd=
> sqrt(sigma_j^2) )),
>                                1)) ))
>       )
>    }
> start.par = c(-0.01,0.01,0.1,1)
> #looping now
> runs=133/20+1 #total data points divided by number od days in each quater+1
> out <- matrix(NA, nrow = runs, ncol = 4,
>  dimnames = list(paste("Quater:", 1:runs, sep = ''),
>  c("al_j", "au_j", "sigma_j", "b_j")))
>
>  for (i in 1:runs) {
>   a[seq(20 * (i - 1) +1, 20 * i), ]

note that this is not what my original code does.  In my code, I
stored the full dataset in an object called "afull", then the object
"a" is assigned as a subect of the rows from afull.  Since the
likelihood function is coded to reference "a", as "a" changes, the
estimates change.  subsetting without assigning the output anywhere
does not actually change "a", so the likelihood function references
the full dataset.  Also, the way the funtion is written, there is no
need to attach "a", and this can be rather dangerous when you are
making changes because when you attach an object, a copy is created
but that is not updated with assignment, so for example:

dat <- data.frame(x = 1:10)
attach(dat)
# look at "x" from attached data frame
x
# now overwrite x in the data frame
dat$x <- rnorm(10)
# compare
dat$x
x

these are no longer the same even though you may be expecting them to
be the same.  To do that you would need to:

## remove copies
detach(dat)
## re-attach() so it now includes the updated version
attach(dat)
x

>   out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[[1]]
>  }
> out
> #results I am getting
>> out
>               al_j       au_j       sigma_j      b_j
> Quater:1 0.04001525 0.06006251 -7.171336e-25 1.049982
> Quater:2 0.04001525 0.06006251 -7.171336e-25 1.049982
> Quater:3 0.04001525 0.06006251 -7.171336e-25 1.049982
> Quater:4 0.04001525 0.06006251 -7.171336e-25 1.049982
> Quater:5 0.04001525 0.06006251 -7.171336e-25 1.049982
> Quater:6 0.04001525 0.06006251 -7.171336e-25 1.049982
> Quater:7 0.04001525 0.06006251 -7.171336e-25 1.049982
>>
>
>
>
>
>
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3647549.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
https://joshuawiley.com/