Least Squares Estimation in Python

print  'calculate trend'
db1 = Database("DAHITI_Waterlevels")
db1.connect()
parameter = ["date","height","std"]
tablename = "tid_%06d" % (int(dahiti_id))
data = db1.select(tablename,parameter,orderby=["date ASC"])				
print 'Points',len(data)
A = np.ones((len(data),2))
b = np.ones((len(data),1))			
for i in range(0,len(data)):
    d = data[i]
    jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0)
    b[i] = float(d['height'])
    A[i][1] = jd
db1.disconnect()

A = np.matrix(A)
b = np.matrix(b)
x = ((A.T * A).I * A.T) * b
v = (A * x - b)
sigma2 = float(((v.T * v) / (b.size - 2)))
Qxx = (A.T*A).I
Sigma_xx = sigma2 * Qxx

error = np.sqrt(float(Sigma_xx[1,1])) * 365.25 * 1000.0
trend = float(x[1])*365.25*1000.0

start_ts = julianDay2000((data[0]['date'][0:4]),(data[0]['date'][5:7]),(data[0]['date'][8:10]),12,0,0)
end_ts = julianDay2000((data[-1]['date'][0:4]),(data[-1]['date'][5:7]),(data[-1]['date'][8:10]),12,0,0)
years = (end_ts - start_ts) / 365.25	

print "Sea Level Trend: %.3f +- %.3f mm/year in %.1f years" % (trend,error,years)
#~ print "   Min. SLA: %.3f mm (%s)" % (start_sla,data[0]['date'])
#~ print "   Max. SLA: %.3f mm (%s)" % (end_sla,data[-1]['date'])

print path+"/www_timeseries_trend.dat"
output = open(path+"/www_timeseries_trend.dat","w")
output.write(julianDayDate(start_ts)['date'][0:10]+' %.3f\n' % (float((x[1] * start_ts + x[0]))))
output.write(julianDayDate(end_ts)['date'][0:10]+' %.3f\n' % (float((x[1] * end_ts + x[0]))))
output.close()

trend =  "Trend: %.1f \261 %.1f mm/year (%.1f years)" % (trend,error,years)
output = open('/tmp/www_timeseries_trend.legend','w')
output.write('G -0.025i\n')
output.write('N 1\n')			
output.write('S 0.12c c 0.1c darkgreen 0.25p 0.3c '+trend+'\n')
output.close()