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 numpyimport pylabimport matplotlib.pyplotimport scipy.optimizefrom 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 Dataxdata = 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 Valuestest=int(raw_input(Enter A Number Assigned Theoretical Model To Test: ))#Theoretical Models To Be Testedif 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 Modelchi_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 Valuesmatplotlib.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()
Comparing data to model and returning a chi squared value
python;graph;statistics;physics
PEP8You say you want to follow PEP8, but you're pretty far from it. A couple of examples:#Recorded Observed Dataxdata = 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 Dataxdata = 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 flowYour 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, poptAnd then you can call the selected command like this:func, popt = commands[test]()So that the chi_squared = line will work.Code organizationYou 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()