loop in optim

classic Classic list List threaded Threaded
14 messages Options
Reply | Threaded
Open this post in threaded view
|

loop in optim

EdBo
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










Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

Joshua Wiley-2
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

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

Re: loop in optim

Joshua Wiley-2
Do you want something like this?

## Your data
afull <- read.table(textConnection("
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"), header = TRUE)
closeAllConnections()

## likelihood 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)

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")[[1]]
}

## Yields
> out
             al_j       au_j       sigma_j      b_j
Run:1  0.04001776 0.06010743  1.092618e-24 1.049971
Run:2  0.04002135 0.06008513 -7.156966e-25 1.049976
Run:3  0.04714390 0.27258724  3.303320e-24 0.948988
Run:4 -0.01000000 0.01000000  1.000000e-01 1.000000

On Mon, Jul 4, 2011 at 8:21 PM, EdBo <[hidden email]> wrote:

> 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
>
> --
> View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3645031.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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

EdBo
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
>





Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

Joshua Wiley-2
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

EdBo
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[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-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, click here.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

Berend Hasselman
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

EdBo
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[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)
>
>
> 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")[[1]]
+ }
> 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.html
Sent 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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

EdBo
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) )-
 
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

Joshua Wiley-2
In reply to this post by EdBo
I don't know what else to say.  Your code looks right to me, and it
all runs.  I would check the value of a at each loop:

for (i in 1:4) {
  a <- afull[seq(20 * (i - 1) +1, 20 * i), ]
  print(a) # so you can see what it is
  out[i, ] <- optim(llik, par = start.par, method = "Nelder-Mead")[[1]]
}

Also, I would try running the code in a clean version of R.  That is,
a version of R without any non-standard packages loaded, or clutter in
the workspace.  You can do this by shutting down R, deleting the old
workspace, and then starting a new session (there are many other ways
to do the same thing, that's just one).

Josh

On Wed, Jul 6, 2011 at 4:34 PM, EdBo <[hidden email]> wrote:

> 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[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)
>>
>>
>> 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")[[1]]
> + }
>> 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.html
> Sent 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-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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

Prof J C Nash (U30A)
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

EdBo
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))[[1]] 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[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

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))[[1]]
                }

## 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[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-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

Joshua Wiley-2
On Thu, Jul 7, 2011 at 9:45 PM, EdBo <[hidden email]> wrote:

> 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))[[1]] but it is giving me an error
> message.

You cannot assume that two different functions will produce the same
sort of output.  Take a look at:

tmp <- optimx(llik, par = c(.01, -.01, .1, -1), method =
"Nelder-Mead", control=list(maximize=TRUE))
str(tmp)

now what you tried to do is equivalent to:

tmp[[1]]

buy if you look at the structure of that

str(tmp[[1]])

its still a list, so...

tmp[[1]][[1]]

or thereabouts should work

>
> 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[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
>
> 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))[[1]]
>                }
>
> ## Yields
> out
>
>
> 2) Is it possible to control estimates of sigma_j so that they are
> never negative?

my knowledge of optimization approaches nil

>
> 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.

simple way would instead of saving your starting parameters in
start.par, add an extra row at the front of "out" to hold them.  Then
inside your loop, use something like:

for(i in 2:5)
out[i, ] <- optimx(llik, par = out[i -1, ]....

so that the previous rows estimates are the starting values.

>
> 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[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-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.
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3653222.html
> Sent 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-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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.