Comparing data to model and returning a chi squared value

Question

This is quite basic but useful to test various (9) different models using one set of data. I have tried to make it clear and use the PEP8 formatting.

I am currently creating a version that can read CSV files to make it even faster. I'd like any suggestions on how I could improve my code or any hints as to what I could add next.

import numpy
import pylab
import matplotlib.pyplot
import scipy.optimize
from scipy.optimize import curve_fit
''' A Program That Determines The Reduced Chi Squared Value For Various Theoretical Models.'''
'''The Best Fit Parameters Are Derived Using Levenberg-Marquardt Algorithm Which Solves The Non-Linear Least Squares Problem.'''
''' Test = 1 Assumes A Linear Model '''
''' Test = 2 Assumes A Quadratic Model '''
''' Test = 3 Assumes A Cubic Model '''
''' Test = 4 Assumes A Quartic Model '''
''' Test = 5 Assumes A Quintic Model '''
''' Test = 6 Assumes A Logarithmic Model '''
''' Test = 7 Assumes A Power Law Model '''
''' Test = 8 Assumes A Exponential Model '''
''' Test = 9 Assumes A Exponential Offset Model '''
#Recorded Observed Data
xdata = numpy.array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
ydata = numpy.array([-10, -7, 1, 5, 7, 7, 6, 7, 10, 15])
xerror = numpy.array([0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3])
yerror = numpy.array([0.1, 0.1, 0.4, 0.3, 0.2, 0.1, 0.1, 0.2, 0.1, 0.1])
#User Prompted Input Values
test=int(raw_input("Enter A Number Assigned Theoretical Model To Test: "))
#Theoretical Models To Be Tested
if test == 1:
 print 'Testing linear model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Linear Model
 def func(x, a, b):
 return a*x + b
 #Constants Of Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant Ax is',("%.2f" %popt[0])
 print 'Constant B is',("%.2f" %popt[1]), '\n'
if test == 2:
 print 'Testing quadratic model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Quadratic Model Polynomial Model
 def func(x, a, b, c):
 return a*x**2 + b*x + c
 #Constants Of Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant Ax^2 is',("%.2f" %popt[0])
 print 'Constant Bx is',("%.2f" %popt[1])
 print 'Constant C is',("%.2f" %popt[2]), '\n' 
if test == 3:
 print 'Testing cubic model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Cubic Model Polynomial Model
 def func(x, a, b, c, d):
 return a*x**3 + b*x**2 + c*x + d
 #Constants Of Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant Ax^3 is',("%.2f" %popt[0])
 print 'Constant Bx^2 is',("%.2f" %popt[1])
 print 'Constant Cx is',("%.2f" %popt[2])
 print 'Constant D is',("%.2f" %popt[3]), '\n'
if test == 4:
 print 'Testing quartic model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Quartic Polynomial Model
 def func(x, a, b, c, d, e):
 return a*x**4 + b*x**3 + c*x**2 + d*x + e
 #Constants Of Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant Ax^4 is',("%.2f" %popt[0])
 print 'Constant Bx^3 is',("%.2f" %popt[1])
 print 'Constant Cx^2 is',("%.2f" %popt[2])
 print 'Constant Dx is',("%.2f" %popt[3])
 print 'Constant E is',("%.2f" %popt[4]), '\n'
if test == 5:
 print 'Testing quintic model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Quintic Polynomial Model
 def func(x, a, b, c, d, e, f):
 return a*x**5 + b*x**4 + c*x**3 + d*x**2 + e*x + f
 #Constants Of Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant Ax^5 is',("%.2f" %popt[0])
 print 'Constant Bx^4 is',("%.2f" %popt[1])
 print 'Constant Cx^3 is',("%.2f" %popt[2])
 print 'Constant Dx^2 is',("%.2f" %popt[3])
 print 'Constant Ex is',("%.2f" %popt[4])
 print 'Constant F is',("%.2f" %popt[5]), '\n'
if test == 6:
 print 'Testing logarithmic model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Logarithmic Model
 def func(x, a, b):
 return a*numpy.log(x)+ b 
 #Constants Of Least Squared Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant A is',("%.2f" %popt[0])
 print 'Constant Bx is',("%.2f" %popt[1]), '\n'
if test == 7:
 print 'Testing power law model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Power Law Model
 def func(x, a, b):
 return a*(x**b) 
 #Constants Of Least Squared Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant A is',("%.2f" %popt[0])
 print 'Constant B is',("%.2f" %popt[1]), '\n'
if test == 8:
 print 'Testing exponential model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Exponential Model
 def func(x, a, b, c):
 return a*numpy.exp(b*x + c)
 #Constants Of Least Squared Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant A is',("%.2f" %popt[0])
 print 'Constant Bx is',("%.2f" %popt[1])
 print 'Constant C is',("%.2f" %popt[2]), '\n'
