Loading Data, Step Function, and Deterministic Variables*

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

rng = np.random.default_rng(6420)

5. Loading Data, Step Function, and Deterministic Variables*#

This example introduces the tracking of deterministic variables, Bayesian \(R^2\), and shows how to recreate the BUGS step function in PyMC.

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.

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)

data.head()
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

The main goal of this exercise is to practice loading the data and implementing a simple model. Let’s first specify the model structure, then write the code. Since our goal is to model the relationship between chemical compounds and perceived taste, we can use a Bayesian linear regression model:

\[\begin{split} \begin{align*} \texttt{taste}_{i} &\sim \mathcal{N}(\mu_{i} , \sigma^{2})\\ \mu_{i} &= \beta_{0} + \beta_{1}\cdot\texttt{Acetic}_{i} + \beta_{2}\cdot\texttt{H2S}_{i} + \beta_{3}\cdot\texttt{Lactic}_{i}\\ \boldsymbol{\beta} &\sim \mathcal{N}(0,\sigma^{2} = 1000^2)\\ \sigma &= \frac{1}{\sqrt{\tau}}\\ \tau &\sim \text{Gamma(0.001, 0.001)}\\ \end{align*} \end{split}\]

This model assumes that the taste score for each observation follows a normal distribution with mean \(\mu_i\) and variance \(\sigma^2\), where \(\mu_i\) is a linear combination of the chemical predictors. We place vague normal priors on the regression coefficients and a weakly informative Gamma prior on the precision parameter \(\tau\) (inverse variance).

It’s often more practical to express the linear model in matrix/vector form rather than writing out each term separately:

\[ \boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta} \]

Here \(\boldsymbol{\eta}\) is the linear predictor, without adding any random noise. If you’re not familiar with the matrix notation for multiple linear regression, \(\mathbf{X}\) and \(\boldsymbol{\beta}\) are:

\[\begin{split} \mathbf{X} = \begin{bmatrix} 1 & \texttt{Acetic}_1 & \texttt{H2S}_1 & \texttt{Lactic}_1 \\ 1 & \texttt{Acetic}_2 & \texttt{H2S}_2 & \texttt{Lactic}_2 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & \texttt{Acetic}_n & \texttt{H2S}_n & \texttt{Lactic}_n \end{bmatrix} \text{ and } \beta = \begin{bmatrix} \text{Coefficient for Intercept}\\ \text{Coefficient for }\texttt{Acetic}\\ \text{Coefficient for }\texttt{H2S}\\ \text{Coefficient for }\texttt{Lactic}\\ \end{bmatrix} \end{split}\]

\(\mathbf{X}\) is the design matrix: which contains all predictor columns along with a column of ones to account for the intercept \(\beta_{0}\). To make this clear in our code, we’ll call it X_aug.

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

with pm.Model() as m:
    # priors
    beta = pm.Normal("beta", mu=0, sigma=1000, shape=X.shape[1] + 1)
    tau = pm.Gamma("tau", 0.001, 0.001)  # precision
    pm.Deterministic("sigma", 1 / pm.math.sqrt(tau))  # sd
    pm.Deterministic("sigma2", 1 / tau)  # variance

    mu = pm.Deterministic("mu", pm.math.dot(X_aug, beta))

    # likelihood
    pm.Normal("taste_score", mu=mu, tau=tau, observed=y)

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

NUTS[nutpie]: [beta, tau]







The var_names argument below can help control the amount of arviz.summary output, depending on which variables are the most important to your analysis. Especially in this case; because mu is going to have a shape of 30, we might want to exclude it from the output.






az.summary(trace, var_names=["beta", "sigma"], ci_prob=0.95, ci_kind="hdi")










  
    

      
      mean
      sd
      hdi95_lb
      hdi95_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      beta[0]
      -28.9
      20.8
      -69
      12
      4638
      6043
      1.00
      0.31
      0.23
    

    

      beta[1]
      0.3
      4.7
      -9
      9.7
      4531
      5625
      1.00
      0.07
      0.053
    

    

      beta[2]
      3.92
      1.31
      1.3
      6.5
      8561
      9567
      1.00
      0.014
      0.011
    

    

      beta[3]
      19.6
      8.8
      2.3
      37
      8392
      9148
      1.00
      0.096
      0.072
    

    

      sigma
      10.43
      1.49
      8
      14
      14160
      11484
      1.00
      0.013
      0.011
    

  





This is how we find the Bayesian \(R^2\) [Gelman et al., 2019]. Note this is a different version from the lectures. This one is both easier to use in PyMC and generally a better metric compared to the lecture version. The lecture version can sometimes result in negative \(R^2\) values.


We saved mu above as a deterministic variable for this purpose. By default, the scale will be the standard deviation, but you can manually define it and switch to variance with the scale_kind argument.






az.bayesian_r2(
    trace,
    pred_mean="mu",
    scale="sigma",
    scale_kind="sd",
    ci_prob=0.95,
    ci_kind="hdi",
)









