Prediction*

import arviz as az
import numpy as np
import pandas as pd
import pymc as pm

rng = np.random.default_rng(6420)

8. Prediction*#

We previously covered Bayesian prediction in Unit 4, Lesson 13. Recall that the core of Bayesian prediction is the posterior predictive distribution

\[ f(x_{n+1}\vert x_{1} , \cdots , x_{n}) = \int_{\Theta}f(x_{n+1} \vert \mathbf{\theta})\pi(\theta \vert x_{1} , \cdots , x_{n})d\theta \]

Think of this as a “weighted average” of the likelihood for new data, with “weights” prescribed by for \(\pi(\theta \vert x_{1} , \cdots , x_{n})\). This is a departure from the “plug-in predictions” you may have encountered in a frequentist course. The result is a distribution of predictions. It should be clear by now that analytical methods are often intractable, and that is no different for the posterior predictive distribution. Fortunately, pymc has some tools that make this easier.

Taste of Cheese*#

Adapted from Unit 6: cheese.odc.

The link in the original .odc file is dead. I downloaded the data from here and have a copy here.

Problem Statement#

As cheddar cheese matures, a variety of chemical processes take place. The taste of matured cheese is related to the concentration of several chemicals in the final product. In a study of cheddar cheese from the Latrobe Valley of Victoria, Australia, samples of cheese were analyzed for their chemical composition and were subjected to taste tests. Overall taste scores were obtained by combining the scores from several tasters.

Can the score be predicted well by the predictors: Acetic, H2S, and Lactic?

data = pd.read_csv("../data/cheese.csv", index_col=0)
X = data[["Acetic", "H2S", "Lactic"]].to_numpy()
# add intercept column to X
X_aug = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1)
y = data["taste"].to_numpy()

# out of sample data for our later predictions
new_obs = np.array([[1.0, 5.0, 7.1, 1.5]])

data.head(5)
taste Acetic H2S Lactic
1 12.3 4.543 3.135 0.86
2 20.9 5.159 5.043 1.53
3 39.0 5.366 5.438 1.57
4 47.9 5.759 7.496 1.81
5 5.6 4.663 3.807 0.99

There are multiple ways to get predictions on out-of-sample data in PyMC. One way is to set up a shared variable using pm.Data in the original model, then using pm.set_data to change to the new observations before calling pm.sample_posterior_predictive. If you’re doing this, it’s useful to tie the shape of the likelihood to the number of rows in your X data. Note the shape=X_data.shape[0] line in the likelihood.

Another way is to do what we call “BUGS-style” posterior predictives. This is just setting up the posterior predictive within the model context, as labeled with a comment below.

with pm.Model() as m:
    # associate data with model
    X_data = pm.Data("X", X_aug)
    y_data = pm.Data("y", y)

    # priors
    beta = pm.Normal("beta", mu=0, sigma=1000, shape=X.shape[1] + 1)
    tau = pm.Gamma("tau", alpha=0.001, beta=0.001)
    sigma = pm.Deterministic("sigma", 1 / pm.math.sqrt(tau))

    mu = pm.math.dot(X_data, beta)

    # likelihood
    pm.Normal(
        "taste_score",
        mu=mu,
        sigma=sigma,
        observed=y_data,
        shape=X_data.shape[0],
    )

    oos_mean_response = pm.Deterministic(
        "oos_mean_response", pm.math.dot(new_obs, beta)
    )
    pm.Normal(
        "oos_prediction", oos_mean_response, sigma=sigma
    )  # BUGS-style post. pred.

    # start sampling
    trace = pm.sample(5000, target_accept=0.95, random_seed=rng)

Hide code cell output

 
NUTS[nutpie]: [beta, tau, oos_prediction]













az.summary(trace, ci_prob=0.95)










  
    

      
      mean
      sd
      eti95_lb
      eti95_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      beta[0]
      -29.4
      20.5
      -71
      11
      4347
      5445
      1.00
      0.31
      0.23
    

    

      beta[1]
      0.5
      4.6
      -8.6
      9.7
      4226
      5232
      1.00
      0.071
      0.054
    

    

      beta[2]
      3.9
      1.3
      1.3
      6.5
      7713
      8903
      1.00
      0.015
      0.012
    

    

      beta[3]
      19.6
      8.9
      1.7
      37
      7314
      8326
      1.00
      0.1
      0.079
    

    

      oos_prediction[0]
      29.9
      11.4
      7.4
      52
      15036
      12078
      1.00
      0.093
      0.069
    

    

      tau
      0.00978
      0.00271
      0.0052
      0.016
      14779
      12317
      1.00
      2.2e-05
      1.7e-05
    

    

      sigma
      10.41
      1.51
      8
      14
      14779
      12317
      1.00
      0.013
      0.011
    

    

      oos_mean_response[0]
      30
      3.76
      22
      37
      5632
      7331
      1.00
      0.05
      0.038
    

  









# single prediction on out-of-sample data (new_obs defined above)
pm.set_data({"X": new_obs}, model=m)
ppc = pm.sample_posterior_predictive(trace, model=m, predictions=True)











Hide code cell output





Sampling: [taste_score]

















az.summary(ppc.predictions, group="predictions", ci_prob=0.95)










  
    

      
      mean
      sd
      eti95_lb
      eti95_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      taste_score[0]
      29.9
      11.1
      8.2
      52
      15777
      17502
      1.00
      0.089
      0.066
    

  









%load_ext watermark
%watermark -n -u -v -iv -p pytensor









Last updated: Sun, 18 Oct 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

pytensor: 3.3.0

arviz : 1.3.0
numpy : 2.4.6
pandas: 3.0.5
pymc  : 6.3.1