loop in optim

Previous Topic Next Topic
 
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
CONTENTS DELETED
The author has deleted this message.
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
CONTENTS DELETED
The author has deleted this message.
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
CONTENTS DELETED
The author has deleted this message.
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
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: loop in optim

Prof J C Nash (U30A)
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

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
CONTENTS DELETED
The author has deleted this message.