bayesian_R2(mean=0.62, hdi_lb=0.44, hdi_ub=0.78)













az.bayesian_r2(
    trace,
    pred_mean="mu",
    scale="sigma2",
    scale_kind="var",
    ci_prob=0.95,
    ci_kind="hdi",
)









bayesian_R2(mean=0.62, hdi_lb=0.44, hdi_ub=0.78)









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 may indicate bias in the results. Increasing target_accept can help if there aren’t any major underlying problems with the parameterization.


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.






Stress, Diet, and Acids#


Adapted from Unit 6: stressacids.odc.


In the study Interrelationships Between Stress, Dietary Intake, and Plasma Ascorbic Acid During Pregnancy conducted at the Virginia Polytechnic Institute and State University, the plasma ascorbic acid levels of pregnant women were compared for smokers versus non-smokers. Thirty-two women in the last three months of pregnancy, free of major health disorders, and ranging in age from 15 to 32 years were selected for the study. Prior to the collection of 20 ml of blood, the participants were told to avoid breakfast, forego their vitamin supplements, and avoid foods high in ascorbic acid content. From the blood samples, the plasma ascorbic acid values of each subject were determined in milligrams per 100 milliliters.




The purpose of this example in lectures was mostly just to show different ways to load data in BUGS. I’m not going to go into that too much, since there are a million ways to prepare your data in Python. In the next cell, I start with the data pasted from stressacids.odc, then use list comprehensions to create one list for smokers and one for nonsmokers.






# fmt: off
plasma = [0.97, 0.72, 1.00, 0.81, 0.62, 1.32, 1.24, 0.99, 0.90, 0.74,
          0.88, 0.94, 1.06, 0.86, 0.85, 0.58, 0.57, 0.64, 0.98, 1.09,
          0.92, 0.78, 1.24, 1.18, 0.48, 0.71, 0.98, 0.68, 1.18, 1.36,
          0.78, 1.64]

smo = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
       1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2]
# fmt: on

nonsmokers = [x for x, y in zip(plasma, smo) if y == 1]
smokers = [x for x, y in zip(plasma, smo) if y == 2]











BUGS step function#


BUGS defines the step function like this:


step(e) is 1 if e >= 0; 0 otherwise.


Keep in mind that in PyMC, step functions are sometimes how they refer to the algorithms used for sampling, as in NUTS or Metropolis. Just different terminology.


We can recreate the BUGS step function with pm.math.switch():


pm.math.switch(e >= 0, 1, 0)










How do I track non-random variables in PyMC?#


One nice thing about BUGS is you can easily track both deterministic and non-deterministic variables while sampling. For PyMC, you can wrap these in pm.Deterministic(). Just make sure to use pm.math functions where possible.






with pm.Model() as m:
    # priors
    tau_nonsmokers = pm.Gamma("tau_nonsmokers", alpha=0.0001, beta=0.0001)
    sigma_nonsmokers = 1 / pm.math.sqrt(tau_nonsmokers)
    mu_nonsmokers = pm.Normal("mu_nonsmokers", mu=0, sigma=100)

    tau_smokers = pm.Gamma("tau_smokers", alpha=0.0001, beta=0.0001)
    sigma_smokers = 1 / pm.math.sqrt(tau_smokers)
    mu_smokers = pm.Normal("mu_smokers", mu=0, sigma=100)

    # likelihood
    pm.Normal(
        "nonsmokers_aa",
        mu=mu_nonsmokers,
        sigma=sigma_nonsmokers,
        observed=nonsmokers,
    )
    pm.Normal(
        "smokers_aa", mu=mu_smokers, sigma=sigma_smokers, observed=smokers
    )

    testmu = pm.Deterministic(
        "test_mu", pm.math.switch(mu_smokers >= mu_nonsmokers, 1, 0)
    )
    r = pm.Deterministic("prec_ratio", tau_nonsmokers / tau_smokers)

    trace = pm.sample(5000, random_seed=rng)











Hide code cell output





NUTS[nutpie]: [tau_nonsmokers, mu_nonsmokers, tau_smokers, mu_smokers]













Notice the default output of arviz.summary is now an 89% equi-tailed credible interval.






az.summary(trace)










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      mu_nonsmokers
      0.9109
      0.0449
      0.84
      0.98
      18996
      14207
      1.00
      0.00033
      0.00025
    

    

      mu_smokers
      0.977
      0.161
      0.73
      1.2
      13218
      11546
      1.00
      0.0015
      0.0013
    

    

      tau_nonsmokers
      22.6
      6.74
      13
      34
      17277
      12408
      1.00
      0.05
      0.04
    

    

      tau_smokers
      6.55
      3.49
      2.1
      13
      13284
      11239
      1.00
      0.029
      0.028
    

    

      test_mu
      0.672
      0.47
      0
      1
      16015
      16015
      1.00
      0.0037
      0.0014
    

    

      prec_ratio
      4.78
      4.1
      1.4
      11
      14091
      11178
      1.00
      0.042
      0.14
    

  









%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