Introduction to

Raul Maldonado, April 2021 BayPiggies

A little bit about me

• Data Analyst @ Autodesk, Teaching Fellow @ Delta Analytics
• You can find me at RaulingAverage.dev
• Enjoy Coffee, Learning, and Running..near the beach

Side remarks

This presentation does not reflect any workings or material from Autodesk. Also, I am not a core-contributor to the PyMC3 project, but a user.

Also, as we continue & recover from the great Panini, please consider helping or contributing to a social causes, nonprofits, or sponsoring Women & BIPOC applicants at your company. #PayItForward

What will we be talking about today?

PyMC3!

But before going into PyMC3, we will briefly cover Bayesian methodology & Bayesian Modeling.

Particularly:

• What is Bayesian Modeling?
• Brief Introduction to Probability
• Getting Started with PyMC3

What is Bayesian Modeling?

• Frequentist Statistics is a method of statistical inference in which conclusions from data is obtained by emphasizing the frequency or proportion of the data. The resulting information is a distribution for some event. The resulting calculations are a single value, point estimate.

• Bayesian Statistics is a method of statistical inference in which Bayes’ theorem is used to update the probability for a hypothesis as more evidence or information becomes available. The resulting calculations are a distribution for some event.

  Leveraging Bayesian Statistics, we can leverage probabilistics programming to create Bayesian (probabilistic) models to solve problems requiring estimated probabilities from resulting probability distributions.

Brief Introduction to Probability

A few core principles of Probabilities:

• are in the interval [0, 1]
• follow the product rule $p(A,B) = p(A|B)p(B)$
  • $p(A, B)$ expression represents the joint probability of A and B
  • $p(A|B)$ expression is used to indicate a conditional probability of A and b
• are Independent when $p(A) = p(B)$

Bayes' Theorem

$$p(\theta | y) = \tfrac{p(y|\theta)p(\theta)}{p(y)}$$

where

• $p(\theta)$: Prior (hypothesis)
• $p(y | \theta)$: Likelihood (data)
• $p(\theta |y)$: Posterior
• $p(y)$: Marginal Likelihood (evidence)

Coin-toss Example

Let's say we have a coin-toss scenario with no understanding of what expected probability.

We want to understand what is the probability that the coin-toss being heads will typically be.

There are several scenarios with a total number of trials $N$.

N= [0,1,2,3,4,8,16,32,50,150]

[0, 1, 2, 3, 4, 8, 16, 32, 50, 150]

Respectively for these $N$ trials, we have a prior probability $p(\theta)$ = 0.35, 35% chance of being heads.

And the likelihood $p(y|\theta)$ of that scenario as well.

Not yet leveraging PyMC3, we use scipy & matplotlib to respectively calculate & visualize our distributions of what the probability of coin-toss being heads would typically be, from an outputted distribution

Data:  [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
𝜃_real:  0.35

Text(0.5, 0.98, 'Coin Toss Scenarios for N Trials')

For the first subplot of 0 trials:

• The uniform (blue) prior represents all the possigble values for bias being equally probable a priori.
• The gaussian-like (orange) prior is centered and concentrated around 0.5
  • So this prior is compatible with information indicating that the coin has more or less about the same chance of landing heads or tails.
  • We could also say this prior is compatible with the belief that most coins are fair.
• The skewed (green) prior puts the most weight on the tail-based outcome

Already, this was enough to sufficiently get a gloss of probabilities. To learn more about Bayesian Statisitcs, we could continue.

Or, as textbook authors say,

"we leave that up to the reader"

Just kidding, we would like to proceed with an actual implemenation of the PyMC3 library.

The Main Takeway - Bayes Theorem

Bayes' Theorem

$$p(\theta | y) = \tfrac{p(y|\theta)p(\theta)}{p(y)}$$

where

• $p(\theta)$: Prior
• $p(y | \theta)$: Likelihood
• $p(\theta |y)$: Posterior
• $p(y)$: Marginal Likelihood (evidence)

Getting Started with PyMC3

Why hasn't Bayesian Modeling been an industry-wide recommendation?

Similar to the growth of Machine Learning & Artificial Intelligence, Probabilistic Programming, implementation of Bayesian Modeling, had a combination of innovation & limitations due to the following:

• computationally challenging
• conceptually challenging to work
• Algorithms advancements like the Metropolis or No-U-Turn Sampler are more recently relevant to for the advancement of the field itself

However, as data technology has become more advanced in the recent years, we can now leverage the opportunity to implement Probabilistic Programming.

What is PyMC3?

PyMC3 is a Python package for Bayesian statistical modeling and Probabilistic Machine Learning. It's associated to one of several Probabilistic Programming Languages (PPL) which focuses on advanced Markov chain Monte Carlo (MCMC) and variational fitting algorithms, and more.

Some features of PyMC3 are the following:

• Intuitive model specification syntax
• Leverage several (Sampling) Algorithms
  • Hamiltonian Monte Carlo Markov Chain
  • No U-Turn Sampler (NUTS)
  • Variational Inference

• Built on top of Theano

  • Theano features
  • Leverage dynamic C compilation for optimized computation
  • Automatic Differentiation

    Note: As of 2021, Theano is now "Aesara". See "The Future of PyMC3, or: Theano is Dead, Long Live Theano " for historical background & next steps

• And more!

PyMC3 Example

Using the coin-toss example, let's say we have flip a coin 4 times.

We have some intuition that the probability from coin tosses are 0.35 and we have some idea of what the probability of the coin tosses would be.

We would like to see the probability of flipping the coin with our prior assumption and 4 attempts.

trials = 4

𝜃_real = 0.35
data = stats.bernoulli.rvs(p=𝜃_real, size =trials)

Trials:  4
Bias (𝜃):  0.35
Data:  [1 0 0 0]

To go about this, we would be leveraging PyMC3 in the following procedure:

1. Instantiate (model) Object
2. Inference
   • Calculate output distribution
3. Predict & Evaluate

1. Instantiate Object

import pymc3 as pm

with pm.Model() as coin_flip_model:
    𝜃 = pm.Beta('𝜃', alpha=1, beta=1)
    y = pm.Bernoulli('y', p = 𝜃, observed=data)

2.Inference

We now perform inference to approximate the posterior distribution.

The options for this inference is either by:

• sampling
• variational inference

with coin_flip_model:
    trace = pm.sample(2_000, random_seed=123, return_inferencedata=True)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
/Users/CloudChaoszero/opt/anaconda3/envs/general/lib/python3.7/site-packages/aesara/graph/fg.py:501: UserWarning: Variable Elemwise{Switch}.0 cannot be replaced; it isn't in the FunctionGraph
  f"Variable {var} cannot be replaced; it isn't in the FunctionGraph"
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [𝜃]

 100.00% [6000/6000 00:04<00:00 Sampling 2 chains, 0 divergences]
Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 19 seconds.

or....

with pm.Model() as coin_flip_model:
    𝜃 = pm.Beta('𝜃', alpha=1, beta=1) # Prior
    y = pm.Bernoulli('y', p = 𝜃, observed=data) # Likelihood
    trace = pm.sample(2_000, random_seed=123, return_inferencedata=True) #Draw Samples

We just did the following:

• Instantiate a model called coin_flip_model
• Set Prior (theta) as a Beta Distribution
• Set the Likelihood (y) as a Bernoulli Distribution
• Lastly, use PyMC3's pm.sample function to create a trace with 2,000 drawn samples from the posterior distribution

Now, for PyMC3 versions 3.9+, we just saw the parameter return_inferencedata=True to return an arviz.InferenceData object, instead of a default MultiTrace object.

trace = pm.sample(draws = 20_000, return_inferencedata=True)

InferenceData has many advantages, compared to a MultiTrace: For example it can be saved/loaded from a file, and can also carry additional (meta)data such as date/version, or posterior predictive distributions.

trace

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:  (chain: 2, draw: 2000)
  Coordinates:
    * chain    (chain) int64 0 1
    * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 1993 1994 1995 1996 1997 1998 1999
  Data variables:
      𝜃        (chain, draw) float64 0.09391 0.1793 0.211 ... 0.2719 0.2625 0.1611
  Attributes:
      created_at:                 2021-04-22T07:56:16.445509
      arviz_version:              0.11.2
      inference_library:          pymc3
      inference_library_version:  3.11.1
      sampling_time:              11.031096935272217
      tuning_steps:               1000

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 2, draw: 2000, y_dim_0: 4)
  Coordinates:
    * chain    (chain) int64 0 1
    * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 1993 1994 1995 1996 1997 1998 1999
    * y_dim_0  (y_dim_0) int64 0 1 2 3
  Data variables:
      y        (chain, draw, y_dim_0) float64 -2.365 -0.09862 ... -0.1757 -0.1757
  Attributes:
      created_at:                 2021-04-22T07:56:16.683197
      arviz_version:              0.11.2
      inference_library:          pymc3
      inference_library_version:  3.11.1

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:             (chain: 2, draw: 2000)
  Coordinates:
    * chain               (chain) int64 0 1
    * draw                (draw) int64 0 1 2 3 4 5 ... 1995 1996 1997 1998 1999
  Data variables: (12/13)
      energy              (chain, draw) float64 5.204 4.924 5.009 ... 3.915 4.358
      perf_counter_diff   (chain, draw) float64 0.0006529 0.0002448 ... 0.0002258
      tree_depth          (chain, draw) int64 2 1 2 2 2 1 2 2 ... 1 2 2 1 1 2 1 1
      process_time_diff   (chain, draw) float64 0.000639 0.000246 ... 0.000227
      diverging           (chain, draw) bool False False False ... False False
      energy_error        (chain, draw) float64 0.4592 -0.2174 ... 0.007773 0.1648
      ...                  ...
      n_steps             (chain, draw) float64 3.0 1.0 3.0 3.0 ... 3.0 1.0 1.0
      step_size           (chain, draw) float64 1.5 1.5 1.5 ... 1.012 1.012 1.012
      max_energy_error    (chain, draw) float64 0.7756 -0.2174 ... 0.007773 0.1648
      step_size_bar       (chain, draw) float64 1.283 1.283 1.283 ... 1.297 1.297
      perf_counter_start  (chain, draw) float64 28.29 28.3 28.3 ... 30.78 30.78
      acceptance_rate     (chain, draw) float64 0.6318 1.0 0.8711 ... 0.9923 0.848
  Attributes:
      created_at:                 2021-04-22T07:56:16.450183
      arviz_version:              0.11.2
      inference_library:          pymc3
      inference_library_version:  3.11.1
      sampling_time:              11.031096935272217
      tuning_steps:               1000

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (y_dim_0: 4)
  Coordinates:
    * y_dim_0  (y_dim_0) int64 0 1 2 3
  Data variables:
      y        (y_dim_0) int32 1 0 0 0
  Attributes:
      created_at:                 2021-04-22T07:56:16.684283
      arviz_version:              0.11.2
      inference_library:          pymc3
      inference_library_version:  3.11.1

  
  

