# Comparing data to model and returning a chi squared value

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:

#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()


### 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():
# ...

# ...

commands = {
1: test_linear_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,
# ...
}

# 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()

• 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. Sep 11 '14 at 20:11
• 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! Sep 11 '14 at 20:18
• 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 Sep 11 '14 at 20:26
• 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? Sep 11 '14 at 20:38
• 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 Sep 12 '14 at 10:49