diff --git a/docs/source/learn/core_notebooks/vector_variables.ipynb b/docs/source/learn/core_notebooks/vector_variables.ipynb
new file mode 100644
index 0000000000..894b5dd4e0
--- /dev/null
+++ b/docs/source/learn/core_notebooks/vector_variables.ipynb
@@ -0,0 +1,995 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "777d3f58",
+ "metadata": {},
+ "source": [
+ "# Working with vector valued variables for multiple groups in PyMC\n",
+ "\n",
+ "This notebook is written as a response to a recurring question on GitHub and Discourse: \n",
+ "*“How do I handle multiple groups of data in PyMC and slice vector random variables correctly?”*\n",
+ "\n",
+ "The user’s example had several groups of observations and a vector of `mu` and `sigma` values, one for each group. \n",
+ "The confusing part was how to turn the group labels into something PyMC can index cleanly, and how to connect those indices to vector valued priors.\n",
+ "\n",
+ "To make the idea clear, I built a very small example with two levels:\n",
+ "\n",
+ "- a category (like Beverage or Snack)\n",
+ "- a family inside each category\n",
+ "\n",
+ "The goal here isn’t to create a big statistical model,\n",
+ "it’s just to show the exact pattern for:\n",
+ "- factorizing group labels \n",
+ "- building the mapping between levels \n",
+ "- creating vector RVs \n",
+ "- and slicing them correctly inside the likelihood \n",
+ "\n",
+ "Once you see this simple version working, the structure becomes easy to reuse in real models.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5269058d",
+ "metadata": {},
+ "source": [
+ "## 1. Setup\n",
+ "\n",
+ "Before jumping into the modeling part, we just import the usual tools: NumPy, pandas, ArviZ, and PyMC. \n",
+ "I also set a random seed so the results don’t wiggle around every time this notebook runs.\n",
+ "\n",
+ "One small thing I like doing , and it helps a lot later is defining named coordinates for the different group levels we’ll work with. \n",
+ "These labels make ArviZ’s output much easier to read because the plots will show actual names instead of axis numbers.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "e8f207b4",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "PyMC version: 5.26.1\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Importing the basic libraries we need\n",
+ "# Nothing special here, just the usual stack for PyMC work\n",
+ "import arviz as az\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "\n",
+ "import pymc as pm\n",
+ "\n",
+ "rng = np.random.default_rng(42)\n",
+ "az.style.use(\"arviz-darkgrid\")\n",
+ "\n",
+ "print(\"PyMC version:\", pm.__version__)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "113aa86d",
+ "metadata": {},
+ "source": [
+ "## 2. A small dataset to illustrate the idea\n",
+ "\n",
+ "The original GitHub issue used five groups of synthetic data. \n",
+ "To keep things intuitive here, I’m using a slightly more “real world sounding” example: categories and families.\n",
+ "\n",
+ "The dataset has three columns:\n",
+ "\n",
+ "- `category` \n",
+ "- `family` (which lives inside a category) \n",
+ "- `sales` (just some numeric values we’ll fit a model to)\n",
+ "\n",
+ "The exact numbers don’t matter what matters is that each row belongs to one family, and each family belongs to one category. \n",
+ "That structure is exactly the situation the user in the issue was dealing with.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "fdbe45c3",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "
category
\n",
+ "
family
\n",
+ "
sales
\n",
+ "
\n",
+ " \n",
+ " \n",
+ "
\n",
+ "
0
\n",
+ "
Beverage
\n",
+ "
Tea
\n",
+ "
15.609434
\n",
+ "
\n",
+ "
\n",
+ "
1
\n",
+ "
Beverage
\n",
+ "
Milk
\n",
+ "
7.920032
\n",
+ "
\n",
+ "
\n",
+ "
2
\n",
+ "
Beverage
\n",
+ "
Soft Drinks
\n",
+ "
21.500902
\n",
+ "
\n",
+ "
\n",
+ "
3
\n",
+ "
Snack
\n",
+ "
Chips
\n",
+ "
8.881129
\n",
+ "
\n",
+ "
\n",
+ "
4
\n",
+ "
Snack
\n",
+ "
Nuts
\n",
+ "
1.097930
\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ],
+ "text/plain": [
+ " category family sales\n",
+ "0 Beverage Tea 15.609434\n",
+ "1 Beverage Milk 7.920032\n",
+ "2 Beverage Soft Drinks 21.500902\n",
+ "3 Snack Chips 8.881129\n",
+ "4 Snack Nuts 1.097930"
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# Tiny example dataset\n",
+ "data = pd.DataFrame(\n",
+ " {\n",
+ " \"category\": [\"Beverage\", \"Beverage\", \"Beverage\", \"Snack\", \"Snack\"],\n",
+ " \"family\": [\"Tea\", \"Milk\", \"Soft Drinks\", \"Chips\", \"Nuts\"],\n",
+ " }\n",
+ ")\n",
+ "\n",
+ "# Pretend we observed some sales numbers (generated from a simple ground truth)\n",
+ "true_sales = {\n",
+ " \"Tea\": 15.0,\n",
+ " \"Milk\": 10.0,\n",
+ " \"Soft Drinks\": 20.0,\n",
+ " \"Chips\": 7.0,\n",
+ " \"Nuts\": 5.0,\n",
+ "}\n",
+ "data[\"sales\"] = [true_sales[f] + rng.normal(0, 2.0) for f in data[\"family\"]]\n",
+ "\n",
+ "data"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "af402ace",
+ "metadata": {},
+ "source": [
+ "## 3. Turning text labels into indices and making the mapping\n",
+ "\n",
+ "The thing that usually trips people up is how to connect the observed labels to vector-valued priors.\n",
+ "\n",
+ "PyMC can only index arrays with integers, so the first step is simply:\n",
+ "- convert the category names into integer codes\n",
+ "- convert the family names into integer codes\n",
+ "\n",
+ "Once we have those, we make a small array that maps each family to its category. \n",
+ "For example: if the 0th family belongs to the 1st category, the mapping array contains a 1 at position 0.\n",
+ "\n",
+ "This mapping array is the heart of the whole trick. \n",
+ "It tells PyMC how the lower level “inherits” from the upper level.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bd665a51",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Category labels: ['Beverage', 'Snack']\n",
+ "Family labels: ['Tea', 'Milk', 'Soft Drinks', 'Chips', 'Nuts']\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "
category
\n",
+ "
family
\n",
+ "
sales
\n",
+ "
cat_code
\n",
+ "
fam_code
\n",
+ "
\n",
+ " \n",
+ " \n",
+ "
\n",
+ "
0
\n",
+ "
Beverage
\n",
+ "
Tea
\n",
+ "
15.609434
\n",
+ "
0
\n",
+ "
0
\n",
+ "
\n",
+ "
\n",
+ "
1
\n",
+ "
Beverage
\n",
+ "
Milk
\n",
+ "
7.920032
\n",
+ "
0
\n",
+ "
1
\n",
+ "
\n",
+ "
\n",
+ "
2
\n",
+ "
Beverage
\n",
+ "
Soft Drinks
\n",
+ "
21.500902
\n",
+ "
0
\n",
+ "
2
\n",
+ "
\n",
+ "
\n",
+ "
3
\n",
+ "
Snack
\n",
+ "
Chips
\n",
+ "
8.881129
\n",
+ "
1
\n",
+ "
3
\n",
+ "
\n",
+ "
\n",
+ "
4
\n",
+ "
Snack
\n",
+ "
Nuts
\n",
+ "
1.097930
\n",
+ "
1
\n",
+ "
4
\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ],
+ "text/plain": [
+ " category family sales cat_code fam_code\n",
+ "0 Beverage Tea 15.609434 0 0\n",
+ "1 Beverage Milk 7.920032 0 1\n",
+ "2 Beverage Soft Drinks 21.500902 0 2\n",
+ "3 Snack Chips 8.881129 1 3\n",
+ "4 Snack Nuts 1.097930 1 4"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "family_to_category mapping (family index → category index): [0 0 0 1 1]\n",
+ "Length (should equal number of families): 5 vs 5\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Factorize categories and families into integer codes\n",
+ "cat_codes, cat_labels = pd.factorize(data[\"category\"])\n",
+ "fam_codes, fam_labels = pd.factorize(data[\"family\"])\n",
+ "\n",
+ "data[\"cat_code\"] = cat_codes\n",
+ "data[\"fam_code\"] = fam_codes\n",
+ "\n",
+ "print(\"Category labels:\", list(cat_labels))\n",
+ "print(\"Family labels:\", list(fam_labels))\n",
+ "display(data)\n",
+ "\n",
+ "# Build mapping\n",
+ "edges = data[[\"fam_code\", \"cat_code\"]].drop_duplicates().sort_values(\"fam_code\")\n",
+ "\n",
+ "family_to_category = edges[\"cat_code\"].to_numpy().astype(\"int64\")\n",
+ "\n",
+ "print(\"\\nfamily_to_category mapping (family index → category index):\", family_to_category)\n",
+ "print(\"Length (should equal number of families):\", len(family_to_category), \"vs\", len(fam_labels))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b6b39354",
+ "metadata": {},
+ "source": [
+ "## 4. Building the PyMC model\n",
+ "\n",
+ "Now that the indices and mapping are ready, the model becomes fairly straightforward.\n",
+ "\n",
+ "We set up:\n",
+ "- one global intercept\n",
+ "- one category effect per category\n",
+ "- one family effect per family\n",
+ "\n",
+ "The important part is that the family effect is centered on the category effect using the mapping array we created earlier. \n",
+ "That’s exactly the pattern the user in the GitHub issue needed but wasn’t sure how to set up.\n",
+ "\n",
+ "For each observation, the expected value is:\n",
+ "\n",
+ "global_mu \n",
+ "+ the effect for its category \n",
+ "+ the effect for its family\n",
+ "\n",
+ "Once the indexing is correct, PyMC takes it from there.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "c7b47307",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$$\n",
+ " \\begin{array}{rcl}\n",
+ " \\text{global\\_mu} &\\sim & \\operatorname{Normal}(10,~10)\\\\\\text{sigma\\_cat} &\\sim & \\operatorname{HalfNormal}(0,~5)\\\\\\text{category\\_effect} &\\sim & \\operatorname{Normal}(\\text{global\\_mu},~\\text{sigma\\_cat})\\\\\\text{sigma\\_fam} &\\sim & \\operatorname{HalfNormal}(0,~3)\\\\\\text{family\\_effect} &\\sim & \\operatorname{Normal}(f(\\text{category\\_effect}),~\\text{sigma\\_fam})\\\\\\text{sigma\\_obs} &\\sim & \\operatorname{HalfNormal}(0,~2)\\\\\\text{sales} &\\sim & \\operatorname{Normal}(f(\\text{family\\_effect}),~\\text{sigma\\_obs})\n",
+ " \\end{array}\n",
+ " $$"
+ ],
+ "text/plain": [
+ " global_mu ~ Normal(10, 10)\n",
+ " sigma_cat ~ HalfNormal(0, 5)\n",
+ "category_effect ~ Normal(global_mu, sigma_cat)\n",
+ " sigma_fam ~ HalfNormal(0, 3)\n",
+ " family_effect ~ Normal(f(category_effect), sigma_fam)\n",
+ " sigma_obs ~ HalfNormal(0, 2)\n",
+ " sales ~ Normal(f(family_effect), sigma_obs)"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "coords = {\n",
+ " \"category\": cat_labels,\n",
+ " \"family\": fam_labels,\n",
+ " \"obs\": np.arange(len(data)),\n",
+ "}\n",
+ "\n",
+ "with pm.Model(coords=coords) as model:\n",
+ " # Global mean\n",
+ " global_mu = pm.Normal(\"global_mu\", mu=10.0, sigma=10.0)\n",
+ "\n",
+ " # Level 0: category level vector variable\n",
+ " sigma_cat = pm.HalfNormal(\"sigma_cat\", sigma=5.0)\n",
+ " category_effect = pm.Normal(\n",
+ " \"category_effect\",\n",
+ " mu=global_mu,\n",
+ " sigma=sigma_cat,\n",
+ " dims=\"category\",\n",
+ " )\n",
+ "\n",
+ " # Level 1: family level vector variable, centered on its category\n",
+ " sigma_fam = pm.HalfNormal(\"sigma_fam\", sigma=3.0)\n",
+ " family_effect = pm.Normal(\n",
+ " \"family_effect\",\n",
+ " mu=category_effect[family_to_category],\n",
+ " sigma=sigma_fam,\n",
+ " dims=\"family\",\n",
+ " )\n",
+ "\n",
+ " # Observation model: each row uses its family's effect\n",
+ " sigma_obs = pm.HalfNormal(\"sigma_obs\", sigma=2.0)\n",
+ " mu = family_effect[fam_codes]\n",
+ "\n",
+ " sales = pm.Normal(\n",
+ " \"sales\",\n",
+ " mu=mu,\n",
+ " sigma=sigma_obs,\n",
+ " observed=data[\"sales\"].values,\n",
+ " dims=\"obs\",\n",
+ " )\n",
+ "\n",
+ "model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2a52eb49",
+ "metadata": {},
+ "source": [
+ "## 5. Sampling\n",
+ "\n",
+ "Here we let PyMC run `pm.sample()` to draw posterior samples. \n",
+ "Because the dataset is tiny, this finishes quickly.\n",
+ "\n",
+ "In real modeling work, I’d look at the diagnostics (R-hat, effective sample size, divergences). \n",
+ "But since this notebook is mainly about *how* to wire up the hierarchical structure, I’m keeping this part simple.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "b8044151",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Initializing NUTS using jitter+adapt_diag...\n",
+ "Multiprocess sampling (4 chains in 4 jobs)\n",
+ "NUTS: [global_mu, sigma_cat, category_effect, sigma_fam, family_effect, sigma_obs]\n",
+ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 47 seconds.\n",
+ "There were 132 divergences after tuning. Increase `target_accept` or reparameterize.\n",
+ "The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details\n",
+ "The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details\n"
+ ]
+ }
+ ],
+ "source": [
+ "%%capture sampling_output\n",
+ "\n",
+ "with model:\n",
+ " idata = pm.sample(\n",
+ " draws=1000,\n",
+ " tune=1000,\n",
+ " target_accept=0.95,\n",
+ " chains=4,\n",
+ " random_seed=42,\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "edfaacb3",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Sampling finished\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "