[R] optim

[ripley@stats.ox.ac.uk]

Do read the help page.  It says:

     `fnscale' An overall scaling to be applied to the value of `fn'
          and `gr' during optimization. If negative, turns the problem
          into a maximization problem. Optimization is performed on
          `fn(par)/fnscale'.

     `parscale' A vector of scaling values for the parameters.
          Optimization is performed on `par/parscale' and these should
          be comparable in the sense that a unit change in any element
          produces about a unit change in the scaled value.

You have not used either, AFAICS.

Also, I doubt if the function value returned by calls to integrate is a 
smooth function of the parameters, so scaling is particularly important 
here, and you may need to supply ndeps too.

Attempting to optimize blindly without supplying derivatives is asking far 
too much of a computer program.


On Fri, 28 Feb 2003, Remigijus Lapinskas wrote:

> Dear all,
> 
> I have a function MYFUN which depends on 3 positive parameters TETA[1],
> TETA[2], and TETA[3]; x belongs to [0,1].
> I integrate the function over [0,0.1], [0.1,0.2] and
> [0.2,0.3] and want to choose the three parameters so that
> these three integrals are as close to, resp., 2300, 4600 and 5800 as
> possible. As I have three equations with three unknowns, I expect the
> exact fit, i.e., the SS (see below) to be zero. However, the optim
> function never gives me what I expect, the minimal SS value(=res$value)
> never comes close to zero, the estimates of the parameters, res$par,
> wildly depends on init etc.
> I would be grateful for any comments on this miserable situation.
> 
> aa <- c(2300,4600,5800)
> init <- c(2.5,8000,0.84) # initial values of parameters
> print(init)
> ###################
> myfun <- function(x,TETA) TETA[2]*(((1-x)^(-1/TETA[3]))-
> 1)^(1/TETA[1])
> ###################
> x <- seq(0,0.3,by=0.01)
> plot(x,myfun(x,init),type="l")
> ###################
> LSS <- function(teta,aa)
> {
> integr <- numeric(3)
>    for(i in 1:3)
>    {integr[i] <- 10*integrate(myfun,
>    lower=(i-1)/10,upper=i/10,TETA=teta)$value
>    }
> SS <- sum((integr-aa)^2) # SS=Sum of Squares
> SS
> }
> ####################
> res <- optim(init,LSS,aa=aa,
> method = "L-BFGS-B",lower=c(0,0,0.5))
> print(res$par)
> print(res$value)
> 
> 
> > source("C:/Program Files/R/integral.R")
> [1]    2.50      7000.00         0.84        # initial
> [1]    2.3487221 6999.9999823    0.5623628   # final
> [1] 75613.05                                 # minSS
> > source("C:/Program Files/R/integral.R")
> [1]     2.5      15000            0.84       # initial
> [1]     2.125804 14999.999747     2.241179   # final
> [1] 50066.35                                 # minSS
> 
> 
> 
> Best regards,
> Remigijus                          mailto:remigijus.lapinskas at maf.vu.lt
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> http://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595