# loop in optim

14 messages
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[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(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
Open this post in threaded view
|

## Re: loop in optim

 Hi Edward, At least for me, your llik() function returns Inf for the starting values specified, so optim() never gets to estimate anything.  You need to alter llik() or find starting parameters that work before worrying about getting the for loop working. Cheers, Josh On Mon, Jul 4, 2011 at 2:34 AM, EdBo <[hidden email]> wrote: > 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[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(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 > > > > > > > > > > > > > -- > View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3643230.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 http://www.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.
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[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) out1 = optim(llik, par=start.par, method="Nelder-Mead") out1
Open this post in threaded view
|

## Re: loop in optim

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[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), ]    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 >
Open this post in threaded view
|

## 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[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. > > ______________________________________________ > [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.
Open this post in threaded view
|

## Re: loop in optim

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
Open this post in threaded view
|

## Re: loop in optim

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[1]; au_j=x[2]; sigma_j=x[3];  b_j=x[4]     sum(na.rm=T,         ifelse(a\$R_j< 0, log(1 / ( sqrt(2*pi) * sigma_j) )-
Open this post in threaded view
|

## Re: loop in optim

Open this post in threaded view
|

## Re: loop in optim

 In reply to this post by EdBo 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[1]; au_j=x[2]; sigma_j=x[3];  b_j=x[4] >     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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.