-
Notifications
You must be signed in to change notification settings - Fork 70
/
example_uncertainty.py
40 lines (33 loc) · 1.52 KB
/
example_uncertainty.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import numpy as np
import pandas as pd
import pastas as ps
ps.set_log_level("ERROR")
# read observations and create the time series model
obs = pd.read_csv("data/head_nb1.csv", index_col=0,
parse_dates=True).squeeze("columns")
ml = ps.Model(obs, name="groundwater head")
# read weather data and create stressmodel
rain = pd.read_csv("data/rain_nb1.csv", index_col=0,
parse_dates=True).squeeze("columns")
evap = pd.read_csv("data/evap_nb1.csv", index_col=0,
parse_dates=True).squeeze("columns")
sm = ps.RechargeModel(prec=rain, evap=evap, rfunc=ps.Exponential,
name='recharge')
ml.add_stressmodel(sm)
# Solve
ml.solve()
#
df = ml.fit.prediction_interval()
inside = (obs > df.loc[obs.index, 0.025]) & (obs < df.loc[obs.index, 0.975])
print('percentage inside:', np.count_nonzero(inside) / len(inside) * 100)
# # Plot some results
axs = ml.plots.results(tmin="2010", tmax="2015", figsize=(10, 6))
axs[0].fill_between(df.index, df.iloc[:, 0], df.iloc[:, 1], color="gray",
zorder=-1, alpha=0.5, label="95% Prediction interval")
axs[0].legend(ncol=3)
df = ml.fit.ci_contribution("recharge", tmin="2010", tmax="2015")
axs[2].fill_between(df.index, df.iloc[:, 0], df.iloc[:, 1], color="gray",
zorder=-1, alpha=0.5, label="95% confidence")
df = ml.fit.ci_step_response("recharge", alpha=0.05, n=1000)
axs[3].fill_between(df.index, df.iloc[:, 0], df.iloc[:, 1], color="gray",
zorder=-1, alpha=0.5, label="95% confidence")