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)

    # start sampling
    trace = pm.sample(5000, target_accept=0.95)
    pm.sample_posterior_predictive(trace, extend_inferencedata=True)
Hide code cell output 
Auto-assigning NUTS sampler...
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 22 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.574
      20.332
      -67.694
      11.618
      0.236
      0.167
      7415.0
      9686.0
      1.0
    

    

      beta[1]
      0.251
      4.608
      -8.735
      9.307
      0.055
      0.039
      6986.0
      8855.0
      1.0
    

    

      beta[2]
      3.936
      1.302
      1.341
      6.459
      0.014
      0.010
      8207.0
      8851.0
      1.0
    

    

      beta[3]
      19.645
      8.873
      1.844
      37.125
      0.091
      0.064
      9620.0
      10322.0
      1.0
    

    

      tau
      0.010
      0.003
      0.005
      0.015
      0.000
      0.000
      10541.0
      10135.0
      1.0
    

    

      sigma
      10.433
      1.514
      7.614
      13.377
      0.015
      0.011
      10541.0
      10135.0
      1.0
    

  











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









r2        0.577258
r2_std    0.075810
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()









mean            30.105300
sd              11.192933
hdi_2.5%         8.161233
hdi_97.5%       52.277533
mcse_mean        0.085700
mcse_sd          0.060667
ess_bulk     17107.966667
ess_tail     18409.700000
r_hat            1.000000
dtype: float64













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









Last updated: Mon Oct 28 2024

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

pytensor: 2.25.5

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