[R] optim() problem

[nashjc at uottawa.ca]

> I am trying to fit a exponential power distribution
>
> y = b/(2*pi*a^2*gamma(2/b))*exp(-(x/a)^b)
>
> to a bunch of data for x and y I have in a table.
>> data
>            x         y
> 1         25        27
> 2         75        59
> 3        125       219
> ...
> 259    12925         1
> 260    12975         0
>
> I know "optim" should do a minimisation, therefor I used as the
> optimisation function
>
> opt.power <- function(val, x, y) {
>    a <- val[1];
>    b <- val[2];
>    sum(y - b/(2*pi*a^2*gamma(2/b))*exp(-(x/a)^b));
> }
>
> I call: (with xm and ym the data from the table)
>
> a1 <- c(0.2, 100)
> opt <- optim(a1, opt.power, method="BFGS", x=xm, y=ym)
>
> but no optimisation of the parameter in a1 takes place.
> Any ideas?

As someone who had a big hand in creating BFGS, I'm guessing (and guessing
can be dangerous) that there is a bad scaling somewhere that simply
prevents progress because the numerical derivatives are not easy to
calculate in such situations. If you send me data off-line, I'll take a
look.

Have you checked that your function evaluates properly? Always a good idea.

For information, I'm working on trying to improve optim's capabilities and
user-friendliness to try to address some of the issues that arise. I
welcome input. In particular, Ravi Varadhan and I would like to develop a
method "GUIDED" to allow users a way to set up the calls so they are
sensible for their problems.

John Nash