In [1]:
%matplotlib inline
import gzip
from pprint import pprint
import matplotlib.pyplot as plt
import numpy as np
import spacepy.toolbox as tb
import spacepy.plot as spp
import pandas as pd
In [2]:
dfile = 'Embudo_stream_ts_raw.txt.gz'
In [3]:
dat = pd.read_csv(dfile, skiprows=30, delimiter='\t', header=None, names=('agency', 'site','datetime', 'TZ',
'discharge', 'd_a', 'gague', 'g_a'))
print(dat[0:20])
In [4]:
index = pd.DatetimeIndex(dat['datetime'] + " " + dat['TZ'], infer_datetime_format=True)
dat.set_index(index, inplace=True)
dat.drop('datetime', axis=1, inplace=True)
dat.drop('TZ', axis=1, inplace=True)
print(dat[:10])
In [5]:
dat['discharge'].plot()
Out[5]:
In [6]:
dat['discharge'].plot()
plt.yscale('log')
In [7]:
daily_max = dat['discharge'].resample('1w',how='max')
dat['discharge'].plot()
daily_max.plot(label='Weekly Max')
plt.yscale('log')
plt.legend()
Out[7]:
In [56]:
print(daily_max.as_matrix().min(), daily_max.as_matrix().max())
daily_max.as_matrix()
data_mcmc = daily_max.as_matrix()[~np.isnan(daily_max.as_matrix())]
plt.hist(data_mcmc, 50)
Out[56]:
In [9]:
import pymc
In [10]:
# Define the priors for the location (xi), scale (alpha) and shape (kappa)
# parameters.
xi = pymc.Uniform('xi', rseed=True, lower=daily_max.min(), upper=daily_max.max(), doc='Location parameter')
@pymc.deterministic(plot=True)
def alpha(xi=xi):
"""Scale parameter"""
return 1./xi
kappa = pymc.Beta('kappa', rseed=True, alpha=5., beta=6., doc='Shape parameter')
# @pymc.data
# @pymc.stochastic_from_data('D', lower=0, upper=1e6, data=daily_max.as_matrix())
def D(value=data_mcmc, location=xi, scale=alpha, shape=kappa, lower=0, upper=1e6):
return pymc.gev_like(value, shape, location, scale)
In [11]:
@pymc.stochastic_from_data
In [12]:
type(daily_max.min())
Out[12]:
In [13]:
model = pymc.MCMC((xi, alpha, kappa, D))
In [14]:
model.sample(iter=50000, verbose=True, burn_till_tuned=True)
In [15]:
pprint(model.stats())
pymc.Matplot.plot(model)
In [16]:
model.variables
Out[16]:
In [68]:
# Define the priors for the location (xi), scale (alpha) and shape (kappa)
# parameters.
scale = pymc.Uniform('scale', lower=0, upper=10, doc='Scale')
shape = pymc.Uniform('shape', lower=0, upper=10, doc='Shape')
weibull = pymc.Weibull('weibull', observed=True, alpha=scale, beta=shape, doc='Weibull', value=data_mcmc)
w2 = pymc.Weibull('weibull', alpha=scale, beta=shape, doc='W2')
In [69]:
model = pymc.MCMC((scale, shape, weibull, w2))
In [70]:
model.sample(iter=100000, verbose=True, burn_till_tuned=True)
In [71]:
pymc.Matplot.plot(model)
In [72]:
pprint(model.stats())
In [83]:
predictive_traces = model.trace("weibull")[:]
plt.hist(predictive_traces, tb.logspace(1e3, 1e7, 100), normed=True)
plt.yscale('log')
plt.xscale('log')
In [93]:
hi,bi = np.histogram(predictive_traces, tb.logspace(1e3, 1e7, 100), normed=True)
In [94]:
import scipy.integrate
ans = []
for ii, b in enumerate(bi[1:],1):
ans.append(scipy.integrate.trapz(hi[:ii],bi[:ii],))
In [98]:
plt.plot(bi[1:], ans)
plt.xscale('log')
plt.yscale('log')
In [ ]: