Reputation: 109
Update: solved! It is producing parameters with the correct signs now, and they do fit the curve. The problem was defining func(a,b,c,x) but curve_fit needs to read x first: func(x,a,b,c). Thanks everyone for all the help! I'll have quantitative analysis when I meet with my boss today :)
Here's some of the new fits: https://i.sstatic.net/8idj0.jpg
(I still get a run-time error:
RuntimeWarning: overflow encountered in power
return a*(math.e**(b*(math.e**(c*x))))
)
Can anyone help me figure out what's wrong with this code? I am new to scipy. I am trying to model bacterial growth with the Gompertz equation, but my code produces a curve_fit that's completely wrong. You can see images of my real data, the model equation, and the fit this code produces in this imgur album Thanks!
Fixed code:
#!/usr/bin/python
from numpy import *
from scipy.optimize import curve_fit
values = numpy.asarray(values)
y = values[:2000//5].astype(numpy.float)
y - y[0] #subtracting blank value
x = numpy.arange(len(y))*5
def Function(x,a,b,c):
#a = upper asymptote
#b = negative = x axis displacement
#c = negative = growth rate
return a*(math.e**(b*(math.e**(c*x))))
parameters, pcov = curve_fit(Function, x, y, p0=[0.1,-1300,-0.0077])
#Graph data and fit to compare
yaj = Function( numpy.asarray(x), parameters[0], parameters[1], parameters[2] )
figure(1, figsize=(8.5,11))
subplot(211)
plot(x,y,'g-')
xlim(min(x),max(x))
ylim(min(y),max(y))
subplot(212)
plot(x,yaj,'r-')
xlim(min(x),max(x))
ylim(min(yaj),max(yaj))
savefig('tempgraph.pdf')
return parameters
Upvotes: 3
Views: 5547
Reputation: 1070
Imports:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
Sample values:
values = np.array('0.400 0.400 0.397 0.395 0.396 0.394 0.392 0.390 0.395 0.393 0.392 0.392 0.390 0.388 0.390 0.388 0.385 0.383 0.388 0.387 0.387 0.387 0.385 0.386 0.387 0.379 0.379 0.378 0.375 0.376 0.374 0.373 0.372 0.368 0.373 0.370 0.371 0.370 0.370 0.370 0.367 0.368 0.368 0.365 0.365 0.366 0.364 0.361 0.361 0.356 0.355 0.357 0.354 0.353 0.350 0.351 0.353 0.355 0.350 0.354 0.352 0.351 0.348 0.348 0.347 0.345 0.346 0.343 0.348 0.346 0.344 0.343 0.342 0.341 0.346 0.346 0.345 0.343 0.348 0.345 0.346 0.342 0.344 0.344 0.340 0.341 0.345 0.345 0.343 0.339 0.343 0.344 0.343 0.346 0.344 0.344 0.345 0.347 0.344 0.344 0.338 0.340 0.343 0.340 0.342 0.336 0.334 0.336 0.337 0.338 0.338 0.343 0.342 0.342 0.336 0.334 0.336 0.330 0.325 0.324 0.323 0.319 0.323 0.322 0.318 0.314 0.314 0.319 0.315 0.316 0.313 0.315 0.314 0.314 0.315 0.313 0.308 0.312 0.311 0.310 0.312 0.311'
' 0.311 0.309 0.309 0.316 0.317 0.312 0.309 0.311 0.308 0.310 0.312'.split('\t'), dtype=float)
Old data preparation:
x=[]
y=[]
x_val = 0
for i in values: #values is a list of several thousand experimental data points
if x_val < 100:
x.append(float(x_val))
y.append(float(i))
x_val += 5
x = np.asarray(x)
y = np.asarray(y)
Easier data prep:
y1 = values[:100//5]
x1 = np.arange(len(y1))*5
Check it is the same:
print np.allclose(y, y1)
print np.allclose(x, x1)
Use numpy to define fit function:
def function(x, a,b,c):
#a = upper asymptote
#b = negative = x axis displacement
#c = negative = growth rate
return a*(np.exp(b*(np.exp(c*x))))
Fit using starting point p0:
pars, pcov = opt.curve_fit(function, x1, y1, p0=[0.1, -10, 0.1])
Draw:
yaj = function(x1, *pars)
plt.figure(1, figsize=(8.5, 11))
plt.plot(x1, y1, 'g-', x1, yaj, 'r-')
plt.xlim(min(x1), max(x1))
plt.ylim(min(y1), max(y1))
plt.show()
Upvotes: 3