[R] loop in optim

[Joshua Wiley]

On Thu, Jul 7, 2011 at 9:45 PM, EdBo <n.bowora at gmail.com> 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]
> <ml-node+3651249-776550566-247474 at n4.nabble.com> 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]]
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
https://joshuawiley.com/