npar = 3

sigma = 0.5

data = 0.5 + 4.*np.cos(theta) + 1.0*np.sin(theta) + 0.5*r.randn(len(theta))

theta = np.linspace(0.,2.*np.pi,51)

nobs = len(theta)

M = np.column_stack([np.ones(nobs), np.cos(theta), np.sin(theta)])

MTM = np.dot(M.transpose(),M)
MTMINV = np.linalg.inv(MTM)
MTY = np.dot(M.transpose(),data)

P = np.dot(MTMINV,MTY)

C = MTMINV * sigma**2

Residuals = data -	np.dot(M,P)
ChiSq =	np.dot(Residuals.transpose(),Residuals)
RMS = m.sqrt(ChiSq/nobs)
print('RMS = %f'%(RMS))

print('a = %f +/- %f ' % (P[0],m.sqrt(C[0,0])))
print('b = %f +/- %f ' % (P[1],m.sqrt(C[1,1])))
print('c = %f +/- %f ' % (P[2],m.sqrt(C[2,2])))

pl.ion()
pl.plot(theta,data,'o')
xx = np.linspace(0.,2.*np.pi,51)
ModelY = P[0] + P[1]*np.cos(xx) + P[2]*np.sin(xx)
pl.plot(xx,ModelY,'k')
pl.xlabel('Theta')
pl.ylabel('Data')
pl.title('Big Boi Curve')
pl.grid('on')