<jupyter_start><jupyter_text>Subsurface Data Analytics  Bayesian Linear Regression for Subsurface Data Analytics in Python  <NAME>, Associate Professor, University of Texas at Austin  [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)  PGE 383 Exercise: Bayesian Linear Regression for Subsurface Modeling in Python Here's a simple workflow, demonstration of Bayesian linear regression for subsurface modeling workflows. This should help you get started with building subsurface models with data analytics and machine learning. Here's some basic details about Bayesian linear regression.  Bayesian UpdatingThe prediction for Bayesian linear regression is distributed:\begin{equation}y \sim N(\beta^{T}X, \sigma^{2} I)\end{equation}We estimate the model parameter distributions through Bayesian updating for infering the model parameters from a prior and likelihood from training data.\begin{equation}p(\beta | y, X) = \frac{p(y,X| \beta) p(\beta)}{p(y,X)}\end{equation}In general for continuous features, we must use a sampling method, such as Markov chain Monte Carlo to sample the posterior. Workflow GoalsLearn the basics of isotonic regression in Python for analysis, modeling and prediction of porosity from density. This includes:* Basic Python workflows and data preparation* Training / fitting a ridge regression model* Checking the model and learning about the impact of hyperparameters Objective Content to support the PGE 383: Subsurface Machine Learning class.   Getting StartedHere's the steps to get setup in Python with the GeostatsPy package:1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). 2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. 3. In the terminal type: pip install geostatspy. 4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. There are examples below with these functions. You can go here to see a list of the available functions, https://git.io/fh4eX, other example workflows and source code.  Import Required PackagesLet's import the GeostatsPy package.<jupyter_code>import os                                                   # to set current working directory 
import numpy as np                                          # arrays and matrix math
import scipy.stats as st                                    # statistical methods
import pandas as pd                                         # DataFrames
import matplotlib.pyplot as plt                             # for plotting
from sklearn.isotonic import IsotonicRegression
from sklearn.model_selection import train_test_split        # train and test split
from sklearn.metrics import explained_variance_score
import seaborn as sns
seed = 73073<jupyter_output><empty_output><jupyter_text>If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  We will need to install a couple of packages that are not available in the standard Anaconda 3.7 build.* pymc3 - streamlines Markov chain Monte Carlo Sampling* arviz - exploratory analysis of Bayesian modelsThis can be accomplished with 'pip install [package_name]' in a terminal window.<jupyter_code>import warnings
warnings.filterwarnings('ignore')
import pymc3 as pm<jupyter_output><empty_output><jupyter_text>Set the working directoryI always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time).  Also, in this case make sure to place the required (see below) data file in this working directory.<jupyter_code>#os.chdir("d:\PGE383")                                       # set the working directory<jupyter_output><empty_output><jupyter_text>Loading DataLet's load the provided dataset. 'GrainSize_por.csv' is available at https://github.com/GeostatsGuy/GeoDataSets. It is a comma delimited file with 20 grain size ($mm$) and porosity (as a percentage) measures from the subsurface. We load the data file with the pandas 'read_csv' function into a data frame we called 'df'.  We preview each with the head function from Pandas DataFrames.<jupyter_code>#df = pd.read_csv("GrainSize_Por.csv")                    # read a .csv file in as a DataFrame
df = pd.read_csv("https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/GrainSize_Por.csv")
df.head()<jupyter_output><empty_output><jupyter_text>Visualize the DataLooks like a reasonable dataset for a linear regression model.<jupyter_code>X_data = df.iloc[:,0]
y_data = df.iloc[:,1]
plt.subplot(111)
plt.scatter(X_data, y_data, color='black', s = 20, alpha = 0.5, label='sample data')
plt.title('Porosity from Grainsize with Training Data'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Least Squares Linear Regression Let's calculate the ordinary least squares (OLS) linear regression model for comparison.* for a non-informative prior our Bayesian linear regression model parameters will be centered on the OLS solution.<jupyter_code>OLS_slope, OLS_intercept, r_value, p_value, std_err = st.linregress(X_data,y_data)

print('The model parameters are, slope (b1) = ' + str(round(OLS_slope,2)) + ', and the intercept (b0) = ' + str(round(OLS_intercept,2)))
plt.subplot(111)
plt.plot(X_data, y_data, 'o', label='sample data', color = 'red', alpha = 0.2, markeredgecolor = 'black')
plt.plot(X_data, OLS_intercept + OLS_slope*X_data, label='model', color = 'black')
plt.title('Porosity vs Grain Size')
plt.xlabel('Grain Size (mm)')
plt.ylabel('Porosity (%)')
plt.legend()
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>The following code was based on the work of <NAME> [blog post](https://towardsdatascience.com/introduction-to-bayesian-linear-regression-e66e60791ea7) and example workflow in [Jupyter Notebook](https://github.com/WillKoehrsen/Data-Analysis/blob/master/bayesian_lr/Bayesian%20Linear%20Regression%20Demonstration.ipynb).https://github.com/WillKoehrsen/Data-Analysis/blob/master/bayesian_lr/Bayesian%20Linear%20Regression%20Demonstration.ipynb Bayesian Linear RegressionIn the following code we take advantage of the pymc3 package's API to build a very readible and compact model.The parts include:* specify the prior models for $\beta_0$, $\beta_1$ and $\sigma^2$.* specify the model estimates* define the distribution of the predictions and the observated predictor features* specify the specific form for the MCMC steps* conduct MCMC samplingNote the default MCMC method in pymc3 is the Hamiltonian Monte Carlo (HMC).* HMC reduces correlation between succesive samples by including proposals for moves to distant states, $L$* No U-Turn Sampler (NUTS) is an extension of HMC that auotmatically tunes L to avoid oscillation<jupyter_code>with pm.Model() as linear_all_data:
    # prior distributions, Naive - set to Gaussian with very large variance
    intercept = pm.Normal('Intercept', mu = 0, sd = 10.0)# model B0 - Intercept Prior
    slope = pm.Normal('slope', mu = 0, sd = 10.0)        # model B1 - Slope Prior 
    sigma = pm.HalfNormal('sigma', sd = 10.0)            # homoscedastic error in standard deviation
    
    # specify the model, linear regression
    mean = intercept + slope * X_data                  # estimate of model
    
    # define the distribution of the predictions from the model and provide the sampled response features
    Y_obs = pm.Normal('Y_obs', mu = mean, sd = sigma, observed = y_data.values[:]) # Observed values
    
    # specify the sampler
    step = pm.NUTS(target_accept = 0.95) # use the No-U-Turn Sampler (NUTS) for Markov Chains Monte Carlo MCMC, addaptive step size

    # Posterior distribution
    linear_trace_all_data = pm.sample(draws = 500, tune = 500, step = step, random_seed = seed)
    #linear_trace_all_data = pm.sample(draws = 500, tune = 500, step = step, random_seed = seed,discard_tuned_samples = False)<jupyter_output>Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, slope, Intercept]
Sampling 4 chains, 0 divergences: 100%|██████████| 4000/4000 [00:03<00:00, 1166.87draws/s]<jupyter_text>Summary Statistics for the Model Parameter Posterior SamplesIt is useful to assess the mean, dispersion and credible intervals for the model parameters.<jupyter_code>def trace_mean(x):
    return pd.Series(np.average(x, 0), name='mean')

def trace_sd(x):
    return pd.Series(np.std(x, 0), name='sd')

def trace_975(x):
    return pd.Series(np.percentile(x, 97.5), name='97.5')

def trace_025(x):
    return pd.Series(np.percentile(x, 2.5), name='2.5')

# Get the expectation / average for the model parameters from sampling the posterior distribution
bayes_exp_intercept = pm.stats.summary(linear_trace_all_data, stat_funcs = [trace_mean]).iloc[0,0]
bayes_exp_slope = pm.stats.summary(linear_trace_all_data, stat_funcs = [trace_mean]).iloc[1,0]

pm.stats.summary(linear_trace_all_data, stat_funcs = [trace_mean, trace_sd, trace_025, trace_975])
pm.stats.summary(linear_trace_all_data, stat_funcs = [trace_mean])<jupyter_output><empty_output><jupyter_text>Diagnostic PlotsFor each of the model parameters, let's look at the results for the 4 Markov chains.* distributions / PDFs based on kernel density estimates (KDE) from the 1000 MCMC states* trace of the samples over each state to observe for burn-in, stuck in local minimums / sufficient variation<jupyter_code>pm.traceplot(linear_trace_all_data)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=2.0, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Posterior DistributionsWe can take our 1000 states from our 4 Markov chains as samples of the posterior distributions of our Bayesian linear regresion model parameters. * once again the distributions are represented as PDFs from kernel density estimation<jupyter_code>pm.plot_posterior(linear_trace_all_data, credible_interval = 0.95)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Visualize the Credible Intervals for the Model ParametersWe can visualize the credible intervals for all of the model parameters* one for each Markov chain<jupyter_code>pm.forestplot(linear_trace_all_data, credible_interval = 0.95)
plt.show()<jupyter_output><empty_output><jupyter_text>Sampling the Bayesian Regression ModelWe can sample from the resulting uncertainty distributions for the model parameters.* we'll sample 1000 times and visualize the reuslting models<jupyter_code>plt.subplot(111)
plt.scatter(X_data, y_data, color='black', s = 20, alpha = 0.5, label='sample data', zorder = 3)
pm.plot_posterior_predictive_glm(linear_trace_all_data, samples = 1000, eval=np.linspace(2, 100, 100), linewidth = 1, color = 'red', alpha = 0.02, label = 'Bayesian Posterior Models',lm = lambda x, sample: sample['Intercept'] + sample['slope'] * x, zorder = 1)
plt.plot(X_data, bayes_exp_intercept + bayes_exp_slope*X_data, label='OLS Model', color = 'black', zorder= 2)
plt.title('Porosity from Grainsize Bayesian Model'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Posterior Prediction Let's pick a specific grain size and apply the above 1000 samples of the model to sample the Bayes posterior prediction.* note we would still need to add in the noise term, $\sigma$, to the model samples<jupyter_code>grain_size = 40
nruns = linear_trace_all_data['slope'].shape[0]

plt.subplot(121)
plt.scatter(X_data, y_data, color='black', s = 20, alpha = 0.5, label='sample data', zorder = 3)
pm.plot_posterior_predictive_glm(linear_trace_all_data, samples = 1000, eval=np.linspace(2, 100, 100), linewidth = 1, color = 'red', alpha = 0.02, label = 'Bayesian Posterior Models',lm = lambda x, sample: sample['Intercept'] + sample['slope'] * x, zorder = 1)

plt.plot(X_data, OLS_intercept + OLS_slope*X_data, label='model', color = 'black', zorder= 2)
plt.vlines(x = grain_size,ymin = 0, ymax = 30, label = 'OLS Prediction',colors = 'black', linestyles='--')
plt.title('Porosity from Grainsize Bayesian Model'); plt.xlabel('Grain Size (mm)'); plt.ylabel('Porosity (%)')
plt.xlim(0,100)#; plt.ylim(0,1500000)
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()

plt.subplot(122)
model_uncert = linear_trace_all_data['Intercept'] + linear_trace_all_data['slope'] * grain_size
sns.kdeplot(model_uncert, label = 'Bayes Model Uncertainty', color = 'blue')
bayes_prediction = linear_trace_all_data['Intercept'] + linear_trace_all_data['slope'] * grain_size + np.random.normal(loc=0,scale=linear_trace_all_data['sigma'],size=nruns)
sns.kdeplot(bayes_prediction, label = 'Bayes Posterior Prediction', color = 'black')
plt.vlines(x = OLS_intercept + OLS_slope * grain_size,ymin = 0, ymax = 2.5, label = 'OLS Prediction',colors = 'red', linestyles='--')
plt.xlim(0,30),plt.ylim(0,1.0)
plt.xlabel('Porosity (%)'), plt.ylabel('Density');
plt.title('Posterior Prediction for Grain Size = ' + str(grain_size));
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.0, wspace=0.2, hspace=0.2)
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10557515878616797, 'fraction_numerical': 0.0273817925194072, 'mean_word_length': 5.142609449501517, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 54, 'https://': 13, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '43981', 'license': 'MIT', 'n_tokens_mistral': 4818, 'n_tokens_neox': 4412, 'n_words': 1400, 'path': 'SubsurfaceDataAnalytics_BayesianRegression.ipynb', 'repo_name': 'ahmeduncc/PythonNumericalDemos', 'size': '14184'}