10. SSVS: Hald*

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

rng = np.random.default_rng(6420)

10. SSVS: Hald*#

Adapted from Haldssvs.odc.

SSVS can be looked at as adding an indicator variable to each covariate, then randomly turning each one on and off as you sample.

data = pd.read_csv("../data/hald_data.csv")
y = data["y"].to_numpy()
X = data.drop("y", axis=1).to_numpy()

X_centered = (X - X.mean(axis=0)) / X.std(axis=0)

_, p = X.shape
draws = 10000
chains = 4

with pm.Model() as m_svss:
    # SVSS prior
    delta = pm.Bernoulli("delta", p=0.5, shape=p)
    alpha = pm.Normal("alpha", 0, tau=0.1, shape=p)
    beta = pm.Deterministic("beta", delta * alpha)

    tau = pm.Gamma("tau", 0.1, 0.1)
    intercept = pm.Normal("intercept", 0, tau=0.001)

    mu = intercept + pm.math.dot(X_centered, beta)
    pm.Normal("likelihood", mu, tau=tau, observed=y)

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

Hide code cell output

 
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>BinaryGibbsMetropolis: [delta]
>NUTS: [alpha, tau, intercept]



Sampling 4 chains for 1_000 tune and 10_000 draw iterations (4_000 + 40_000 draws total) took 9 seconds.















az.summary(trace, var_names="delta", kind="stats")










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
    

  
  
    

      delta[0]
      1
      0.013
      1
      1
    

    

      delta[1]
      0.93
      0.25
      0
      1
    

    

      delta[2]
      0.42
      0.49
      0
      1
    

    

      delta[3]
      0.94
      0.23
      0
      1
    

  





Delta values, representing how often each feature was selected for the model, match the BUGS results. But what about looking at how often each model was selected? I really don’t want to code that quadruple for-loop (the model variable in the BUGS example) or try to vectorize it.


Let’s look in the trace. We can view the delta variable across all chains and samples.






trace.posterior.delta

























<xarray.DataArray 'delta' (chain: 4, draw: 10000, delta_dim_0: 4)> Size: 1MB
array([[[1, 1, 0, 1],
        [1, 1, 0, 1],
        [1, 1, 1, 1],
        ...,
        [1, 1, 1, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1]],

       [[1, 1, 0, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1],
        ...,
        [1, 1, 0, 1],
        [1, 1, 0, 1],
        [1, 1, 0, 1]],

       [[1, 1, 0, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1],
        ...,
        [1, 1, 0, 1],
        [1, 1, 0, 1],
        [1, 1, 0, 1]],

       [[1, 1, 0, 1],
        [1, 1, 0, 1],
        [1, 1, 0, 1],
        ...,
        [1, 1, 0, 1],
        [1, 1, 1, 1],
        [1, 1, 0, 1]]], shape=(4, 10000, 4))
Coordinates:
  * chain        (chain) int64 32B 0 1 2 3
  * draw         (draw) int64 80kB 0 1 2 3 4 5 ... 9994 9995 9996 9997 9998 9999
  * delta_dim_0  (delta_dim_0) int64 32B 0 1 2 3




We just need a way to count each unique combination of predictors. I didn’t see an obvious way to do that in xarray, so I’m copying the delta values to a NumPy array, then reshaping to include all draws for each chain along dimension 0, and finally getting the unique vectors of length p along with their counts, which we can use to calculate the probability of each model being selected.






rows = chains * draws

deltas = trace.posterior.delta.to_numpy()
models, counts = np.unique(
    deltas.reshape((rows, p)), axis=0, return_counts=True
)

for i in np.argsort(counts)[::-1]:
    print(f"{models[i]}: prob={counts[i]/rows:.3f}")









[1 1 0 1]: prob=0.525
[1 1 1 1]: prob=0.352
[1 0 1 1]: prob=0.056
[1 1 0 0]: prob=0.040
[1 1 1 0]: prob=0.015
[1 0 0 1]: prob=0.012
[0 1 1 1]: prob=0.000
[0 1 1 0]: prob=0.000
[0 0 1 1]: prob=0.000









The most chosen models are the same as the BUGS lecture version: intercept + x0 + x1 + x3 is selected over 50% of the time, intercept + x0 + x1 + x2 + x3 about a third of the time, and so on.






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









Last updated: Tue, 10 Nov 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.2