Skip to content

Commit 3d5a97f

Browse files
committed
Making all pm.sample return_inferencedata=False when needed
1 parent f2409c3 commit 3d5a97f

File tree

2 files changed

+869
-388
lines changed

2 files changed

+869
-388
lines changed

examples/case_studies/bayesian_ab_testing.ipynb

+853-368
Large diffs are not rendered by default.

myst_nbs/case_studies/bayesian_ab_testing.myst.md

+16-20
Original file line numberDiff line numberDiff line change
@@ -5,10 +5,6 @@ jupytext:
55
format_name: myst
66
format_version: 0.13
77
jupytext_version: 1.13.7
8-
kernelspec:
9-
display_name: Python 3 (ipykernel)
10-
language: python
11-
name: python3
128
---
139

1410
```{code-cell} ipython3
@@ -19,12 +15,12 @@ import arviz as az
1915
import matplotlib.pyplot as plt
2016
import numpy as np
2117
import pandas as pd
22-
import pymc3 as pm
23-
import pymc3.math as pmm
18+
import pymc as pm
19+
import pymc.math as pmm
2420
2521
from scipy.stats import bernoulli, expon
2622
27-
print(f"Running on PyMC3 v{pm.__version__}")
23+
print(f"Running on PyMC v{pm.__version__}")
2824
```
2925

3026
```{code-cell} ipython3
@@ -76,13 +72,13 @@ $$y_A \sim \mathrm{Binomial}(n = N_A, p = \theta_A), y_B \sim \mathrm{Binomial}(
7672

7773
With this, we can sample from the joint posterior of $\theta_A, \theta_B$.
7874

79-
You may have noticed that the Beta distribution is the conjugate prior for the Binomial, so we don't need MCMC sampling to estimate the posterior (the exact solution can be found in the VWO paper). We'll still demonstrate how sampling can be done with PyMC3 though, and doing this makes it easier to extend the model with different priors, dependency assumptions, etc.
75+
You may have noticed that the Beta distribution is the conjugate prior for the Binomial, so we don't need MCMC sampling to estimate the posterior (the exact solution can be found in the VWO paper). We'll still demonstrate how sampling can be done with PyMC though, and doing this makes it easier to extend the model with different priors, dependency assumptions, etc.
8076

8177
Finally, remember that our outcome of interest is whether B is better than A. A common measure in practice for whether B is better than is the _relative uplift in conversion rates_, i.e. the percentage difference of $\theta_B$ over $\theta_A$:
8278

8379
$$\mathrm{reluplift}_B = \theta_B / \theta_A - 1$$
8480

85-
We'll implement this model setup in PyMC3 below.
81+
We'll implement this model setup in PyMC below.
8682

8783
```{code-cell} ipython3
8884
@dataclass
@@ -127,7 +123,7 @@ We illustrate these points with prior predictive checks.
127123

128124
+++
129125

130-
Note that we can pass in arbitrary values for the observed data in these prior predictive checks. PyMC3 will not use that data when sampling from the prior predictive distribution.
126+
Note that we can pass in arbitrary values for the observed data in these prior predictive checks. PyMC will not use that data when sampling from the prior predictive distribution.
131127

