[SciPy-User] efficiency of the simplex routine: R (optim) vs scipy.optimize.fmin

Thu May 9 16:31:25 EDT 2013

On Thu, May 9, 2013 at 4:20 PM, Arnaldo Russo <arnaldorusso at gmail.com> wrote:
> Hi,
>
> Sorry for crossposting, but I'm trying to replace my "optim" function inside
> R, converting it to Python using scipy.optimize.fmin and leastsq.
> In Python, my results are aproximately  the same, when the function is
> setted to run with the right parameters.
> I have read that fmin needs as a result a single value, so I have to np.sum
> the return value and with leastsq I don't have to.
> But, when compared with R results of the same function ( "optim" function),
> the results are much different.
> Example below:
>
> # Inside python
> import numpy as np
> from scipy.optimize import leastsq, fmin
>
> def pp_min(params, x, y):
>     alpha = params[0]
>     Beta = params[1]
>     gama  = params[2]
>
>     model =  ( (y-gama*(1-np.exp(-alpha*x/gama))*np.exp(-Beta*x/gama))**2)
>
>     return model
>
>
> params = (0.4, 1.5 , 80)
> x = np.array([0, 1, 8, 16, 64, 164, 264, 364, 464, 564, 664, 764, 864, 964,
> 1064, 1164])
> y = np.array([0., 0.2436, 1.6128, 2.8224, 6.72, 15.1536, 22.176, 30.576,
> 31.1808, 40.2696, 47.4096, 41.7144, 61.6896, 56.6832, 62.5632, 63.5544])
>
> pp = leastsq(pp_min, params, args=(x, y))

for leastsq the function should not be squared, taking square and sum
is part of the algorithm.
drop the **2 in pp_min and it should work

from the docstring:

x = arg min(sum(func(y)**2,axis=0))
         y

Josef

>
> In [2]: pp
> Out[2]: (array([  8.48990490e-02,  -1.56537197e-02,   6.51505458e+01]), 1)
>
>
> # inside R
> x <- c(0, 1, 8, 16, 64, 164, 264, 364, 464, 564, 664, 764, 864, 964, 1064,
> 1164)
>
> y <- c(0.0, 0.2436, 1.6128, 2.8224, 6.72, 15.1536, 22.176, 30.576, 31.1808,
> 40.2696, 47.4096, 41.7144, 61.6896, 56.6832, 62.5632, 63.5544)
>
> dat<-as.data.frame(cbind(x, y))
> names(dat)<-c("x","y")
> attach(dat)
>
> pp_min<-function(params,data=dat)
> {  alpha<-params[1]
>   Beta<-params[2]
>   gama<-params[3]
>   return( sum( (y-gama*(1-exp(-alpha*x/gama))*exp(-Beta*x/gama))^2))
>  }
>
> pp <- optim(par=c(0.4, 1.5 , 80), fn=pp_min)
>
> Out:
>> pp
> $par
> [1]   0.09157204   0.02129695 148.89173924
>
> The third argument is specially much different.
> Any reasonable explanation?
>
> Thank you.
>
>
> ---
> Arnaldo D'Amaral Pereira Granja Russo
> Lab. de Estudos dos Oceanos e Clima
> Instituto de Oceanografia - FURG
>
>
>
>
> 2012/7/23 servant mathieu <servant.mathieu at gmail.com>
>>
>> Hi Denis,
>> Thanks for your response. For the fmin function in scipy, I took the
>> default ftol and stol values. I'm just trying to minize a chi square between
>> observed experimental data and simulated data. I've done this in python and
>> R with the Nelder-Mead algorithm, with exactly the same starting values.
>> While the solutions produced by R and python are not very different, R
>> systematicaly produces a lower chi-square after the same amount of
>> iterations. This may be related to ftol and stol, but I don't know which
>> value I should give to these parameters....
>> Cheers,
>> Mat
>>
>> 2012/7/20 denis <denis-bz-gg at t-online.de>
>>>
>>> Hi Mathieu,
>>>   (months later) two differences among implementations of Nelder-Mead:
>>> 1) the start simplex: x0 +- what ? It's common to take x0 + a fixed
>>> (user-specified) stepsize in each dimension. NLOpt takes a "walking
>>> simplex", don't know what R does
>>>
>>> 2) termination: what ftol, xtol did you specify ? NLOpt looks at
>>> fhi - flo: fhi changes at each iteration, flo is sticky.
>>>
>>> Could you post a testcase similar to yours ?
>>> That would sure be helpful.
>>>
>>> cheers
>>>    -- denis
>>>
>>>
>>> On 24/05/2012 10:15, servant mathieu wrote:
>>> > Dear scipy users,
>>> > Again a question about optimization.
>>> >   I've just compared the efficiency of the simplex routine in R
>>> > (optim) vs scipy (fmin), when minimizing a chi-square. fmin is faster
>>> > than optim, but appears to be less efficient. In R, the value of the
>>> > function is always minimized step by step (there are of course some
>>> > exceptions) while there is lot of fluctuations in python. Given that
>>> > the
>>> > underlying simplex algorithm is supposed to be the same, which
>>> > mechanism
>>> > is responsible for this difference? Is it possible to constrain fmin so
>>> > it could be more rigorous?
>>> > Cheers,
>>> > Mathieu
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>