Pyro is a deep probabilistic programming framework based on PyTorch. TensorFlow has its own PPL branch with an Edward taste and there is the inevitable PyMC3 as well but Pyro feels very natural and the API more direct than the aforementioned.
The example below is simplistic but shows what I mean by ‘natural’.
|
1 2 3 4 5 6 7 8 9 10 |
import torch import pyro import pyro.distributions as dist import numpy as np import matplotlib.pyplot as plt import seaborn as sea sea.set() sea.set %matplotlib inline |
You get immediate results with Pyro or should one say ‘eager execution’:
|
1 2 3 4 |
normal = dist.Normal(0, 1) x = normal.sample() print(f"sample: {x}") print(f"log probability: {normal.log_prob(x)}") |
Defining a variable is different than PyMC3. A normally distribute variable ‘size’ would be:
|
1 2 |
x = pyro.sample("size", normal) print(x) |
By combining probabilistic variables you can create a simplistic weather model liks so:
|
1 2 3 4 5 6 7 |
def weather(): cloudy = pyro.sample('cloudy', dist.Bernoulli(0.3)) cloudy = 'cloudy' if cloudy.item() == 1.0 else 'sunny' mean_temp = {'cloudy': 55.0, 'sunny': 75.0}[cloudy] scale_temp = {'cloudy': 10.0, 'sunny': 15.0}[cloudy] temp = pyro.sample('temp', dist.Normal(mean_temp, scale_temp)) return cloudy, temp.item() |
The logic is straightforward and you can get a few states:
|
1 2 |
for _ in range(3): print(weather()) |
Let’s model the sales of ice-cream in function of the weather:
|
1 2 3 4 5 |
def ice_cream_sales(): cloudy, temp = weather() expected_sales = 200. if cloudy == 'sunny' and temp > 80.0 else 50. ice_cream = pyro.sample('ice_cream', dist.Normal(expected_sales, 10.0)) return ice_cream |
Unlike other PPL frameworks you can get an idea directly from this definition
|
1 |
print(f"Possible sales of ice-cream: {ice_cream_sales():.2f}") |
The key-idea behind PPL is to have a marginal distribution of the relevant variable. Marginal means, in simple terms, getting a probability distribution with all de inner-dependencies integrated out.
To do this in Pyro you first define the posterior:
|
1 |
posterior = pyro.infer.Importance(ice_cream_sales, num_samples=100) |
and feed it to the EmpiricalMarginal:
|
1 |
marginal = pyro.infer.EmpiricalMarginal(posterior.run()) |
To see it in a histogran use something like this:
|
1 2 3 4 |
plt.hist([marginal().item() for _ in range(1000)], range=(5.0, 100.0)) plt.title("P(sales)") plt.xlabel("sales") plt.ylabel("#"); |
Considering the weather model and how the sales depends on it, the most likely value is around 50.
Both the way the story is assembled and how concrete the result is strikes me as fantastically clear. Bayesian reasoning and PPL can be abstract at first but this example brings it all together for me.
If we look at a linear regression in PyMC3 versus Pyro you get something like the following.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 |
%matplotlib inline import pymc3 as pm import numpy as np import pandas as pd import scipy.stats as stats import matplotlib.pyplot as plt import seaborn as sns palette = 'muted' sns.set_palette(palette) sns.set_color_codes(palette) np.set_printoptions(precision=2) pd.set_option('display.precision', 2) np.random.seed(1) N = 100 alfa_real = 2.5 beta_real = 0.9 eps_real = np.random.normal(0, 0.5, size=N) x = np.random.normal(10, 1, N) y_real = alfa_real + beta_real * x y = y_real + eps_real plt.figure(figsize=(10, 5)) plt.subplot(1, 2, 1) plt.plot(x, y, 'b.') plt.xlabel('$x$', fontsize=16) plt.ylabel('$y$', fontsize=16, rotation=0) plt.plot(x, y_real, 'k') plt.subplot(1, 2, 2) sns.kdeplot(y) plt.xlabel('$y$', fontsize=16) plt.tight_layout() with pm.Model() as model: alpha = pm.Normal('alpha', mu=0, sd=10) beta = pm.Normal('beta', mu=0, sd=1) epsilon = pm.HalfCauchy('epsilon', 5) mu = pm.Deterministic('mu', alpha + beta * x) y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y) start = pm.find_MAP() step = pm.Metropolis() trace = pm.sample(10000, step, start) pm.traceplot(trace) |
The same in Pyro look like this:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
import torch import pyro from pyro.distributions import Normal, HalfCauchy from pyro.infer.mcmc import NUTS, MCMC from pyro.infer import EmpiricalMarginal import seaborn as sns import matplotlib.pyplot as plt def model(X, y): b_est = pyro.sample("beta", Normal(0,1)) a_est = pyro.sample("alpha", Normal(0,1)) s2_est = pyro.sample("variance", HalfCauchy(1)) y_hat = pyro.sample("yhat", Normal(loc=b_est * X + a_est , scale=s2_est), obs=y) return y_hat X = torch.tensor(range(100)).float() b =torch.tensor(3.0).float() a =torch.tensor(2.0).float() noise = torch.randn((100,)).float() y = X * b + a + noise nuts_kernel = NUTS(model, adapt_step_size=True) mcmc_run = MCMC(nuts_kernel, num_samples=1000, warmup_steps=400).run(X, y) posterior_beta = EmpiricalMarginal(mcmc_run, "beta") posterior_alpha = EmpiricalMarginal(mcmc_run, "alpha") print(posterior_beta.mean) print(posterior_alpha.mean) fig, ax = plt.subplots(1,2,figsize=(15,5)) sns.distplot(posterior_beta.get_samples_and_weights()[0], ax=ax[0]) sns.distplot(posterior_alpha.get_samples_and_weights()[0], ax=ax[1]) fig.show() |
