import arviz as az
import numpy as np
import pandas as pd
import pymc as pm
import matplotlib.pyplot as plt

8. Prediction*#

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()

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 
with pm.Model() as m:
    # associate data with model (this makes prediction easier)
    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])

    # start sampling
    trace = pm.sample(5000, target_accept=0.95)
    pm.sample_posterior_predictive(trace, extend_inferencedata=True)
Hide code cell output 
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta, tau]



Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 12 seconds.
Sampling: [taste_score]

















az.summary(trace, hdi_prob=0.95)













  
    

      
      mean
      sd
      hdi_2.5%
      hdi_97.5%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      beta[0]
      -28.978
      21.005
      -68.730
      14.372
      0.234
      0.166
      8043.0
      9468.0
      1.0
    

    

      beta[1]
      0.347
      4.758
      -9.408
      9.370
      0.055
      0.040
      7613.0
      9108.0
      1.0
    

    

      beta[2]
      3.915
      1.321
      1.359
      6.570
      0.014
      0.010
      9099.0
      10491.0
      1.0
    

    

      beta[3]
      19.644
      9.006
      1.539
      37.053
      0.089
      0.065
      10208.0
      10434.0
      1.0
    

    

      tau
      0.010
      0.003
      0.005
      0.015
      0.000
      0.000
      10815.0
      10486.0
      1.0
    

    

      sigma
      10.484
      1.540
      7.704
      13.522
      0.015
      0.011
      10815.0
      10486.0
      1.0
    

  











y_pred = trace.posterior_predictive.stack(sample=("chain", "draw"))[
    "taste_score"
].values.T
az.r2_score(y, y_pred)









r2        0.576145
r2_std    0.076283
dtype: float64









Results are pretty close to OpenBUGS:





mean


sd


MC_error


val2.5pc


median


val97.5pc


start


sample








beta0


-29.75


20.24


0.7889


-70.06


-29.75


11.11


1000


100001




beta1


0.4576


4.6


0.189


-8.716


0.4388


9.786


1000


100001




beta2


3.906


1.291


0.02725


1.345


3.912


6.47


1000


100001




beta3


19.79


8.893


0.2379


2.053


19.88


37.2


1000


100001




tau


0.009777


0.002706


2.29E-05


0.00522


0.009528


0.01575


1000


100001







PyMC gives some warnings about the model unless we increase the target_accept parameter of pm.sample, while BUGS doesn’t. This is because PyMC uses more diagnostics to check if there are any problems with its exploration of the parameter space. Divergences indicate bias in the results. BUGS will happily run this model without reporting any problems, but it doesn’t mean that there aren’t any.


For further reading, check out Diagnosing Biased Inference with Divergences.


It looks like there are multiple ways to get predictions on out-of-sample data in PyMC. The easiest 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.






# single prediction on out-of-sample data
new_obs = np.array([[1.0, 5.0, 7.1, 1.5]])
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]













The default behavior now is to give one prediction per y-value. The professor often asks for a single prediction based on the new data; the equivalent here would be to take the mean of the predicted values.






az.summary(ppc.predictions, hdi_prob=0.95)













  
    

      
      mean
      sd
      hdi_2.5%
      hdi_97.5%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      taste_score[0]
      29.961
      11.334
      7.922
      52.989
      0.084
      0.059
      18324.0
      19608.0
      1.0
    

  











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









Last updated: Sun Mar 09 2025

Python implementation: CPython
Python version       : 3.12.7
IPython version      : 8.29.0

pytensor: 2.26.4

pymc      : 5.19.1
numpy     : 1.26.4
arviz     : 0.20.0
matplotlib: 3.9.2
pandas    : 2.2.3