docs: add exploratory statistical significance tutorial

docs/Tutorial_Statistical_Significance.ipynb

@@ -0,0 +1,153 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Exploratory statistical significance for a matrix profile\n",
+    "\n",
+    "A matrix profile reports nearest-neighbor distances, but a small distance is not automatically evidence that a motif is meaningful. This tutorial shows one reproducible, exploratory way to compare an observed minimum profile distance with a null distribution made from circular block permutations. It is a teaching aid, not a universal significance cutoff."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Set up a series with a repeated motif\n",
+    "\n",
+    "The fixed random seed makes the example reproducible. The injected waveform is repeated at two known locations, while the rest of the series contains a smooth background and noise."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import stumpy\n",
+    "\n",
+    "rng = np.random.default_rng(2025)\n",
+    "n, m = 240, 24\n",
+    "time = np.arange(n)\n",
+    "series = 0.35 * np.sin(time / 18) + 0.08 * rng.normal(size=n)\n",
+    "motif = 0.8 * np.sin(np.linspace(0, 2 * np.pi, m))\n",
+    "for start in (58, 154):\n",
+    "    series[start : start + m] += motif\n",
+    "\n",
+    "plt.figure(figsize=(12, 3))\n",
+    "plt.plot(series, color=\"#3366aa\")\n",
+    "plt.axvspan(58, 58 + m, alpha=0.2, color=\"#dd8844\")\n",
+    "plt.axvspan(154, 154 + m, alpha=0.2, color=\"#dd8844\")\n",
+    "plt.title(\"Synthetic series with a repeated motif\")\n",
+    "plt.xlabel(\"Sample\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Measure the observed motif distance\n",
+    "\n",
+    "We use the minimum finite matrix-profile value as a deliberately simple test statistic. The profile index gives the location of the strongest candidate, but the value alone does not provide a p-value."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "profile = stumpy.stump(series, m)[:, 0]\n",
+    "finite_profile = profile[np.isfinite(profile)]\n",
+    "observed_distance = float(np.min(finite_profile))\n",
+    "observed_index = int(np.nanargmin(profile))\n",
+    "print(f\"candidate subsequence: {observed_index}\")\n",
+    "print(f\"observed minimum distance: {observed_distance:.4f}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Build a block-permutation null distribution\n",
+    "\n",
+    "Independent shuffling would destroy all local dependence. Instead, we split the series into circular blocks, rotate the block order, and rotate each block internally. This preserves short-range structure approximately while breaking the original motif alignment. The block size is a modeling choice, not a theorem."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def circular_block_permutation(values, block_size, random_generator):\n",
+    "    \"\"\"Return one reproducible circular block permutation.\"\"\"\n",
+    "    blocks = [values[i : i + block_size] for i in range(0, len(values), block_size)]\n",
+    "    random_generator.shuffle(blocks)\n",
+    "    rotated = [np.roll(block, random_generator.integers(len(block))) for block in blocks]\n",
+    "    return np.concatenate(rotated)[: len(values)]\n",
+    "\n",
+    "block_size = m\n",
+    "n_surrogates = 24\n",
+    "surrogate_minima = np.empty(n_surrogates)\n",
+    "for i in range(n_surrogates):\n",
+    "    surrogate = circular_block_permutation(series, block_size, rng)\n",
+    "    surrogate_profile = stumpy.stump(surrogate, m)[:, 0]\n",
+    "    surrogate_minima[i] = np.nanmin(surrogate_profile)\n",
+    "\n",
+    "# +1 avoids reporting an exact zero from a small finite sample.\n",
+    "p_value = (1 + np.count_nonzero(surrogate_minima <= observed_distance)) / (n_surrogates + 1)\n",
+    "percentile = 100 * (1 - p_value)\n",
+    "print(f\"empirical lower-tail p-value: {p_value:.3f}\")\n",
+    "print(f\"observed distance is above {percentile:.1f}% of the null minima\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure(figsize=(8, 3))\n",
+    "plt.hist(surrogate_minima, bins=8, color=\"#99bbee\", edgecolor=\"white\")\n",
+    "plt.axvline(observed_distance, color=\"#cc5533\", linewidth=2, label=\"observed minimum\")\n",
+    "plt.xlabel(\"minimum matrix-profile distance\")\n",
+    "plt.ylabel(\"surrogate count\")\n",
+    "plt.title(\"Observed statistic against the block-permutation null\")\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Check sensitivity to the number of surrogates\n",
+    "\n",
+    "The estimate is intentionally coarse with 24 surrogates. Repeating the calculation with more surrogates should make the histogram and percentile less variable, at a higher runtime cost. For a real analysis, choose the count before inspecting the result and report the seed, block size, and number of surrogates."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretation and limitations\n",
+    "\n",
+    "* Overlapping subsequences are dependent, so this is not a collection of independent tests.\n",
+    "* Searching every profile index is a multiple-comparisons problem; using the minimum statistic partly reflects that search but does not solve every inference issue.\n",
+    "* The null distribution depends on the block size and permutation scheme. Different choices answer different questions.\n",
+    "* The empirical p-value is a descriptive comparison, not a guarantee that the motif is causal, novel, or useful. Validate any threshold on domain-specific data."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"},
+  "language_info": {"name": "python", "version": "3.10"}
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/tutorials.rst

@@ -27,3 +27,4 @@ Tutorials
     Tutorial_Multidimensional_Motif_Discovery.ipynb
     Tutorial_Annotation_Vectors.ipynb
     Tutorial_Shapelet_Discovery.ipynb
+    Tutorial_Statistical_Significance.ipynb

notebooks/Tutorial_Statistical_Significance.ipynb

@@ -0,0 +1 @@
+../docs/Tutorial_Statistical_Significance.ipynb