import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import *
from scipy import stats
import smooth
import glob

modelruns=1000.
modlen=500
calcslope=[]

for kk in np.arange(modelruns):
    dF=np.random.normal(size=(modlen+18))*20
    dF=smooth.smoothListGaussian(dF,degree=9)
    dR=np.random.normal(size=(modlen))*1
    lmd=6
    c=14*12. # heat capacity
    t=[0.];dtdt=[]

    for ii in range(modlen):
        dtdt1=(dF[ii]+dR[ii]-lmd*t[-1])/c
        t.append(t[-1]+dtdt1);dtdt.append(dtdt1)

    t=np.array(t);dtdt=np.array(dtdt)
    # t=(t[1:]+t[:-1])/2
    t=t[:-1]

    calcslope.append(stats.linregress(t,dR-lmd*t)[0])

print 'avg. bias = %4.2f, avg. std. = %4.2f' %
((-np.average(calcslope)/lmd-1)*100,(np.std(calcslope)/lmd)*100)