import numpy as np
import scipy.optimize
from pymcmcstat.MCMC import MCMC
import matplotlib.pyplot as plt
# Initialize MCMC object
mcstat = MCMC()
# Next, create a data structure for the observations and control
# variables. Typically one could make a structure |data| that
# contains fields |xdata| and |ydata|.
ndp = 7
x = np.array([28, 55, 83, 110, 138, 225, 375]) # (mg / L COD)
x = x.reshape(ndp, 1) # enforce column vector
y = np.array([0.053, 0.060, 0.112, 0.105, 0.099, 0.122, 0.125]) # (1 / h)
y = y.reshape(ndp, 1) # enforce column vector
# data structure
mcstat.data.add_data_set(x, y)
def modelfun(x, theta):
return theta[0]*x/(theta[1] + x)
def ssfun(theta,data):
return sum((data.ydata[0] - modelfun(data.xdata[0], theta))**2)
# Calculate initial covariance matrix
def residuals(p, x, y):
return y - modelfun(x, p)
theta0, ssmin = scipy.optimize.leastsq(
residuals,
x0=[0.15, 100],
args=(mcstat.data.xdata[0].reshape(ndp,),
mcstat.data.ydata[0].reshape(ndp,)))
n = mcstat.data.n[0] # number of data points in model
p = len(theta0) # number of model parameters (dof)
ssmin = ssfun(theta0, mcstat.data) # calculate the sum-of-squares error
mse = ssmin/(n-p) # estimate for the error variance
J = np.array([[x/(theta0[1]+x)], [-theta0[0]*x/(theta0[1]+x)**2]])
J = J.transpose()
J = J.reshape(n, p)
tcov = np.linalg.inv(np.dot(J.transpose(), J))*mse
# Define model parameters, simulation options, and model settings.
mcstat.parameters.add_model_parameter(
name='$\mu_{max}$',
theta0=theta0[0],
minimum=0)
mcstat.parameters.add_model_parameter(
name='$K_x$',
theta0=theta0[1],
minimum=0)
mcstat.simulation_options.define_simulation_options(
nsimu=int(5.0e3),
updatesigma=1,
qcov=tcov)
mcstat.model_settings.define_model_settings(
sos_function=ssfun,
sigma2=0.01**2)
# Run simulation
mcstat.run_simulation()
# Extract results and print statistics
results = mcstat.simulation_results.results
names = results['names']
chain = results['chain']
s2chain = results['s2chain']
names = results['names'] # parameter names
mcstat.chainstats(chain, results)
Sampling these parameters: name start [ min, max] N( mu, sigma^2) $\mu_{max}$: 0.15 [ 0.00e+00, inf] N( 0.00e+00, inf) $K_x$: 49.05 [ 0.00e+00, inf] N( 0.00e+00, inf) [-----------------100%-----------------] 5000 of 5000 complete in 1.3 sec --------------------- name : mean std MC_err tau geweke $\mu_{max}$: 0.1606 0.0322 0.0027 28.6882 0.9409 $K_x$ : 72.5970 51.2930 4.4053 35.0766 0.7246 ---------------------
Define function for generating intervals, setup calculations, and generate.
def predmodelfun(data, theta):
return modelfun(data.xdata[0], theta)
mcstat.PI.setup_prediction_interval_calculation(
results=results,
data=mcstat.data,
modelfunction=predmodelfun)
mcstat.PI.generate_prediction_intervals(nsample=500,
calc_pred_int=False)
def format_plot():
plt.xlabel('x (mg/L COD)', fontsize=20)
plt.xticks(fontsize=20)
plt.ylabel('y (1/h)', fontsize=20)
plt.yticks(fontsize=20)
plt.title('Predictive envelopes of the model', fontsize=20);
Generating credible/prediction intervals: Interval generation complete
Available inputs: (Defaults in Parantheses)
plot_pred_int
: Flag to include PI on plot. (True
)adddata
: Flag to include data on plot. (False
)addlegend
: Flag to include legend on plot. (True
)figsizeinches
: Specify figure size in inches [Width, Height]. ([7,5]
)model_display
: Model display settings. (See below)data_display
: Data display settings. (See below)interval_display
: Interval display settings. (See below)Default display options:
interval_display = {'linestyle': ':', 'linewidth': 1, 'alpha': 0.5, 'edgecolor': 'k'}
model_display = {'linestyle': '-', 'marker': '', 'color': 'r', 'linewidth': 2, 'markersize': 5, 'label': 'model', 'alpha': 1.0}
data_display = {'linestyle': '', 'marker': '.', 'color': 'b', 'linewidth': 1, 'markersize': 5, 'label': 'data', 'alpha': 1.0}
mcstat.PI.plot_prediction_intervals()
format_plot()
data_display = dict(
marker='o',
color='k',
markersize=10)
mcstat.PI.plot_prediction_intervals(adddata=True, data_display=data_display)
format_plot()
model_display = dict(
linestyle='-.',
linewidth=3,
color='k',
marker='o',
markersize=10)
mcstat.PI.plot_prediction_intervals(adddata=True,
model_display=model_display)
format_plot()
interval_display = dict(
linestyle='-',
linewidth=3,
alpha=0.75,
edgecolor='k')
mcstat.PI.plot_prediction_intervals(adddata=True,
interval_display=interval_display)
format_plot()