import arviz as az
+import bambi as bmb
+import numpy as np
+import pandas as pd
+import pymc as pm
+import pytensor.tensor as pt
+
+from matplotlib import pyplot as plt
Horseshoe Prior
+In this example, we will use the Horseshoe Prior (Carvalho et al., 2009) to model a large number of variables, with only a few slopes being significantly different from zero.
+Here is what we did:
+-
+
- We defined an intercept. +
- We defined a vector of 50 betas, 5 of which were drawn from a normal(5,1) distribution, and then assigned a random sign. +
- We created the design matrix with normal(0,1) entries and set \(\sigma\) to 1. +
- We calculated the deterministic means \(\mu\) using the intercept and the design matrix multiplied by the betas. +
- We simulated 100 response variables (observations) from a normal distribution with mean \(\mu\) and standard deviation \(\sigma\). +
Next, we proceeded with the Bayesian estimation of the model. We proposed the horseshoe prior, for which the following parameters were calculated:
+\[\mu_i = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + ... + \beta_p x_{pi}\]
+\[y_i \sim N(\mu_i, \sigma^2)\]
+\[\alpha \sim N(0,1)\]
+\[\beta_j \sim N(0,\lambda_j^2 \tau^2)\]
+\[\lambda_j \sim C^+(0,1)\]
+\[\tau \sim T^+(df=3)\]
+\[\sigma^2 \sim N^+(0,1)\]
+= 50
+ D = 5
+ D0
+= 123456789 # for reproducibility
+ SEED
+= np.random.default_rng(SEED)
+ rng
+= rng.uniform(-3, 3) # simulate an intercept
+ INTERCEPT
+= np.zeros(D)
+ COEF # Simulate the slopes for significant variables
+= rng.choice([-1, 1], size=D0) * rng.normal(5, 1, size=D0)
+ COEF[:D0]
+= 100
+ N = rng.normal(size=(N, D))
+ X = 1.
+ SIGMA # Simulate the data
+= INTERCEPT + X.dot(COEF) + rng.normal(0, SIGMA, size=N) y
Here we create the dataframe and the term name for the set of variables, to define the formula.
+= pd.DataFrame(X)
+ df = [f"x{i}" for i in range(X.shape[1])]
+ df.columns "y"] = y df[
= "c(" + ", ".join([f"x{i}" for i in range(X.shape[1])]) + ")"
+ term_name = f"y ~ {term_name}"
+ formula formula
'y ~ c(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49)'
+Finally, we call the Horseshoe prior and create the model
+= {
+ priors "Horseshoe"),
+ term_name: bmb.Prior(
+ }= bmb.Model(formula, df, priors=priors)
+ model "predictors"})
+ model.set_alias({term_name:
+
+ model.build() model.graph()
= model.fit(target_accept = 0.95, chains=2) idata
Auto-assigning NUTS sampler...
+Initializing NUTS using jitter+adapt_diag...
+Multiprocess sampling (2 chains in 2 jobs)
+NUTS: [sigma, Intercept, predictors_tau, predictors_lam, predictors_raw]
+Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 444 seconds.
+There were 64 divergences after tuning. Increase `target_accept` or reparameterize.
+Chain 0 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
+Chain 1 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
+We recommend running at least 4 chains for robust computation of convergence diagnostics
+= {
+ priors "Horseshoe", tau_nu = 3, lam_nu = 3),
+ term_name: bmb.Prior(
+ }= bmb.Model(formula, df, priors=priors)
+ model "predictors"})
+ model.set_alias({term_name:
+
+ model.build() model.graph()
= model.fit(target_accept = 0.97, chains=2) idata
Auto-assigning NUTS sampler...
+Initializing NUTS using jitter+adapt_diag...
+Multiprocess sampling (2 chains in 2 jobs)
+NUTS: [sigma, Intercept, predictors_tau, predictors_lam, predictors_raw]
+Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 293 seconds.
+There were 29 divergences after tuning. Increase `target_accept` or reparameterize.
+We recommend running at least 4 chains for robust computation of convergence diagnostics
+= az.plot_forest(
+ ax,
+ idata, =["predictors"],
+ var_names={"predictors_dim": range(D0)},
+ coords='ridgeplot',
+ kind=False,
+ ridgeplot_truncate=0.5,
+ ridgeplot_alpha=0.95,
+ hdi_prob=True,
+ combined=(8, 6)
+ figsize
+ )-1], ax.get_yticks(), c='C1', label="Actual value");
+ ax.scatter(COEF[:D0][::r"$\beta_i$");
+ ax.set_xlabel(=None, top=1.55 * ax.get_yticks().max())
+ ax.set_ylim(bottomrange(D0)[::-1]);
+ ax.set_yticklabels(r"$i$");
+ ax.set_ylabel(='upper center');
+ ax.legend(loc"Posterior distribution of nonzero coefficients"); ax.set_title(
%load_ext watermark
+%watermark -n -u -v -iv -w
Last updated: Thu Aug 22 2024
+
+Python implementation: CPython
+Python version : 3.11.9
+IPython version : 8.24.0
+
+pandas : 2.2.2
+numpy : 1.26.4
+bambi : 0.14.1.dev12+g64e57423.d20240730
+arviz : 0.18.0
+matplotlib: 3.8.4
+pymc : 5.16.1
+pytensor : 2.23.0
+
+Watermark: 2.4.3
+
+