[SciPy-user] Using SciPy/NumPy optimization

[Brandon Nuttall]

Folks,

I've gotten further, but still am not there. I have the following code:

import numpy as np
from scipy.optimize import fmin_tnc

class hyp_func:
     def __init__(self,*args):
         self.x, self.y = args
         self.avg = y.mean()
     def rmsd(self,*args):
         """A function suitable for passing to the fmin() minimizers
         """
         a, b, c = args[0]
         sse = 0.0
         sst = 0.0
         # minimize the RMSD
         for i in range(len(x)):
             y = a*(1.0+b*c*x[i])**(-1.0/b)
             sse += (self.y[i]-y)**2
             sst += (self.avg-y)**2
         # this won't work because I don't know what x should be or if it 
is related
         diff = [0.0,0.0] # this is just an initial value and won't get any 
results
         # diff[0] = (-a*c-a*c*c*b*?????)**((-1-b)/b)
         # diff[1] = ((-1-b)/b)*(-a*c*c*b)*(-a*c-a*c*c*b*?????)**(((-1-b)/b)-1)
         return (sse / sst, diff)

# fake data for testing
a, b, c = [1.25,0.75,0.25]
x = []
y = []
for i in range(1,61):
     x.append(float(i))
     y.append(a*(1.0+b*c*float(i))**(-1.0/b))

x = np.array(x)
y = np.array(y)

myfunc=hyp_func(x,y)
params0 = [1.0,0.5,0.5]

retcode, nfeval, optimal = fmin_tnc(myfunc.rmsd,params0,
                                     bounds=[(None,None),(0.001,5.0), 
(-1.0,1.0)])

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'm not sure these derivatives are particularly informative about which way 
to go to minimize the RMSD. At the moment, the code fails with a 
SystemError (error return without exception set) in the call to 
moduleTNC.minimize(). So, any suggestions on my next step?

FWIW: In tnc.py, in test(), the call to fmin_tnc() has a typo in the 
keyword arguments. "maxnfeval=" should be "maxfun=". When running tnc.py to 
execute the tests, example() and test1fg() run. The test2fg() detects an 
infeasible condition (apparently, I guess that is what was to be tested); 
Python raises an exception and terminates without running the other tests.


Brandon C. Nuttall

BNUTTALL at UKY.EDU                         Kentucky Geological Survey
(859) 257-5500                                     University of Kentucky
(859) 257-1147 (fax)                              228 Mining & Mineral 
Resources Bldg
http://www.uky.edu/KGS/home.htm       Lexington, Kentucky 40506-0107