fix: LaTeX rendering, Exercise C worked example, Python 3.11+ requirements, Sphinx docs

.travis.yml

@@ -2,9 +2,8 @@
 
 language: python
 python:
-  - 3.8
-  - 3.7
-  - 3.6
+  - 3.11
+  - 3.12
 
 # Command to install dependencies, e.g. pip install -r requirements.txt --use-mirrors
 install: pip install -U tox-travis
@@ -25,4 +24,4 @@ deploy:
   on:
     tags: true
     repo: MacroMagic/bayesian_feature_selection
-    python: 3.8
+    python: 3.11

demo/horseshoe_mcmc_feature_selection_demo.ipynb

@@ -25,32 +25,32 @@
    "source": [
     "## 1) Mathematical intuition: why the horseshoe induces sparsity\n",
     "\n",
-    "For feature coefficients \\(\\beta_j\\), the (regularized) horseshoe hierarchy is:\n",
+    "For feature coefficients $\\beta_j$, the (regularized) horseshoe hierarchy is:\n",
     "\n",
-    "\\[\n",
+    "$$\n",
     "\\tau \\sim \\text{Half-Cauchy}(0, \\tau_0), \\quad\n",
     "\\lambda_j \\sim \\text{Half-Cauchy}(0, 1), \\quad\n",
     "c^2 \\sim \\text{Inverse-Gamma}(1,1)\n",
-    "\\]\n",
-    "\\[\n",
+    "$$\n",
+    "$$\n",
     "\\tilde{\\lambda}_j = \\sqrt{\\frac{c^2\\lambda_j^2}{c^2 + \\tau^2\\lambda_j^2}}, \\quad\n",
     "\\beta_j \\mid \\tau,\\lambda_j,c^2 \\sim \\mathcal{N}(0,\\tau^2\\tilde{\\lambda}_j^2)\n",
-    "\\]\n",
+    "$$\n",
     "\n",
     "For a Gaussian GLM:\n",
-    "\\[\n",
+    "$$\n",
     "\\eta_i = \\alpha + x_i^\\top\\beta, \\quad\n",
     "y_i \\sim \\mathcal{N}(\\eta_i, \\sigma^2)\n",
-    "\\]\n",
+    "$$\n",
     "\n",
     "### Key sparsity mechanism\n",
-    "- **Global shrinkage** \\(\\tau\\): pushes all coefficients toward zero.\n",
-    "- **Local shrinkage** \\(\\lambda_j\\): lets a few relevant features escape shrinkage.\n",
-    "- **Slab regularization** \\(c^2\\): avoids unrealistically huge coefficients.\n",
+    "- **Global shrinkage** $\\tau$: pushes all coefficients toward zero.\n",
+    "- **Local shrinkage** $\\lambda_j$: lets a few relevant features escape shrinkage.\n",
+    "- **Slab regularization** $c^2$: avoids unrealistically huge coefficients.\n",
     "\n",
-    "A useful view is the shrinkage factor \\(\\kappa_j = 1/(1 + \\tau^2\\lambda_j^2)\\):\n",
-    "- Noise features often have \\(\\kappa_j \\approx 1\\) (strong shrinkage to near zero).\n",
-    "- True signals often have smaller \\(\\kappa_j\\), preserving non-zero effects.\n"
+    "A useful view is the shrinkage factor $\\kappa_j = 1/(1 + \\tau^2\\lambda_j^2)$:\n",
+    "- Noise features often have $\\kappa_j \\approx 1$ (strong shrinkage to near zero).\n",
+    "- True signals often have smaller $\\kappa_j$, preserving non-zero effects.\n"
    ]
   },
   {
@@ -349,10 +349,10 @@
     "We focus on MCMC inference (`method='mcmc'`) for full posterior uncertainty.\n",
     "\n",
     "Rule-of-thumb initialization for global scale:\n",
-    "\\[\n",
+    "$$\n",
     "\\tau_0 \\approx \\frac{p_0}{p-p_0}\\cdot\\frac{1}{\\sqrt{n}}\n",
-    "\\]\n",
-    "where \\(p_0\\) is expected number of relevant features.\n"
+    "$$\n",
+    "where $p_0$ is expected number of relevant features.\n"
    ]
   },
   {
@@ -1533,6 +1533,152 @@
     "print(f\"  • The rule-of-thumb τ₀ ≈ {scale_global:.4f} balances the two\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "exercise_c_worked_md",
+   "metadata": {},
+   "source": [
+    "## 7c) Worked example — Exercise C: correlated predictors\n",
+    "\n",
+    "When predictors are highly correlated, the horseshoe posterior can spread mass\n",
+    "across multiple correlated features instead of concentrating it on a single signal.\n",
+    "This is a fundamental identifiability challenge—the data cannot distinguish which\n",
+    "correlated feature carries the true effect.\n",
+    "\n",
+    "**Setup**: 4 blocks of 5 correlated features each (within-block $\\rho = 0.9$),\n",
+    "plus 30 independent noise features.  One feature per block has a true nonzero\n",
+    "coefficient.\n",
+    "\n",
+    "**What to watch**:\n",
+    "- `method='beta'`: selects features with consistently large coefficients.  When\n",
+    "  multiple correlated features share posterior mass the *average* coefficient per\n",
+    "  feature is smaller, so inclusion probabilities can drop below the threshold.\n",
+    "- `method='lambda'`: uses the local shrinkage parameter $\\lambda_j$.  Correlated\n",
+    "  features tend to have elevated $\\lambda_j$ together, giving broader (sometimes\n",
+    "  over-inclusive) selection.\n",
+    "- `method='both'`: averages the two inclusion probabilities, balancing the\n",
+    "  specificity of beta-evidence with the sensitivity of lambda-evidence.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "exercise_c_worked_code",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ── Exercise C: correlated predictors ─────────────────────────────────────────\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "from bayesian_feature_selection import HorseshoeGLM, InferenceConfig\n",
+    "\n",
+    "np.random.seed(7)\n",
+    "\n",
+    "# ---- build correlated feature matrix ----------------------------------------\n",
+    "n          = 200\n",
+    "n_blocks   = 4      # number of correlated blocks\n",
+    "block_size = 5      # features per block\n",
+    "n_noise    = 30     # independent noise features\n",
+    "rho        = 0.9    # within-block correlation\n",
+    "\n",
+    "blocks = []\n",
+    "true_beta_c = np.zeros(n_blocks * block_size + n_noise)\n",
+    "for b in range(n_blocks):\n",
+    "    # All features in a block share a common factor + independent noise\n",
+    "    common = np.random.normal(size=(n, 1))\n",
+    "    indep  = np.random.normal(size=(n, block_size))\n",
+    "    block  = np.sqrt(rho) * common + np.sqrt(1 - rho) * indep\n",
+    "    blocks.append(block)\n",
+    "    # First feature in each block is the true signal\n",
+    "    true_beta_c[b * block_size] = 1.5\n",
+    "\n",
+    "X_noise = np.random.normal(size=(n, n_noise))\n",
+    "X_c = np.hstack(blocks + [X_noise])  # shape (200, 50)\n",
+    "p_c = X_c.shape[1]\n",
+    "\n",
+    "noise_sd_c = 1.0\n",
+    "y_c = X_c @ true_beta_c + np.random.normal(scale=noise_sd_c, size=n)\n",
+    "\n",
+    "# ---- fit horseshoe MCMC (short chain for demo speed) -----------------------\n",
+    "scale_global_c = (n_blocks / max(p_c - n_blocks, 1)) / np.sqrt(n)\n",
+    "scale_global_c = float(np.clip(scale_global_c, 0.05, 1.0))\n",
+    "\n",
+    "cfg_c = InferenceConfig(\n",
+    "    method=\"mcmc\",\n",
+    "    num_warmup=300,\n",
+    "    num_samples=500,\n",
+    "    num_chains=2,\n",
+    "    use_gpu=False,\n",
+    ")\n",
+    "model_c = HorseshoeGLM(family=\"gaussian\", scale_global=scale_global_c)\n",
+    "model_c.fit(X_c, y_c, config=cfg_c)\n",
+    "\n",
+    "# ---- compare three selection methods --------------------------------------\n",
+    "methods = ['beta', 'lambda', 'both']\n",
+    "results_c = {}\n",
+    "for m in methods:\n",
+    "    imp = model_c.get_feature_importance(method=m, threshold=0.5)\n",
+    "    selected = set(imp.loc[imp['selected'], 'feature_idx'].tolist())\n",
+    "    true_signals = set(range(0, n_blocks * block_size, block_size))  # indices 0,5,10,15\n",
+    "    tp = len(selected & true_signals)\n",
+    "    fp = len(selected - true_signals)\n",
+    "    fn = len(true_signals - selected)\n",
+    "    prec = tp / (tp + fp) if (tp + fp) > 0 else 0.0\n",
+    "    rec  = tp / (tp + fn) if (tp + fn) > 0 else 0.0\n",
+    "    f1   = 2 * prec * rec / (prec + rec) if (prec + rec) > 0 else 0.0\n",
+    "    results_c[m] = dict(selected=sorted(selected), TP=tp, FP=fp, FN=fn,\n",
+    "                        precision=prec, recall=rec, F1=f1)\n",
+    "\n",
+    "print(f\"True signal indices: {sorted(true_signals)}\")\n",
+    "print()\n",
+    "df_c = pd.DataFrame(results_c).T[['TP','FP','FN','precision','recall','F1']]\n",
+    "print(df_c.round(3).to_string())\n",
+    "\n",
+    "# ---- visualize inclusion probabilities by method --------------------------\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharey=True)\n",
+    "col_key = {'beta': 'beta_inclusion_prob', 'lambda': 'lambda_inclusion_prob',\n",
+    "           'both': 'combined_inclusion_prob'}\n",
+    "\n",
+    "for ax, m in zip(axes, methods):\n",
+    "    imp = model_c.get_feature_importance(method=m, threshold=0.5)\n",
+    "    prob_col = col_key[m]\n",
+    "    if prob_col not in imp.columns:\n",
+    "        # Fall back to whichever inclusion column exists\n",
+    "        prob_col = [c for c in imp.columns if 'inclusion' in c][0]\n",
+    "    probs = np.zeros(p_c)\n",
+    "    for _, row in imp.iterrows():\n",
+    "        probs[int(row['feature_idx'])] = row[prob_col]\n",
+    "\n",
+    "    colors = ['gold' if i in true_signals else 'steelblue' for i in range(p_c)]\n",
+    "    ax.bar(range(p_c), probs, color=colors, alpha=0.8, width=0.8)\n",
+    "    ax.axhline(0.5, color='red', linestyle='--', linewidth=1, label='threshold=0.5')\n",
+    "    ax.set_title(f\"method='{m}'\")\n",
+    "    ax.set_xlabel('Feature index')\n",
+    "    if ax == axes[0]:\n",
+    "        ax.set_ylabel('Inclusion probability')\n",
+    "    ax.legend(fontsize=8)\n",
+    "\n",
+    "fig.suptitle(\n",
+    "    'Exercise C — Correlated predictors: inclusion probabilities by method\\n'\n",
+    "    '(Gold bars = true signal features, dashed = selection threshold)',\n",
+    "    fontsize=11\n",
+    ")\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "# ---- discussion summary ---------------------------------------------------\n",
+    "print(\"\\nDiscussion:\")\n",
+    "print(\"  * beta  : requires the coefficient to be consistently large. \"\n",
+    "      \"Correlated features dilute the posterior so some signals may be missed.\")\n",
+    "print(\"  * lambda: less sensitive to exact coefficient magnitude; detects \"\n",
+    "      \"'unshrunk' features but can be less specific.\")\n",
+    "print(\"  * both  : averaging rewards features that show evidence in both \"\n",
+    "      \"criteria, often giving a balanced precision/recall tradeoff.\")\n",
+    "print(\"\\n  In high-correlation settings, model averaging or spike-and-slab \"\n",
+    "      \"priors that explicitly model group structure may outperform the horseshoe.\")\n"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "c8343096",

docs/conf.py

@@ -19,8 +19,31 @@
 #
 import os
 import sys
+from unittest.mock import MagicMock
+
 sys.path.insert(0, os.path.abspath('../src'))
 
+# Mock heavy C-extension dependencies so Sphinx can import the package without
+# requiring a full JAX / NumPyro installation in the docs build environment.
+_MOCK_MODULES = [
+    'jax', 'jax.numpy', 'jax.random', 'jax.scipy', 'jax.scipy.special',
+    'jaxlib',
+    'numpyro', 'numpyro.distributions', 'numpyro.infer', 'numpyro.primitives',
+    'numpyro.util',
+    'optax',
+    'arviz',
+    'sklearn', 'sklearn.preprocessing', 'sklearn.model_selection',
+    'pandas',
+    'numpy',
+    'matplotlib', 'matplotlib.pyplot', 'matplotlib.ticker',
+    'seaborn',
+    'yaml',
+    'typer',
+    'rich', 'rich.console', 'rich.progress',
+]
+for _mod_name in _MOCK_MODULES:
+    sys.modules.setdefault(_mod_name, MagicMock())
+
 import bayesian_feature_selection
 
 # -- General configuration ---------------------------------------------
@@ -31,7 +54,18 @@
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.viewcode', 'sphinx.ext.napoleon']
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.viewcode',
+    'sphinx.ext.napoleon',
+]
+
+# Mock imports for autodoc (complements the sys.modules mocking above)
+autodoc_mock_imports = [
+    'jax', 'jaxlib', 'numpyro', 'optax', 'arviz',
+    'sklearn', 'pandas', 'numpy', 'matplotlib', 'seaborn',
+    'yaml', 'typer', 'rich',
+]
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']

pyproject.toml

@@ -13,12 +13,11 @@ authors = [
 maintainers = [
   {name = "Hong Mu", email = "hm761@nyu.edu"}
 ]
+requires-python = ">=3.11"
 classifiers = [
     "Development Status :: 3 - Alpha",
     "Intended Audience :: Science/Research",
     "License :: OSI Approved :: MIT License",
-    "Programming Language :: Python :: 3.9",
-    "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
 ]

tox.ini

@@ -1,11 +1,10 @@
 [tox]
-envlist = py36, py37, py38, flake8
+envlist = py311, py312, flake8
 
 [travis]
 python =
-    3.8: py38
-    3.7: py37
-    3.6: py36
+    3.12: py312
+    3.11: py311
 
 [testenv:flake8]
 basepython = python