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:
print 'Testing quadratic model', '\n'

#Scipy Optimise cuve_fit Model Produces Expected Values Using A Quadratic Model Polynomial Model

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

StackExchange Code Review Q#62661, answer score: 4