132128
```{code-cell} ipython3
133129
weak_prior = ConversionModelTwoVariant(BetaPrior(alpha=100, beta=100))
@@ -139,12 +135,12 @@ strong_prior = ConversionModelTwoVariant(BetaPrior(alpha=10000, beta=10000))
139135

140136
```{code-cell} ipython3
141137
with weak_prior.create_model(data=[BinomialData(1, 1), BinomialData(1, 1)]):
142-
weak_prior_predictive = pm.sample_prior_predictive(samples=10000)
138+
weak_prior_predictive = pm.sample_prior_predictive(samples=10000, return_inferencedata=False)
143139
```
144140

145141
```{code-cell} ipython3
146142
with strong_prior.create_model(data=[BinomialData(1, 1), BinomialData(1, 1)]):
147-
strong_prior_predictive = pm.sample_prior_predictive(samples=10000)
143+
strong_prior_predictive = pm.sample_prior_predictive(samples=10000, return_inferencedata=False)
148144
```
149145

150146
```{code-cell} ipython3
@@ -204,9 +200,9 @@ def run_scenario_twovariant(
204200
generated = generate_binomial_data(variants, true_rates, samples_per_variant)
205201
data = [BinomialData(**generated[v].to_dict()) for v in variants]
206202
with ConversionModelTwoVariant(priors=weak_prior).create_model(data):
207-
trace_weak = pm.sample(draws=5000, return_inferencedata=True, cores=1, chains=2)
203+
trace_weak = pm.sample(draws=5000, cores=1, chains=2)
208204
with ConversionModelTwoVariant(priors=strong_prior).create_model(data):
209-
trace_strong = pm.sample(draws=5000, return_inferencedata=True, cores=1, chains=2)
205+
trace_strong = pm.sample(draws=5000, cores=1, chains=2)
210206
211207
true_rel_uplift = true_rates[1] / true_rates[0] - 1
212208
@@ -310,7 +306,7 @@ def run_scenario_bernoulli(
310306
generated = generate_binomial_data(variants, true_rates, samples_per_variant)
311307
data = [BinomialData(**generated[v].to_dict()) for v in variants]
312308
with ConversionModel(priors).create_model(data=data, comparison_method=comparison_method):
313-
trace = pm.sample(draws=5000, return_inferencedata=True, chains=2, cores=1)
309+
trace = pm.sample(draws=5000, chains=2, cores=1)
314310
315311
n_plots = len(variants)
316312
fig, axs = plt.subplots(nrows=n_plots, ncols=1, figsize=(3 * n_plots, 7), sharex=True)
@@ -492,7 +488,7 @@ data = [
492488

493489
```{code-cell} ipython3
494490
with RevenueModel(c_prior, mp_prior).create_model(data, "best_of_rest"):
495-
revenue_prior_predictive = pm.sample_prior_predictive(samples=10000)
491+
revenue_prior_predictive = pm.sample_prior_predictive(samples=10000, return_inferencedata=False)
496492
```
497493

498494
```{code-cell} ipython3
@@ -553,7 +549,7 @@ def run_scenario_value(
553549
with RevenueModel(conversion_rate_prior, mean_purchase_prior).create_model(
554550
data, comparison_method
555551
):
556-
trace = pm.sample(draws=5000, return_inferencedata=True, chains=2, cores=1)
552+
trace = pm.sample(draws=5000, chains=2, cores=1)
557553
558554
n_plots = len(variants)
559555
fig, axs = plt.subplots(nrows=n_plots, ncols=1, figsize=(3 * n_plots, 7), sharex=True)
@@ -655,9 +651,9 @@ There are many other considerations to implementing a Bayesian framework to anal
655651
* How do we plan the length and size of A/B tests using power analysis, if we're using Bayesian models to analyse the results?
656652
* Outside of the conversion rates (bernoulli random variables for each visitor), many value distributions in online software cannot be fit with nice densities like Normal, Gamma, etc. How do we model these?
657653

658-
Various textbooks and online resources dive into these areas in more detail. [Doing Bayesian Data Analysis](http://doingbayesiandataanalysis.blogspot.com/) by John Kruschke is a great resource, and has been translated to PyMC3 here: https://github.com/JWarmenhoven/DBDA-python.
654+
Various textbooks and online resources dive into these areas in more detail. [Doing Bayesian Data Analysis](http://doingbayesiandataanalysis.blogspot.com/) by John Kruschke is a great resource, and has been translated to PyMC here: https://github.com/JWarmenhoven/DBDA-python.
659655

660-
We also plan to create more PyMC3 tutorials on these topics, so stay tuned!
656+
We also plan to create more PyMC tutorials on these topics, so stay tuned!
661657

662658
---
663659

@@ -669,5 +665,5 @@ Author: [Cuong Duong](https://github.com/tcuongd) (2021-05-23)
669665

670666
```{code-cell} ipython3
671667
%load_ext watermark
672-
%watermark -n -u -v -iv -w -p theano,xarray
668+
%watermark -n -u -v -iv -w -p aesara,xarray
673669
```

0 commit comments

Comments
 (0)