3. Predict & EvaluatePredict & Evaluate

Visualizing with ArviZ

Using the InferenceData object, we can leverage ArviZ to visualize our results.

ArviZ is a Python package for exploratory analysis of Bayesian models. Includes functions for posterior analysis, data storage, sample diagnostics, model checking, and comparison.

The goal is to provide backend-agnostic tools for diagnostics and visualizations of Bayesian inference in Python

We can plot the "trace" object using az.plot_trace

import arviz as az
az.plot_trace(trace)

array([[<AxesSubplot:title={'center':'𝜃'}>,
        <AxesSubplot:title={'center':'𝜃'}>]], dtype=object)

Or maybe we want a summary table for the trace, similar to Pandas pd.DataFrame().describe()

az.summary(trace)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
𝜃	0.33	0.181	0.016	0.642	0.005	0.003	1563.0	2202.0	1.0

Maybe we want to observe the aforementioned information in a different form, seen by arviZ's plot_posterior of the trace object

az.plot_posterior(trace)

<AxesSubplot:title={'center':'𝜃'}>

A summary

So, we went through the following:

• What is Bayesian Modeling?
• Brief Introduction to Probability
• Getting Started with PyMC3

However, there is so much more to learn about PyMC3! In the meantime, you can learn more about PyMC3 via their documentation...

or wait for another exciting BAyPIGgies Meetup! (S/O to Jeff Fischer & team)

With that being said

Thank you,

• RaulingAverage.dev
• Twitter: @RaulingAverage
• Notebook @ CloudChaoszero/Presentations/BayPiggies/Survival-Analysis

Resources

• Introduction to PyMC3 by Tung N. Nguyen
• PyMC3 Documentation
• Arviz Documentation
• For consulting services: There is a consultancy group called PyMC3 Labs
  • PyMC3 labs to help answer tough Business questions that require leveraging Bayesian Modeling.