# loop in optim Classic List Threaded 14 messages Reply | Threaded
Open this post in threaded view
|

## loop in optim

 Hi May you help me correct my loop function. I want optim to estimates al_j; au_j; sigma_j;  b_j by looking at 0 to 20, 21 to 40, 41 to 60 data points. The final result should have 4 columns of each of the estimates AND 4 rows of each of 0 to 20, 21 to 40, 41 to 60. ###MY code is n=20 runs=4 out=matrix(0,nrow=runs) llik = function(x)    {     al_j=x; au_j=x; sigma_j=x;  b_j=x     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(pnorm(au_j,mean=b_j*a\$R_m,sd=sqrt(sigma_j^2))-                            pnorm(au_j,mean=b_j*a\$R_m,sd=sqrt(sigma_j^2)))))                )    } start.par = c(0, 0, 0.01, 1) out1 = optim(llik, par=start.par, method="Nelder-Mead") for (i in 1: runs) {  index_start=20*(i-1)+1  index_end= 20*i  out[i]=out1[index_start:index_end] } out Thank you in advance Edward UCT ####My data R_j R_m -0.0625 0.002320654 0 -0.004642807 0.033333333 0.005936332 0.032258065 0.001060848 0 0.007114057 0.015625 0.005581558 0 0.002974794 0.015384615 0.004215271 0.060606061 0.005073116 0.028571429 -0.006001279 0 -0.002789594 0.013888889 0.00770633 0 0.000371663 0.02739726 -0.004224228 -0.04 0.008362539 0 -0.010951605 0 0.004682924 0.013888889 0.011839993 -0.01369863 0.004210383 -0.027777778 -0.04658949 0 0.00987272 -0.057142857 -0.062203157 -0.03030303 -0.119177639 0.09375 0.077054642 0 -0.022763619 -0.057142857 0.050408775 0 0.024706076 -0.03030303 0.004043701 0.0625 0.004951088 0 -0.005968731 0 -0.038292548 0 0.013381097 0.014705882 0.006424728 -0.014492754 -0.020115626 0 -0.004837891 -0.029411765 -0.022054654 0.03030303 0.008936428 0.044117647 8.16925E-05 0 -0.004827246 -0.042253521 0.004653096 -0.014705882 -0.004222151 0.029850746 0.000107267 -0.028985507 -0.001783206 0.029850746 -0.006372981 0.014492754 0.005492374 -0.028571429 -0.009005846 0 0.001031683 0.044117647 0.002800551
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 CONTENTS DELETED The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 Hi I have re-worked on my likelihood function and it is now working(#the code is below#).  May you help me correct my loop function.  I want optim to estimates al_j; au_j; sigma_j;  b_j by looking at 0 to 20,  21 to 40, 41 to 60 data points.  The final result should have 4 columns of each of the estimates AND 4 rows  of each of 0 to 20, 21 to 40, 41 to 60. #likelihood function a=read.table("D:/hope.txt",header=T) attach(a) a llik = function(x)    {     al_j=x; au_j=x; sigma_j=x;  b_j=x     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) out1 = optim(llik, par=start.par, method="Nelder-Mead") out1
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 CONTENTS DELETED The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 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. 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; au_j=x; sigma_j=x;  b_j=x      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), ]    out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[]  } 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 >
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 CONTENTS DELETED The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 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).From: "Joshua Wiley-2 [via R]" <[hidden email]> Date: Tue, 5 Jul 2011 20:28:09 -0700 (PDT)To: EdBo<[hidden email]>Subject: Re: loop in optim On Tue, Jul 5, 2011 at 7:07 PM, EdBo <[hidden email]> 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; au_j=x; sigma_j=x; Â b_j=x > Â  Â  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")[] > Â } > 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. > >______________________________________________ > [hidden email] 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/______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland 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, click here.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 EdBo wrote 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). You don't need to "attach" the data by using attach(). You read the data into an object afull and then select the part you need and store that in object a. BTW: shouldn't the for (i in 1:4) be for (i in 1:3) if I understand the original question correctly? Berend
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 I am sorry if I sound stupid but I am not able to correct the error even after running this code. > afull=read.table("D:/hope.txt",header=T) > llik = function(x) +   { +    al_j=x; au_j=x; sigma_j=x;  b_j=x +    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) > > > out <- matrix(NA, nrow = 4, ncol = 4, dimnames = list( +  paste("Run:", 1:4, sep = ''), +  c("al_j", "au_j", "sigma_j", "b_j"))) > > ## Estimate parameters based on rows 0-20, 21-40, 41-60 of 'afull' > for (i in 1:4) { +  a <- afull[seq(20 * (i - 1) +1, 20 * i), ] +  out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[] + } > out            al_j      au_j       sigma_j      b_j Run:1 0.1088116 0.1621605 -1.554167e-24 0.969153 Run:2 0.1088116 0.1621605 -1.554167e-24 0.969153 Run:3 0.1088116 0.1621605 -1.554167e-24 0.969153 Run:4 0.1088116 0.1621605 -1.554167e-24 0.959875 On 6 July 2011 11:46, Berend Hasselman [via R] <[hidden email]> wrote: > EdBo wrote: > 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). > > You don't need to "attach" the data by using attach(). > You read the data into an object afull and then select the part you need and > store that in object a. > > BTW: shouldn't the for (i in 1:4) be for (i in 1:3) if I understand the > original question correctly? > > Berend > > ________________________________ > If you reply to this email, your message will be added to the discussion > below: > http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3648171.html> To unsubscribe from loop in optim, click here. -- View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3650297.htmlSent from the R help mailing list archive at Nabble.com.         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 In reply to this post by EdBo I have one last theoretical question, I did not adjust my code prior so that it maximise the likehood function. I googled that to make optim maximise you multiply fn by -1. In my code, would that be the same as saying "-sum" on the "sum" part of my code (see below)? llik = function(x)    {     al_j=x; au_j=x; sigma_j=x;  b_j=x     sum(na.rm=T,         ifelse(a\$R_j< 0, log(1 / ( sqrt(2*pi) * sigma_j) )-
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 In reply to this post by EdBo CONTENTS DELETED The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 In reply to this post by EdBo CONTENTS DELETED The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 Hi Josh, 1) I got hold of the library(optimx) but surprising I seem to be failing to use it. On the previous code, that was working, I replaced the optim function with: out[i, ] <- optimx(llik, par = start.par, method = "Nelder-Mead", control=list(maximize=TRUE))[] but it is giving me an error message. Error in out[i, ] <- optimx(llik, par = start.par, method = "Nelder-Mead",  :   incorrect number of subscripts on matrix The full code is as follows: afull=read.table("D:/a.txt",header=T) afull library(optimx) #likilihood function llik = function(x)    {     al_j=x; au_j=x; sigma_j=x;  b_j=x     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 out <- matrix(NA, nrow = 4, ncol = 4, dimnames = list(paste("Qtr:", 1:4, sep = ''), c("al_j", "au_j", "sigma_j", "b_j"))) ## Estimate parameters based on rows 0-20, 21-40, 41-60 of a for (i in 1:4) {                 index_start=20*(i-1)+1                 index_end= 20*i   a=afull[index_start:index_end,]                 out[i, ] <- optimx(llik, par = start.par, method = "Nelder-Mead", control=list(maximize=TRUE))[]                 } ## Yields out 2) Is it possible to control estimates of sigma_j so that they are never negative? 3) I had not thought of using previous estimates as starting values of the next loop because I did not know I could do that but I think it will help quicken the maximisation. Edward UCT On 7 July 2011 13:40, John C Nash [via R] <[hidden email]> wrote: > 2 comments below. > > On 07/07/2011 06:00 AM, [hidden email] wrote: >> Date: Wed, 6 Jul 2011 20:39:19 -0700 (PDT) >> From: EdBo <[hidden email]> >> To: [hidden email] >> Subject: Re: [R] loop in optim >> Message-ID: <[hidden email]> >> Content-Type: text/plain; charset=us-ascii >> >> I have one last theoretical question, I did not adjust my code prior so >> that >> it maximise the likehood function. I googled that to make optim maximise >> you >> multiply fn by -1. >> >> In my code, would that be the same as saying "-sum" on the "sum" part of >> my >> code (see below)? >> >> llik = function(x) >>    { >>     al_j=x; au_j=x; sigma_j=x;  b_j=x >>     sum(na.rm=T, >> >>         ifelse(a\$R_j< 0, log(1 / ( sqrt(2*pi) * sigma_j) )- >> > The optimx package has a "maximize" control because we felt that the fnscale > approach, > while perfectly correct, is not comfortable for users and is not standard > across other > optimization tools. Note that this package is undergoing a fairly extensive > overhaul at > the moment (the development version is on R-forge in the project > 'optimizer') to include > some safeguards on functions that return NaN etc. as well as a number of > other changes -- > hopefully improvements. > > A second comment on this looping: Why do you not use the parameters from the > last > estimation as the starting parameters for the next? Unless you are expecting > very extreme > changes over the moving window of data, this should appreciably speed up the > optimization. > > John Nash > > ______________________________________________ > [hidden email] 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-tp3643230p3651249.html> To unsubscribe from loop in optim, click here.
Reply | Threaded
Open this post in threaded view
|

## Re: loop in optim

 CONTENTS DELETED The author has deleted this message.