import arviz as az
import pymc as pm
from pymc.math import switch, ge
import pandas as pd
import numpy as np

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

This example introduces data containers, tracking of deterministic variables, 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)
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()
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, nuts_sampler="numpyro", target_accept=0.95)
    pm.sample_posterior_predictive(trace, extend_inferencedata=True)

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.853
      20.516
      -69.087
      12.323
      0.229
      0.163
      8059.0
      9626.0
      1.0
    

    

      beta[1]
      0.318
      4.615
      -8.586
      9.608
      0.053
      0.038
      7612.0
      8683.0
      1.0
    

    

      beta[2]
      3.913
      1.304
      1.373
      6.510
      0.013
      0.010
      9372.0
      10963.0
      1.0
    

    

      beta[3]
      19.693
      9.117
      1.982
      37.957
      0.095
      0.068
      9319.0
      10143.0
      1.0
    

    

      tau
      0.010
      0.003
      0.005
      0.015
      0.000
      0.000
      10086.0
      10429.0
      1.0
    

    

      sigma
      10.435
      1.492
      7.818
      13.410
      0.015
      0.011
      10086.0
      10429.0
      1.0
    

  











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









r2        0.576385
r2_std    0.075948
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.






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#


I think this is the first time we’ve seen the 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 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)






We should also probably use pm.math.ge for greater than or equal, as well, so:


pm.math.switch(ge(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
    plasma_aa_ns = pm.Normal(
        "nonsmokers_aa",
        mu=mu_nonsmokers,
        sigma=sigma_nonsmokers,
        observed=nonsmokers,
    )
    plasma_aa_s = pm.Normal(
        "smokers_aa", mu=mu_smokers, sigma=sigma_smokers, observed=smokers
    )

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

    trace = pm.sample(5000)










Hide code cell output




Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [tau_nonsmokers, mu_nonsmokers, tau_smokers, mu_smokers]
/Users/aaron/miniforge3/envs/pymc_macos15/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()





/Users/aaron/miniforge3/envs/pymc_macos15/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()







Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 3 seconds.















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
    

  
  
    

      mu_nonsmokers
      0.912
      0.045
      0.824
      1.002
      0.000
      0.000
      16807.0
      13465.0
      1.0
    

    

      mu_smokers
      0.976
      0.163
      0.659
      1.311
      0.002
      0.001
      12810.0
      9964.0
      1.0
    

    

      tau_nonsmokers
      22.597
      6.639
      10.794
      36.012
      0.050
      0.035
      17077.0
      13026.0
      1.0
    

    

      tau_smokers
      6.546
      3.471
      0.933
      13.314
      0.028
      0.020
      13354.0
      9401.0
      1.0
    

    

      test_mu
      0.666
      0.472
      0.000
      1.000
      0.004
      0.003
      16412.0
      16412.0
      1.0
    

    

      prec_ratio
      4.813
      4.152
      0.690
      11.976
      0.043
      0.031
      13374.0
      11204.0
      1.0
    

  











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









Last updated: Thu Nov 14 2024

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

pytensor: 2.26.0

numpy : 1.26.4
pymc  : 5.18.0
arviz : 0.20.0
pandas: 2.2.3