Create B-spline basis with 10 knots, degree 3 (cubic)

B = pmx.BSplineBasis(n_knots=10, degree=3) basis_matrix = B.build(x) # shape: (100, n_basis)

with pm.Model() as spline_model: # Coefficients for each basis function w = pm.Normal("w", mu=0, sigma=1, shape=basis_matrix.shape[1]) mu = pm.math.dot(basis_matrix, w) sigma = pm.HalfNormal("sigma", sigma=1) y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs)


### Penalized Splines with Smoothing Priors

Prevent overfitting by penalizing roughness via random walk priors on coefficients:

```python
with pm.Model() as penalized_spline:
    # Smoothing parameter controls wiggliness
    tau = pm.HalfCauchy("tau", beta=1)

    # Random walk prior on spline coefficients (penalizes second differences)
    w = pm.GaussianRandomWalk("w", sigma=tau, shape=basis_matrix.shape[1])

    mu = pm.math.dot(basis_matrix, w)
    sigma = pm.HalfNormal("sigma", sigma=1)
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs)

ZeroSumNormalBasis

For identifiable models where basis coefficients must sum to zero:

# Ensures sum-to-zero constraint on spline coefficients
B_zs = pmx.ZeroSumNormalBasis(n_knots=10, degree=3)

See references/splines.md for detailed spline reference.

Distributional Regression (GAMLSS-style)

Model ALL distribution parameters as functions of covariates, not just the mean.

Basic Pattern

with pm.Model() as dist_reg:
    # Model both mean AND variance as functions of covariates
    # Location model
    beta_mu = pm.Normal("beta_mu", 0, 1, shape=X.shape[1])
    mu = pm.math.dot(X, beta_mu)

    # Scale model (log link to ensure positivity)
    beta_sigma = pm.Normal("beta_sigma", 0, 0.5, shape=Z.shape[1])
    log_sigma = pm.math.dot(Z, beta_sigma)
    sigma = pm.math.exp(log_sigma)

    y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs)

Common Distributional Patterns

# Beta regression: both mu and phi vary with covariates
with pm.Model() as beta_reg:
    # Mean model (logit link)
    beta_mu = pm.Normal("beta_mu", 0, 1, shape=X.shape[1])
    mu = pm.math.sigmoid(pm.math.dot(X, beta_mu))

    # Precision model (log link)
    beta_phi = pm.Normal("beta_phi", 0, 0.5, shape=Z.shape[1])
    phi = pm.math.exp(pm.math.dot(Z, beta_phi))

    # Reparameterize: alpha = mu * phi, beta = (1-mu) * phi
    y = pm.Beta("y", alpha=mu * phi, beta=(1 - mu) * phi, observed=y_obs)

# Negative binomial: mu and alpha both depend on covariates
with pm.Model() as nb_dist:
    beta_mu = pm.Normal("beta_mu", 0, 1, shape=X.shape[1])
    mu = pm.math.exp(pm.math.dot(X, beta_mu))

    beta_alpha = pm.Normal("beta_alpha", 0, 0.5, shape=Z.shape[1])
    alpha = pm.math.exp(pm.math.dot(Z, beta_alpha))

    y = pm.NegativeBinomial("y", mu=mu, alpha=alpha, observed=y_obs)

See references/distributional.md for more patterns.

R2D2M2CP Prior

A prior for regression coefficients that reasons about the proportion of variance explained (R2) and decomposes it across predictors.

with pm.Model() as r2d2_model:
    # R2D2M2CP: specify expected R2 and concentration
    beta, r2 = pmx.R2D2M2CP(
        "beta",
        X,                      # Design matrix
        y_obs,                  # Observed response
        r2=0.5,                 # Prior mean for R2
        variance_concentration=None,  # Equal concentration (default)
    )
    mu = pm.math.dot(X, beta)
    sigma = pm.HalfNormal("sigma", sigma=1)
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs)

When to Use R2D2M2CP vs Horseshoe

Aspect	R2D2M2CP	Horseshoe
Interpretability	Reasons about R2 directly	Reasons about shrinkage
Sparsity	Soft shrinkage	Strong sparsity (near-zero coefficients)
Best for	Dense signals, moderate p	Sparse signals, large p
Tuning	Prior on R2 + concentration	Global/local shrinkage scales
Implementation	`pmx.R2D2M2CP()`	Manual or pmx helper

For general prior specification guidance and elicitation workflows, see the prior-elicitation skill.

Regularized Horseshoe Prior

For sparse regression where most coefficients are expected to be near zero:

with pm.Model() as horseshoe_model:
    n, p = X.shape

    # Expected number of relevant predictors
    p0 = 5
    # Global shrinkage — controls overall sparsity
    tau0 = p0 / (p - p0) / np.sqrt(n)
    tau = pm.HalfCauchy("tau", beta=tau0)

    # Local shrinkage — per-coefficient
    lam = pm.HalfCauchy("lam", beta=1, shape=p)

    # Regularized: slab component prevents unbounded coefficients
    c2 = pm.InverseGamma("c2", alpha=2, beta=1)
    lam_tilde = lam * pm.math.sqrt(c2 / (c2 + tau**2 * lam**2))

    beta = pm.Normal("beta", mu=0, sigma=tau * lam_tilde, shape=p)
    mu = pm.math.dot(X, beta)
    sigma = pm.HalfNormal("sigma", sigma=1)
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs)

See references/r2d2_horseshoe.md for detailed comparison.

Marginalization of Discrete Parameters

pmx.marginalize() analytically integrates out discrete latent variables, avoiding the need for specialized samplers.

with pm.Model() as mixture:
    # Mixture weights
    w = pm.Dirichlet("w", a=np.ones(3))

    # Component means
    mu = pm.Normal("mu", mu=[-5, 0, 5], sigma=1, shape=3)
    sigma = pm.HalfNormal("sigma", sigma=1, shape=3)

    # Discrete component assignment (will be marginalized)
    comp = pm.Categorical("comp", p=w, shape=len(y_obs))

    # Likelihood
    y = pm.Normal("y", mu=mu[comp], sigma=sigma[comp], observed=y_obs)

    # Marginalize out the discrete variable
    pmx.marginalize(model=mixture, rvs_to_marginalize=["comp"])

    # Now sample with standard NUTS — no need for CategoricalGibbsMetropolis
    idata = pm.sample()

Supported Distributions for Marginalization

Categorical assignments in mixture models
Bernoulli indicators in spike-and-slab
Discrete latent states in hidden Markov models

Benefits

Enables NUTS sampling for models with discrete parameters
Often dramatically improves sampling efficiency
Removes label switching issues in mixture models

Laplace Approximation

Fast approximate inference via second-order Taylor expansion at the MAP:

with pm.Model() as model:
    beta = pm.Normal("beta", 0, 1, shape=5)
    mu = pm.math.dot(X, beta)
    sigma = pm.HalfNormal("sigma", sigma=1)
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs)

    # Laplace approximation
    idata_laplace = pmx.fit_laplace(model)

When to Use Laplace

Scenario	Use Laplace?
Quick approximate inference	Yes
Large datasets, simple model	Yes — fast and accurate
Posterior is near-Gaussian	Yes — Laplace is exact for Gaussian posteriors
Multimodal posterior	No — Laplace finds only one mode
Heavy-tailed posteriors	No — underestimates tail uncertainty
Hierarchical models with few groups	No — posterior is often non-Gaussian
Model comparison (for screening)	Yes — fast screening before full MCMC

Comparing Laplace to MCMC

# Quick Laplace fit for model development
idata_laplace = pmx.fit_laplace(model)

# Full MCMC for final inference
idata_mcmc = pm.sample()

# Compare: Laplace posteriors should be similar to MCMC
# if the approximation is good
import arviz as az
az.plot_forest(
    [idata_laplace, idata_mcmc],
    model_names=["Laplace", "MCMC"],
)

Laplace is excellent for rapid iteration during model development. Use full MCMC for final inference and reporting.

pymc-extras

PyMC-Extras (pmx) — Advanced Extensions

Splines

BSplineBasis