You are right Joshua.
I changed the code because I failed to understand how you attached the full
data set. How you made the data part of your code.
I am new to R so I am used to one way of attaching data(the way I redone it).
-----Original Message-----
From: "Joshua Wiley-2 [via R]" <ml-node+3647622-640541925-247...@n4.nabble.com>
Date: Tue, 5 Jul 2011 20:28:09
To: EdBo<n.bow...@gmail.com>
Subject: Re: loop in optim
On Tue, Jul 5, 2011 at 7:07 PM, EdBo <n.bow...@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@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/
______________________________________________
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.
_______________________________________________
If you reply to this email, your message will be added to the discussion below:
http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3647622.html
To unsubscribe from loop in optim, visit
http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3643230&code=bi5ib3dvcmFAZ21haWwuY29tfDM2NDMyMzB8LTEwOTc4OTY3Mw==
--
View this message in context:
http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3647640.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
______________________________________________
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.
