title: Fit an integral function with parametric limit to data with Python (Debye Model)
body: <p>I am trying to fit resistivity vs temperature data to Bloch-Gruneisen formula for resistivity in metals:</p>
  
  <p><a href="http://en.wikipedia.org/wiki/Electrical_resistivity_and_conductivity#Metals" rel="nofollow">function</a></p>
  
  <p>as you can see there is an integral function with a parametric limit. I don't know how to implement an algorithm to run a least squares fit.
  I came up with:</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  import pylab as pl
  import scipy as sp
  from scipy.optimize import leastsq
  
  #retrieve data from file
  data = pl.loadtxt('salita.txt')
  Temp = data[:, 0]
  Res = data[:, 2]
  
  def debye_func(p, T, r):
      rho0, AD, TD = p
      coeff = AD*np.power(T, 5)/np.power(TD, 4)
      f = np.power(x^5)/np.power(np.sinh(x), 2) #function to integrate
      err_debye = r - rho0 - coeff * #integral???
      return err_debye
  
  p0 = sp.array([0.0001 , 0.00001, 50])
  
  plsq = leastsq(debye_func, p0, args=(Temp, Res))
  
  print plsq
  </code></pre>
  
  <p>Ideas on how could I write it?</p>
  
  <p>EDIT: my code has become:</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  import pylab as pl
  import scipy as sp
  from scipy.optimize import leastsq
  from scipy.integrate import quad
  
  #retrieve data from file
  data = pl.loadtxt('salita.txt')
  Temp = data[:, 0]
  Res = data[:, 2]
  
  def debye_integrand(x):
      return np.power(x, 5)/np.power(np.sinh(x), 2)
  
  def debye_func(p, T, r):
      rho0, AD, TD = p
      coeff = AD*np.power(T, 5)/np.power(TD, 4)
      err_debye = r - rho0 - coeff * quad(debye_integrand, 0, TD/(2*T))
      return err_debye
  
  p0 = sp.array([0.0001 , 0.00001, 50])
  
  plsq = leastsq(debye_func, p0, args=(Temp, Res))
  
  print plsq
  </code></pre>
  
  <p>Now I get a ValueError:</p>
  
  <pre><code>Traceback (most recent call last):
    File "debye.py", line 24, in &lt;module&gt;
      plsq = leastsq(debye_func, p0, args=(Temp, Res))
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/minpack.py", line 348, in leastsq
      m = _check_func('leastsq', 'func', func, x0, args, n)[0]
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/minpack.py", line 14, in _check_func
      res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
    File "debye.py", line 19, in debye_func
      err_debye = r - rho0 - coeff * quad(debye_integrand, 0, TD/(2*T))
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 247, in quad
      retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 296, in _quad
      if (b != Inf and a != -Inf):
  ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
  </code></pre>
  
  <p>I think that means I'm providing leastsq an argument it can't take, but I don't know how to modify my code.</p>
  
  <p>EDIT2: I solved my equation analytically with Maxima and I got</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  import pylab as pl
  import scipy as sp
  from scipy.optimize import leastsq
  from scipy.integrate import quad
  from scipy.special import zetac
  from mpmath import polylog
  
  #retrieve data from file
  data = pl.loadtxt('salita.txt')
  Temp = data[:, 0]
  Res = data[:, 2]
  
  def debye_integrand(x):
      return np.power(x, 5)/np.power(np.sinh(x), 2)
  
  def debye_func(p, T, r, integral):
      rho0, AD, TD = p
      coeff = AD*np.power(T, 5)/np.power(TD, 4)
      den = np.exp(TD/T) -1
      m1 = 5*((TD/(2*T))**4)*np.log(np.exp(TD/(2*T)+1)*(np.exp(TD/T)-1)+120*polylog(5, np.exp(TD/(T))*(1-np.exp(TD/(2*T)))
      m2 = 120*(TD/(2*T))*polylog(4, np.exp(TD/(2*T)))*(np.exp(np.exp(TD/T))-1)+60*((TD/(2*T))**2)*polylog(3, np.exp(TD/(2*T))*(1-np.exp((TD/(2*T)))
      m3 = 20*((TD/(2*T))**3)*polylog(2, np.exp(TD/(2*T))*(np.exp(TD/T)-1)+120**polylog(5, -np.exp(TD/(2*T)))*(1-np.exp(TD/T))
      m4 = 120*(TD/(2*T))*polylog(4, -np.exp(TD/(2*T)))*(np.exp(TD/T)-1)+60*((TD/(2*T))**2)*polylog(3, -np.exp(TD/(2*T)))*(1-np.exp(TD/T))
      m5 = 20*((TD/(2*T))**3)*polylog(2, -np.exp(TD/(2*T)))*(np.exp(TD/T)-1) -2*((TD/(2*T))**5)*np.exp(TD/T)
      m6 = 5*((TD/(2*T))**4)*np.log(1-np.exp(TD/(2*T))*(np.exp(TD/T)-1)
      zeta = 15.0*zetac(5)/2
  
      integral = (m1+m2+m3+m4+m5+m6)/den +zeta
  
      err_debye = r - rho0 - coeff * integral
      return err_debye
  
  #initizalied with Einstein model fit
  p0 = sp.array([0.00001 , 0.0000001, 70.0])  
  
  plsq = leastsq(debye_func, p0, args=(Temp, Res))
  
  print plsq
  </code></pre>
  
  <p>It says <code>SyntaxError: invalid syntax</code> in m2. I tried to do it with loops in the numerical way, but I didn't succeed.</p>
  
  <p>My .txt file is <a href="http://pastebin.com/tvzcdxYA" rel="nofollow">here</a>, if you want to try. First column is temperature, third one is resistivity.</p>