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] # print (d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]) jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0) # print jd 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()
