nobs = len(theta)
M = np.column_stack((np.ones(nobs),np.cos(theta),np.sin(theta)))
MTM = np.dot(np.transpose(M),M)
MTMINV = np.linalg.inv(MTM)
MTY = np.dot(np.transpose(M),data)
P = np.dot(MTMINV,MTY)
print('Solution = ',P)

Residuals = data - np.dot(M,P)
ChiSq = np.dot(Residuals.transpose(),Residuals)
RMS = math.sqrt(ChiSq/nobs)

print('RMS = %f'%(RMS))

C = MTMINV * ChiSq/(nobs-npar)

# print solution
print('Constant = %f +/- %f ' % (P[0],math.sqrt(C[0,0])))
print('Linear   = %f +/- %f ' % (P[1],math.sqrt(C[1,1])))
print('Quadratic= %f +/- %f ' % (P[2],math.sqrt(C[2,2])))