[DOC] Add self-similarity and Hurst exponent example to ANM Lec 2

lectures/lec2/Time_Domain_Features.ipynb

@@ -408,6 +408,74 @@
     "plt.tight_layout()\n",
     "plt.show()\n"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Self-Similarity and Hurst Exponent\n",
+    "\n",
+    "Self-similarity means similar statistical behavior after rescaling time and amplitude.\n",
+    "This example generates a fractional Brownian motion (fBm) and estimates H from\n",
+    "the scaling of its increment standard deviation. For fBm increments, H > 0.5\n",
+    "indicates persistence, H = 0.5 no correlation, and H < 0.5 anti-persistence.\n",
+    "H is a scaling exponent, not a self-similarity score. This estimator is for\n",
+    "fBm-like paths; arbitrary or stationary signals require checking the model first.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "metadata": {},
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Generate a fractional Brownian motion (fBm) with unit time steps\n",
+    "def generate_self_similar(n=8192, H=0.75, seed=42):\n",
+    "    if not isinstance(n, int) or n < 16 or not 0 < H < 1:\n",
+    "        raise ValueError(\"Use integer n >= 16 and 0 < H < 1.\")\n",
+    "    k = np.arange(n + 1, dtype=float)\n",
+    "    covariance = 0.5 * (np.abs(k - 1)**(2*H) - 2*k**(2*H) + (k + 1)**(2*H))\n",
+    "    eigenvalues = np.fft.fft(np.r_[covariance, covariance[-2:0:-1]]).real\n",
+    "    if eigenvalues.min() < -1e-10:\n",
+    "        raise ValueError(\"Invalid covariance embedding; increase n or adjust H.\")\n",
+    "    noise = np.random.default_rng(seed).standard_normal(2 * n)\n",
+    "    increments = np.fft.ifft(np.sqrt(np.maximum(eigenvalues, 0)) * np.fft.fft(noise)).real[:n]\n",
+    "    return np.r_[0, np.cumsum(increments)]\n",
+    "\n",
+    "# For fBm-like paths: std(x[t + lag] - x[t]) is proportional to lag**H\n",
+    "def estimate_hurst(x):\n",
+    "    x = np.asarray(x, dtype=float)\n",
+    "    if x.ndim != 1 or len(x) < 256 or not np.all(np.isfinite(x)):\n",
+    "        raise ValueError(\"Use a finite 1-D fBm-like path with at least 256 points.\")\n",
+    "    lags = np.unique(np.geomspace(2, min(256, len(x) // 16), 20).astype(int))\n",
+    "    fluctuation = np.array([np.std(x[lag:] - x[:-lag]) for lag in lags])\n",
+    "    if np.any(fluctuation <= 0):\n",
+    "        raise ValueError(\"The series must have nonzero fluctuations.\")\n",
+    "    H, intercept = np.polyfit(np.log(lags), np.log(fluctuation), 1)\n",
+    "    return H, lags, fluctuation, np.exp(intercept) * lags**H\n",
+    "\n",
+    "# Generate data and estimate H\n",
+    "true_H = 0.75\n",
+    "self_similar_series = generate_self_similar(H=true_H)\n",
+    "estimated_H, lags, fluctuation, fitted = estimate_hurst(self_similar_series)\n",
+    "print(f\"True H: {true_H:.2f}; Estimated H: {estimated_H:.3f}\")\n",
+    "\n",
+    "# Plot the time series and its scaling relationship\n",
+    "fig, axs = plt.subplots(2, 1, figsize=(15, 10))\n",
+    "axs[0].plot(self_similar_series, color='blue')\n",
+    "axs[0].set(title=f'Self-Similar Time Series (H = {true_H})', xlabel='Time', ylabel='Value')\n",
+    "axs[1].loglog(lags, fluctuation, 'o', label='Increment standard deviation')\n",
+    "axs[1].loglog(lags, fitted, '--', label=f'Estimated H = {estimated_H:.3f}')\n",
+    "axs[1].set(title='Hurst Exponent Estimation', xlabel='Time Lag', ylabel='Fluctuation')\n",
+    "axs[1].legend()\n",
+    "for ax in axs:\n",
+    "    ax.grid(True)\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n"
+   ],
+   "execution_count": null,
+   "outputs": []
   }
  ],
  "metadata": {