SciPy LeastSq Failed to Provide Best Solution

Question

I tried to use Leastsq to fit a very simple curve. However, its solutions are not optimized. Could anyone give me some suggestion? Below is my code:

from scipy import optimize
import numpy as np

hl_obs = np.array([10.0, 23.0, 20.0])
ph=np.array([5.0,7.0,9.0])

tp=60

def residuals(k_abn, ph, tp, hl_obs):
        hr=np.log(2)/(hl_obs*24.0)
        ph_adj=6013.79/(tp+273.15) + 23.6521*np.log10(tp+273.15)-64.7013
        err = peval(k_abn, ph, ph_adj)-hr
        return err

def peval(k_abn, ph, ph_adj):
    temp= k_abn[0]*np.power(10,-ph) + k_abn[1] + k_abn[2]*np.power(10,(-ph_adj + ph))
    return temp

k_abn =np.array([1, 0, 0])

from scipy.optimize import leastsq
p,ier = leastsq(residuals, k_abn, args=(hl_obs, ph, tp), maxfev=2000000)
print p, ier

From EXCEL's solver, I know the solution should be k_abn=[165, 0.001237578, 2.14]. Once I fed Excel's solution to the function peval, it generated the right answer...

peval([165,0.001238,2.14], 5.0, 13.01573)=0.002888113
peval([165,0.001238,2.14], 7.0, 13.01573)=0.001255701
peval([165,0.001238,2.14], 9.0, 13.01573)=0.001444057

In addition, I tried to increase the precision using epsfcn=np.finfo(np.float32).eps Can anyone give me some suggestions? Thanks!

Answer 1

If you get args in the right order:

p,ier = leastsq(residuals, k_abn, args=(ph, tp, hl_obs), maxfev=2000000)

You get the result:

[  1.65096852e+02   1.23712405e-03   2.14392540e+00]