When to scale covariance by reduced chi square for error estimation

estimating_error_from_fit.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

x_values = np.array([-1.12163752, -1.01270045, -0.90458398, -0.79520288, -0.68504665,
        0.51738399,  0.62679271,  0.73598531,  0.84455173,  0.95383917,
       -1.12527093, -1.01615466, -0.90649584, -0.79919983, -0.68616358,
        0.51407783,  0.62394869,  0.73268909,  0.84366404,  0.9526824])
        
x_errors = np.array([ 0.00045032,  0.00045026,  0.00045033,  0.00045039,  0.00045036,
        0.00045031,  0.00045032,  0.00045027,  0.00045026,  0.00045025,
        0.00045019,  0.00045032,  0.00045028,  0.00045022,  0.00045025,
        0.00045027,  0.0004502 ,  0.00045015,  0.00045015,  0.00045014])
        
y_values = np.array([-23428.52636038, -19121.554785  , -15267.93044101, -11832.62588902,
        -8795.16114067,  -5037.41297316,  -7376.38054883, -10135.69508081,
       -13311.54188623, -16971.49358459, -23560.83904409, -19243.4881944 ,
       -15325.97613736, -11926.26280754,  -8796.70719322,  -4999.47538712,
        -7315.88408177, -10048.5216945 , -13279.96166291, -16942.68912419])
         
y_errors = np.array([ 0.64854273,  0.50100175,  0.45147562,  0.47937886,  0.46338331,
        0.5306017 ,  0.55723888,  0.47414199,  0.68312683,  0.53722473,
        0.53450656,  0.55133882,  0.47488376,  0.55920694,  0.55848547,
        0.66562946,  0.49291723,  0.54251019,  0.74092278,  0.50016698])



# The data is fitted to a parabola with an offset
# the curvature A is assumed to be known and the 
# offset in x is 0
# Curvature is fixed
A = -18.577e3

def model(x, A, c):
    return (A * x**2 + c)

# Function used in least square fit to find residuals
# since the x_errors are much larger then the y errors
# we include the x_errors by making a linear approximation
# to the function over the range of each poin + error bar
# so that the error in y due to the error in x is given
# by the gradient at that point
def residuals(p, x, y, y_error, x_error, A):
    c = p
    adj_error = np.sqrt(y_error**2 + (2 * A * x * x_error)**2)
    err = np.divide(y - model(x, A, c), adj_error)
    return err

def chi_2(p, x, y, y_error, x_error):
    return (residuals(p, x, y, y_error, x_error, A)**2).sum()

def red_chi_2(p, x, y, y_error, x_error):
    return chi_2(p, x, y, y_error, x_error) / (len(x_values) - len(p))


xerr_scale = [1, 1e2, 1e-2]

for scale in xerr_scale:
    scaled_x_errors = x_errors * scale
    #Intial guess for offset
    z0=[-70]
    
    #Fit the data using least square method
    results = leastsq(residuals, z0, 
                args=(x_values, y_values, y_errors, scaled_x_errors, A),
                full_output=1)
                
    redchi2 = red_chi_2(results[0], x_values, y_values, 
                                    y_errors, scaled_x_errors)
                
    err_lit = np.sqrt(results[1][0])
    err_scipy = np.sqrt(results[1][0] * redchi2)
                
                
    print "Errors scaled by:\t", scale
    print "Offset:\t\t\t", np.round(results[0][0], 4)
    print "Chi 2:\t\t\t", np.round(chi_2(results[0], x_values, y_values, y_errors,
                                scaled_x_errors), 4)
    print "Red. Chi 2:\t\t", np.round(redchi2, 4)
    print "Error literature\t", np.round(err_lit[0], 4)
    print "Error scipy\t\t", np.round(err_scipy[0], 4)
    print ""
    

                
plt.errorbar(x_values, y_values, yerr=y_errors, xerr=x_errors, fmt=".")
x_range = np.linspace(x_values.min(), x_values.max(), 100)
plt.plot(x_range, model(x_range, A, results[0]))
plt.show()