From a given data set, I set out to complete a task which is below
- Fit the data of the previous exercise to fit Eq. (8.18) using the SciPy function
scipy.optimize.curve_fit()
. Plot the data as symbols and the fit as a line on linear and on log-log axes in two separate plots in the same figure window. Compare your results to those of the previous exercise.
Here is the data I had.
I attempted the problem with the following Python script.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(K, t, p):
return (K * np.power(t, p))
def func1(c, m, x):
return (m * x + c)
def LineFitWt(x, y, dy):
"""Fit to straight line.
Inputs: x, y, and dy (y-uncertainty) arrays.
Ouputs: slope and y-intercept of best fit to data.
"""
dy2 = dy ** 2
norm = (1. / dy2).sum()
xhat = (x / dy2).sum() / norm
yhat = (y / dy2).sum() / norm
slope = ((x - xhat) * y / dy2).sum() / ((x - xhat) * x / dy2).sum()
yint = yhat - slope * xhat
dy2_slope = 1. / ((x - xhat) * x / dy2).sum()
dy2_yint = dy2_slope * (x * x / dy2).sum() / norm
return slope, yint, np.sqrt(dy2_slope), np.sqrt(dy2_yint)
def redchisq(x, y, dy, slope, yint):
chisq = (((y - yint - slope * x) / dy) ** 2).sum()
return chisq / float(x.size - 2)
# Extract data from the text file for t, s and ds
t, s, ds = np.loadtxt("growthdata.txt", skiprows=3, unpack=True)
print ("time =", t)
print("size = ", s)
print ("Uncertainty =", ds)
# Set X and Y for the relevant axis
X = np.log(t)
Y = np.log(s)
dY = np.log(ds)
# Call LineFitWt function to calculate the gradient(m), y-intercept(c) and uncertainties in these(dm, dc)
m, c, dm, dc = LineFitWt(X, Y, dY)
rchisq = redchisq(X, Y, dY, m, c)
SciPy_Fit = curve_fit(func, t, s, p0=([0.54, 5.83]))
# Assign values for scipy fit
p_s = SciPy_Fit[0][0]
K_s = np.exp(SciPy_Fit[0][1])
# Calculate straight line properties for scipy parameters
m_s = p_s
c_s = np.log(K_s)
print ("m = ", m)
print ("c = ", c)
print ("dm = ",dm)
print ("dc = ", dc)
print ("Reduced chi square = ", rchisq)
print ("SciPy_Fit values = ", SciPy_Fit)
print ("SciPy Fit p = ", p_s)
print ("SciPy Fit K = ", K_s)
# Calculate the values for p and K
p = m
K = np.exp(c)
print ("p = ", p)
print ("K = ", K)
# Calculate values for custom fit points (y = mx + c)
Xext = 0.05*(X.max()-X.min())
Xfit = (np.array([X.min()-Xext, X.max()+Xext]))
Yfit= (c+m*Xfit)
# Calculate points for log log graph using Custom fit
Y_custom = func(K, t, p)
# Calculate values for SciPy y = mx + c fit
X_scipy = X
Y_scipy = (m_s * X + c_s)
# Calculate values for log-log plot using SciPy parameters
Y_scipy_loglog = func(K_s, t, p_s)
# Assign a figure object to plot on.
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(X, Y,"x", label="Data")
plt.errorbar(X, Y, yerr=dY, zorder=-1, label="Unc in s")
plt.plot(Xfit, Yfit,"+--", label="Custom Fit")
plt.plot(X_scipy, Y_scipy, "D--", label="SciPy Fit")
plt.text(-1, -2, 'Custom fit m={0:0.4f}, c={1:0.4f}'.format(m,c))
plt.text(-1, -5, 'SciPy Fit m={0:0.4f}, c={1:0.4f}'.format(m_s,c_s))
plt.xlabel("ln s")
plt.ylabel("ln r")
plt.legend()
plt.plot()
plt.subplot(2,1,2)
plt.loglog(t, s, "x", label="Data")
plt.errorbar(t, s, yerr=ds, label="Uncertainty in s", zorder=-1)
plt.loglog(t, Y_custom, "--", label="Custom fit")
plt.loglog(t, Y_scipy_loglog, "D--", label="SciPy Fit")
plt.xlabel("t")
plt.ylabel("s")
plt.text(0.1, 1, 'Custom fit K={0:0.4f}, p={1:0.4f}'.format(K, p))
plt.text(0.1, 10, 'SciPy Fit K={0:0.4f}, p={1:0.4f}'.format(K_s, p_s))
plt.legend()
plt.tight_layout()
plt.show()
I then obtained the following graphs
Here are the printed outputs.
time = [ 0.12 0.18 0.42 0.9 2.1 6. 18. 42. ]
size = [ 115. 130. 202. 335. 510. 890. 1700. 2600.]
Uncertainty = [10. 12. 14. 18. 20. 30. 40. 50.]
m = 0.5419106669494728
c = 5.837596806432137
dm = 0.5518341454711854
dc = 1.0280703524098318
Reduced chi square = 0.0002705479986781181
SciPy_Fit values = (array([ 0.67373601, -10.71638416]), array([[-3.48433902e+11, 1.40195650e+13], [ 1.40195650e+13, -5.64090353e+14]]))
SciPy Fit p = 0.6737360069383264
SciPy Fit K = 2.217856767946089e-05
p = 0.5419106669494728
K = 342.9541642911326
What am I missing here? Why is the SciPy fit so far from the data and the custom fit?