[SciPy-user] Using SciPy/NumPy optimization
Alok Singhal
as8ca@virginia....
Wed Feb 28 10:50:24 CST 2007
On 28/02/07: 10:56, Brandon Nuttall wrote:
> myfunc.rmsd() now returns the root mean square deviation which is to be
> minimized. However, in looking at the example code in
> scipy/optimize/tnc.py, I find that myfunc.rmsd() needs to return a second
> argument, g, the gradient of the function. It looks like this needs to be
> the first and second derivative of my function which if my [really, really]
> rusty calculus should be:
>
> y' = (-a*c-a*c*c*b*x)**((-1-b)/b)
> y'' = ((-1-b)/b)*(-a*c*c*b)*(-a*c-a*c*c*b*x)**(((-1-b)/b)-1)
>
I think the second return value should be a list of (first, partial)
derivatives with respect to the parameters being estimated, i.e., a,
b, c in your case. The example in optimize/tnc.py defines the
function as pow(x[0],2.0)+pow(abs(x[1]),3.0), and the derivatives g[0]
and g[1] are defined as 2.0*x[0] and
sign(x[1])*3.0*pow(abs(x[1]),2.0), which are derivatives of f with
respect to x[0] and x[1], the parameters. I calculated those
derivatives (and I think did it right :-) ), you can see them in an
earlier message by me in this thread.
> I'm not sure these derivatives are particularly informative about which way
> to go to minimize the RMSD.
Yeah, that is why it makes more sense to calculate the derivatives wrt
the parameters being estimated.
Here is another way to do the fitting (adapted from scipy cookbook):
import numpy as np
from scipy import rand
from scipy.optimize import leastsq
import pylab
n = 151
x = np.mgrid[1:10:n*1j]
def f(p, x):
return p[0]*(1.0+p[1]*p[2]*x)**(-1.0/p[1])
def errf(p, x, y):
return f(p, x) - y
abc = [15.0, 2.5, 0.3]
y = f(abc, x) + rand(n) - 0.5
abc_guess = [30.0, 4.0, 1.0]
abc1, success = leastsq(errf, abc_guess[:], args = (x, y))
print success # should be 1
print abc1 # I get [ 14.71976079 2.4536336 0.27560377]
pylab.plot(x, y, x, f(abc1, x))
pylab.show()
-Alok
--
Alok Singhal * *
Graduate Student, dept. of Astronomy * * *
University of Virginia
http://www.astro.virginia.edu/~as8ca/ * *
More information about the SciPy-user
mailing list