<jupyter_start><jupyter_text>Example linear regression model using Bayesian inferenceA linear regression model is defined as$$y_t \sim N(\beta_0 + \beta_1 x, \sigma^2)$$where $y_t$ is the modelled quantity of interest (TSS), x is the predictor variable (NTU), $\beta_0$, $\beta_1$, and $\sigma$ are unknown quantities that we wish to infer. $N(\mu,\sigma^2)$ represents a normal distribution with a mean, $\mu$ and standard deviation, $\sigma$. The linear regression model is the mean process. The inference is performed using an MCMC sampler with `pymc3`.PyMC3 uses Bayesian statistics to find distributions on model parameters i.e.,$$p(parameters | data) \propto p(data | parameters) \times p(parameters)$$In words this says the probability of the parameters given the data (the posterior distribution) is proportional to the probability of the data given the model parameters (the likelihood) times the probability of the parameters (the prior belief).<jupyter_code>import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import pickle
import seaborn as sns

from theano import shared
%matplotlib notebook
# Load the data
infile = './In-situ_LISST_2031_RS2019.pkl'
with open(infile,'rb') as f:
    data = pickle.load(f)
data.head()
# Scatter plot of the data
x = data['Bottle_conc']
y = data['LISST_conc']

plt.figure()
plt.plot(x,y,'o')
# Use seaborn to make a fancy "stats" plot of the data
g = sns.jointplot('Bottle_conc', 'LISST_conc', data=data, kind="reg",
                  xlim=(0, 50), ylim=(0, 200), color="c", height=7)


# Convert the import data from numpy to theano.tensor types
yin = shared(y.values)
xin = shared(x.values)

with pm.Model() as model:
    
    # Define the priors on the parameters i.e., p(parameters)
    beta_0 = pm.TruncatedNormal('beta_0',mu=0,sigma=1,lower=0) # Intercept with lower bound
    beta_1 = pm.Normal('beta_1',3,3) # Slope 
    sigma = pm.HalfNormal('sigma',5)
    
    # DEfine the mean function (linear regression)
    mu = beta_0 + beta_1*xin
    #mu = beta_1*xin # No intercept model

    
    # Define the model for the likelihood i.e., p(data | parameters) (this is where we input the data)
    y_t = pm.Normal('y_t',mu, sigma, observed=yin)
    #y_t = pm.Gamma('y_t',mu=mu, sd=sigma, observed=yin)

    # Do the inference
    trace = pm.sample(500, tune=1000)

#plots on the left show the posterior distributions of each unknown parameter
# Plots on the right are trace plots - they should jump around when the mcmc sampling is working
pm.traceplot(trace)
# Summary prints out the details about each parameter
# n_eff and Rhat indicate how well the mcmc went
pm.summary(trace)
pm.plot_posterior(trace)
# Get some estimates of new data
xnew = np.linspace(0,50,100)

# Set the value of xin to do a prediction on new values 
xin.set_value(xnew)
with model:
    ppc = pm.sample_posterior_predictive(trace,samples=500)
ppc['y_t'].shape
# Plot the mean and 99% credible intervals of the prediction
y_mean = np.mean(ppc['y_t'],axis=0)
y_min = np.percentile(ppc['y_t'],0.5,axis=0)
y_max = np.percentile(ppc['y_t'],99.5,axis=0)

plt.figure()
plt.plot(xnew,y_mean,'k--')
plt.fill_between(xnew,y_min,y2=y_max,color='0.5')
plt.plot(x,y,'ro')
plt.xlabel('Bottle')
plt.ylabel('LISST')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.115255273616631, 'fraction_numerical': 0.020483032711708957, 'mean_word_length': 4.612349914236707, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '218857', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1152, 'n_tokens_neox': 1091, 'n_words': 373, 'path': 'bayesian_inference_example.ipynb', 'repo_name': 'IOMRC/notebooks', 'size': '3271'}