I want to adjust the parameters of a model to a given set of data.
I'm trying to use scipy
's function curve_fit
in Sage, but I keep getting
TypeError: only length-1 arrays can be converted to Python scalars
Here´s my code:
from numpy import cos,exp,pi
f = lambda x: exp( - 1 / cos(x) )
import numpy as np
def ang(time): return (time-12)*pi/12
def temp(x,maxtemp):
cte=(273+maxtemp)/f(0)**(1/4)
if 6<x and x<18:
return float(cte*f(ang(x))**(1/4)-273)
else:
return -273
lT=list(np.linspace(15,40,1+24*2))
lT=[float(num) for num in lT] #list of y data
ltimes=np.linspace(0,24,6*24+1)[1:]
ltimes=list(ltimes) #list of x data
u0=lT[0]
def u(time,maxtemp,k): #the function I want to fit to the data
def integ(t): return k*exp(k*t)*temp(t,maxtemp)
return exp(-k*time)*( numerical_integral(integ, 0, time)[0] + u0 )
import scipy.optimize as optimization
print optimization.curve_fit(u, ltimes, lT,[1000,0.0003])
scipy.optimize.curve_fit
expects the model function to be vectorized: that is, it must be able to receive an array (ndarray, to be precise), and return an array of values. You can see the problem right away by adding a debug printout
def u(time,maxtemp,k):
print time % for debugging
def integ(t): return k*exp(k*t)*temp(t,maxtemp)
return exp(-k*time)*( numerical_integral(integ, 0, time)[0] + u0 )
The output of print will be the entire array ltimes
you are passing to curve_fit
. This is something numerical_integral
is not designed to handle. You need to give it values one by one.
Like this:
def u(time,maxtemp,k):
def integ(t): return k*exp(k*t)*temp(t,maxtemp)
return [exp(-k*time_i)*( numerical_integral(integ, 0, time_i)[0] + u0) for time_i in time]
This will take care of the “only length-1 arrays can be converted" error. You will then have another one, because your lists ltimes and lT are of different length, which doesn't make sense since lT is supposed to be the target outputs for inputs ltimes. You should revise the definitions of these arrays to figure out what size you want.