Looking at BIC

gistfile1.txt

import numpy as np

NumPoints = 100

def computeFit(numPoints, order):
    xvals = np.arange(numPoints)*3.
    ygiven = xvals*1.5 + 100.
    vecs = []
    for i in range(order+1):
        vecs.append(xvals**i)

    vecs = np.asarray(vecs).T

    result = np.linalg.lstsq(vecs, ygiven)
    fitCoeffs = result[0]

    yfit = np.zeros(len(xvals))
    for i in range(order+1):
        yfit += fitCoeffs[i]*xvals**i
    chi2 = np.sum((ygiven-yfit)**2)
    BIC = chi2+vecs.size*np.log(numPoints)

    return result, xvals, ygiven, yfit, chi2, BIC

print "numPoints  order  chi2   BIC"
for numPoints in 1e3, 1e4, 1e5:
    doBreak = False # set true when residuals first > 10, then solve once more
    for order in range(1, 10):
        result, xvals, ygiven, yfit, chi2,BIC = computeFit(numPoints, order)
        yresid = yfit - ygiven
        yresidRms = np.sqrt(np.mean(np.square(yresid)))
        #print "%6i  %6i   %10.2f    %s    %.2f    %.2f" % (numPoints, order, yresidRms, ", ".join("%6.2g" % (val,) for val in result[0]), chi2, BIC)
        print "%6i  %6i  %10.2f  %10.2f" % (numPoints, order, chi2, BIC)
        if doBreak:
            break
        elif yresidRms > 10:
            doBreak = True
    print

outputs:

numPoints order chi2 BIC
1000 1 0.00 13815.51
1000 2 0.00 20723.27
1000 3 0.00 27631.02
1000 4 0.03 34538.80
1000 5 404857.56 446304.10
1000 6 5484088.56 5532442.85

10000 1 0.00 184206.81
10000 2 0.00 276310.21
10000 3 0.03 368413.65
10000 4 68959528.56 69420045.58
10000 5 5949378116.54 5949930736.96

100000 1 0.00 2302585.09
100000 2 0.00 3453877.64
100000 3 63483467.20 68088637.39
100000 4 16964906007298.82 16964911763761.55