Python least squares with scipy.integrate.quad

Question

I;m having trouble decyphering an error message for my code to find some parameters to a complicated least squares fit in two parameters (eps and sig).

from pylab import *
import scipy
import numpy as np

from scipy import integrate, optimize

# Estimate parameters with least squares fit   
T = [90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300]
B = [-0.2221, -0.18276, -0.15348, -0.13088, -0.11293, -0.09836, -0.086301, -0.076166, -0.067535, -0.060101, -0.053636, -0.047963, -0.04295, -0.038488, -0.034494, -0.030899, -0.027648, -0.02469, -0.022, -0.019534, -0.017268, -0.015181]


def funeval(Temp,eps,sig):
    return -2.*np.pi*scipy.integrate.quad( lambda x: np.exp(4.*eps/Temp*((sig/x)**6.-(sig/x)**12.)*(x**2)) ,0.0,Inf )[0]

def residuals(p,y,Temp):
    eps,sig = p
    err = y-(funeval(Temp,eps,sig) )
    return err

print funeval(90.,0.001, 0.0002) 

plsq = scipy.optimize.leastsq(residuals, [0.00001, 0.0002], args=(B, T))

The funeval gives a reasonable float but when I run the code it returns:

error: Supplied function does not return a valid float.

The error doesnt seem sensitive to the initial conditions. I'm new to python so any help or guides to help would be much appreciated. Thanks.

Answer 1

With funeval(90.,0.001, 0.0002), Temp is a singular value; however, when you are calling scipy.optimize you are passing the entire T array to funeval causing scipy.integrate to crash.

A quick fix would be to do something like:

def funeval(Temp,eps,sig):
    out=[]
    for T in Temp:
        val = scipy.integrate.quad( lambda x: np.expm1( ((4.*eps)/T)* ((sig/x)**12.-(sig/x)**6.)* (x**2.) ), 0.0, np.inf )[0]
        out.append(val)
    return np.array(out)

def residuals(p,y,Temp):
    eps,sig = p
    err = y-(funeval(Temp,eps,sig) )
    return err

print funeval([90],0.001, 0.0002)

plsq = scipy.optimize.leastsq(residuals, [0.00001, 0.0002], args=(B, T))
(array([  3.52991175e-06,   9.04143361e-02]), 1)

This, unfortunately, does not converge very well. Can you explain what you are trying to do?