if test == 9:
 print 'Testing exponetial offset model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Exponential Offset Model
 def func(x, a, b, c, d):
 return a*numpy.exp(b*x + c) + d 
 #Constants Of Least Squared Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant A is',("%.2f" %popt[0])
 print 'Constant Bx is',("%.2f" %popt[1])
 print 'Constant C is',("%.2f" %popt[2])
 print 'Constant D is',("%.2f" %popt[3]), '\n'
#Derived Chi Squared Value For This Model
chi_squared = numpy.sum(((func(xdata, *popt)-ydata)/xerror)**2)
reduced_chi_squared = (chi_squared)/(len(xdata)-len(popt))
print 'The degrees of freedom for this test is', len(xdata)-len(popt)
print 'The chi squared value is: ',("%.2f" %chi_squared)
print 'The reduced chi squared value is: ',("%.2f" %reduced_chi_squared)
#Observed Values Are Plotted With Expected Values
matplotlib.pyplot.figure()
matplotlib.pyplot.scatter(xdata, ydata, s=0)
matplotlib.pyplot.plot(xdata,func(xdata, *popt), label='Theoretical Model')
matplotlib.pyplot.errorbar(xdata, ydata, xerr=xerror, yerr=yerror, linestyle="", color='red')
matplotlib.pyplot.xlabel(' Observed Values For x ')
matplotlib.pyplot.ylabel(' f(x)')
matplotlib.pyplot.legend(loc='lower right')
matplotlib.pyplot.show()

Answer 1

PEP8

You say you want to follow PEP8, but you're pretty far from it. A couple of examples:

#Recorded Observed Data
xdata = numpy.array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
test=int(raw_input("Enter A Number Assigned Theoretical Model To Test: "))
print 'Constant Ax is',("%.2f" %popt[0])
reduced_chi_squared = (chi_squared)/(len(xdata)-len(popt))

This is how it should be, look carefully to find the differences on each line:

# Recorded Observed Data
xdata = numpy.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
test = int(raw_input("Enter A Number Assigned Theoretical Model To Test: "))
print 'Constant Ax is', ("%.2f" % popt[0])
reduced_chi_squared = chi_squared / (len(xdata) - len(popt))

You have any of these kind of issues. You can use the pep8 command line tool to point out all the violations. Install with pip install pep8.

Logical flow

Your many if test == 1, if test == 2 statements are mutually exclusive, so you really should use elif:

if test == 1:
 # ...
elif test == 2:
 # ...
# ...

Even better, put the content of each of those if blocks into functions, and then make a dictionary of those functions, for example:

def test_linear_model():
 # ...
def test_quadratic_model():
 # ...
commands = {
 1: test_linear_model,
 2: test_quadratic_model,
 # ...
}

So that you can run the right commands with:

commands[test]()

To use the func and popt defined in each function, return them! For example:

def test_linear_model():
 print 'Testing linear model', '\n'
 #Scipy Optimise cuve_fit Model Produces Expected Values Using A Linear Model
 def func(x, a, b):
 return a*x + b
 #Constants Of Theoretical Model
 popt, pcov = curve_fit(func, xdata, ydata)
 print 'Constant Ax is',("%.2f" %popt[0])
 print 'Constant B is',("%.2f" %popt[1]), '\n'
 return func, popt

And then you can call the selected command like this:

func, popt = commands[test]()

So that the chi_squared = line will work.

Code organization

You should not just dump your code in the global namespace. You should have all functionality inside functions. The main method that triggers everything can be in a main function, which you can start like this:

def main():
 commands = {
 1: test_linear_model,
 2: test_quadratic_model,
 # ...
 }
 # User Prompted Input Values
 test = int(raw_input("Enter A Number Assigned Theoretical Model To Test: "))
 func, popt = commands[test]()
 #Derived Chi Squared Value For This Model
 chi_squared = numpy.sum(((func(xdata, *popt) - ydata) / xerror) ** 2)
 reduced_chi_squared = chi_squared / (len(xdata) - len(popt))
 print 'The degrees of freedom for this test is', len(xdata) - len(popt)
 print 'The chi squared value is: ', ("%.2f" % chi_squared)
 print 'The reduced chi squared value is: ', ("%.2f" % reduced_chi_squared)
 # ... 
if __name__ == '__main__':
 main()

• user52931: Thanks for the feedback! I was looking at using the functions but I cannot see how I could incorporate a general version of the graph/chi squared calculations as I'm calling different functions.
• janos: Note sure if I understand correctly. I added to the Logical flow section, I hope it helps, let me know if you need something more!
• user52931: Sorry it was poorly worded. Right now chi/graph is defined outside of the functions. To calculate those for the functions they'd need to be defined after the function and repeated for each function. so def f(): Stuff chi = graph.. for each test
• janos: I noticed you need both func and popt. I revised again the middle part, and added one more example at the end. See what I mean?
• user52931: Hey, I tried to follow the organization you have but it would keep giving me indentation errors. I followed the codes logical flow but dont see why it doesnt like the indentation of main > test > function. I pasted my progress below, sorry if this sounds stupid. pastebin.com/02qHC4ra