218 lines
92 KiB
Plaintext
218 lines
92 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "7608f93c",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# 16 — Robust drift estimation"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 1,
|
||
|
|
"id": "5281ce54",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-05-12T10:16:20.462286Z",
|
||
|
|
"iopub.status.busy": "2026-05-12T10:16:20.461994Z",
|
||
|
|
"iopub.status.idle": "2026-05-12T10:16:21.152375Z",
|
||
|
|
"shell.execute_reply": "2026-05-12T10:16:21.150556Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"from optimizr import _core as opt\n",
|
||
|
|
"plt.rcParams['figure.figsize'] = (7, 4)\n",
|
||
|
|
"plt.rcParams['figure.dpi'] = 110\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "c351399a",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Synthetic stationary process with 5 % outliers"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"id": "08c7105d",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-05-12T10:16:21.156097Z",
|
||
|
|
"iopub.status.busy": "2026-05-12T10:16:21.155724Z",
|
||
|
|
"iopub.status.idle": "2026-05-12T10:16:21.211614Z",
|
||
|
|
"shell.execute_reply": "2026-05-12T10:16:21.210552Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"observation length = 5001\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"rng = np.random.default_rng(7)\n",
|
||
|
|
"true_a, true_b = 1.0, -0.5\n",
|
||
|
|
"dt, n = 0.01, 5000\n",
|
||
|
|
"x = [0.0]\n",
|
||
|
|
"for k in range(n):\n",
|
||
|
|
" if k % 20 == 0:\n",
|
||
|
|
" eps = rng.uniform(-2.0, 2.0)\n",
|
||
|
|
" else:\n",
|
||
|
|
" eps = rng.uniform(-0.1, 0.1)\n",
|
||
|
|
" x.append(x[-1] + (true_a + true_b * x[-1]) * dt + eps * np.sqrt(dt))\n",
|
||
|
|
"x = np.array(x)\n",
|
||
|
|
"print('observation length =', len(x))\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 3,
|
||
|
|
"id": "4fc130d8",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-05-12T10:16:21.214664Z",
|
||
|
|
"iopub.status.busy": "2026-05-12T10:16:21.214284Z",
|
||
|
|
"iopub.status.idle": "2026-05-12T10:16:21.507788Z",
|
||
|
|
"shell.execute_reply": "2026-05-12T10:16:21.506501Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvYAAAGtCAYAAAB9QDCJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAABDrAAAQ6wFQlOh8AAC+zUlEQVR4nOzdd5gT1foH8O9MspvtvbBshwUWpJelL1WqiiJ2Eexee7lesAOKgCK2K+pVbFyUa/cnICBVkCZFeoeFZXvPluxukjm/PyaTTTZ9N9lkk/fzPD6yk0nmZCaZvHPmPe/hGGMMhBBCCCGEkHaNd3cDCCGEEEIIIa1HgT0hhBBCCCFegAJ7QgghhBBCvAAF9oQQQgghhHgBCuwJIYQQQgjxAhTYE0IIIYQQ4gUosCeEEEIIIcQLUGBPCCGEEEKIF6DAnhBCCCGEEC9AgT0hhBBCCCFegAJ7QtopjuMwe/bsNt/u6NGjkZaW1ubbdaZt27aB4zh88cUX7m6Kw+bNmweO45CTk2PX+s56r95w3N3N3Hc2LS0No0ePdsn2HHltd51PPBntE9IeUWBPiBnFxcWYO3cuevXqhbCwMISGhqJTp0644YYbsGLFijZpQ2VlJebNm4dt27a1yfYMzZs3Dz///HObb5e0zN9//4158+bZHewT+23btg3z5s1DZWWlu5tCnCwnJwfz5s3D33//7e6mEOI0cnc3gBBPc/nyZWRlZaG0tBQzZszA/fffD39/f1y4cAE7d+7EO++8g3vvvdfl7aisrMT8+fMBwGU9epbMnz8fs2bNwvXXX2/y2MaNG8EYa9P2OFt2djZUKhX8/Pzc3RSHvfjii5g7dy4UCoV+2d9//4358+dTr7oLbNu2DfPnz8fs2bMRERHRqtdSqVSQyWTOaZiTeXLbXCUnJwfz589HWloa+vbta/K4L+4T0v5RYE9IM2+++SaKiorwzjvv4IknnjB5vLCw0A2t8hz+/v7ubkKLKZVKhIWFged5BAQEuLs5LSKXyyGX06m7PfLkz5wnt81daJ+Q9ohScQhp5uzZswCAcePGmX28Q4cO+n8/+uij4DgOJ06cMFmvsbERsbGxGDhwoH6ZlPN65swZTJs2DeHh4QgJCcGUKVNw7tw5/XpffPEF0tPTAYi95xzHgeM4s72x+/btw9ixYxESEoKIiAjceuutKC4uNtueN954A71790ZgYCDCwsIwfvx4/PHHH/p1pHxsAPjyyy/125WWAZZzrS9evIj7778fqampUCgUiI+Px4QJE/D777+b3Y+GTp48idtuuw3JyclQKBSIi4vDsGHD8OmnnxqtxxjDJ598gqysLAQHByM4OBjDhg0zmzYk5cdu27YNo0ePRlhYGPr06WP0PpvnnTvy+uvXr8fYsWMRFxeHgIAAJCUlYfLkydixY4fV9/r111+D4zj8+uuvRst79OgBjuOwcuVKo+XDhg1Dp06d9H83z7GfPXs27r77bgDAmDFj9MfLXG7wypUr0bt3bwQEBCAxMRHPP/88tFqt1fY2V1hYiJkzZyI6OhqBgYHIzs7G/v37za77ww8/YNSoUQgLC0NgYCD69etnckwB8S7Qbbfdhs6dO+s/m9nZ2Sb76MUXXwTHcfjzzz/Nbq9Lly5ITU2FIAgOfzfNGT16tP6uWXp6un7fzps3DwBQXV2Nl156CUOGDEFsbCz8/f2RlpaGRx99FOXl5Sav50jO9qFDhzBjxgzExcXB398fnTp1wty5c1FXV2ey7oEDBzB+/HgEBwcjMjISN954o8NpWebaJi2z5xzzxRdfgOM4bN26Fe+88w66du0KhUKB9PR0LFu2zOw29+7di2uuuQZRUVEICAhAZmYmXn31VTQ2NurX+fTTT8FxHFatWmX2NaT3rVQqAQCnTp3CI488gp49eyI8PByBgYHo1asXli5davRZnzdvHsaMGQMAuPvuu/XH1vDuqKXjtWrVKgwePFh/jhgyZAhWr15tsp50rrT3O7Nq1SoMHToUUVFRCAwMREpKCqZPn272M0yIRYwQYuThhx9mANjTTz/N1Gq11XWPHDnCALCnnnrK5LHVq1czAOzjjz/WL0tNTWUZGRksPj6e3X///ezDDz9kTz/9NPPz82M9evRgWq2WMcbY+fPn2dtvv80AsBtuuIGtXLmSrVy5kv3000/61wLA+vbty6KiotiTTz7JPvroI3b//fczjuPYxIkTjdqiVqvZuHHjmFwuZzNnzmQffPABe/PNN1mfPn2YTCZjv/76K2OMscLCQrZy5UoGgI0cOVK/3ZUrV+pfa9SoUSw1NdXo9Q8cOMAiIiKYv78/e+CBB9jy5cvZkiVL2A033MD+9a9/Wd2HpaWlLD4+nsXGxrIXX3yRrVixgi1dupTNmjWL3XnnnUbrzpo1i3Ecx66//nr2zjvvsHfeeYdlZ2czAOzDDz80WhcAu+qqq1hwcDB77LHH2Mcff8yWLl3KGGNs69atDAD7/PPPW/T627dvZzKZjPXq1Yu98cYbbMWKFWzhwoVs6tSp7N1337X6fgsLCxkA9vjjj+uX5eXlMQCM53l211136ZdXVVUxuVzO7r//fv2yV155hQFgFy9eZIwxtmvXLvbAAw8wAOz555/XH69du3YZvdchQ4aw5ORkNn/+fLZ8+XI2fvx4BoAtWrTIanslo0aNYjExMSwjI4Pdcsst7IMPPmAvv/wyCw0NZTExMUypVBqt//LLLzMAbMyYMWzp0qXsgw8+YNdddx0DwObMmWO07m233cbGjBnDXn75Zfaf//yHLVy4kHXt2pUBYKtXr9avd/bsWQaA3XfffSbt27FjBwPAXnzxRcaY499NczZu3MhuuOEGBoC9/fbb+n17+PBhxhhjJ0+eZHFxceyhhx5iy5YtY8uXL2ezZ89mcrmc9enThzU2Nhq9HgA2a9Yso2Wpqals1KhRRst+++03plAoWEZGBluwYAH7+OOP2cMPP8z8/PzYiBEjjM5LBw4cYEFBQSw0NJTNmTOH/fvf/2bXXXcdS01NZTExMSavbYm5tjlyjvn888/1n7MePXqw119/nb333nts4MCBDAD75ptvTN6jn58fi46OZi+88AJ7//332eTJkxkANmnSJP25sKqqigUFBbHx48ebtPny5cuM53mj88SHH37Iunfvzp599lm2fPlytmzZMv1n/R//+Id+vcOHD7Pnn3+eAWAPPPCA/thu3LjR6j556aWXGADWq1cvtmTJErZ48WLWs2dPBoAtXLjQaF1HvjP//e9/GQA2fPhw9s4777BPP/2UzZ8/n40aNYr9+OOPVo4cIcYosCekmfPnz7Pw8HAGgMXFxbEbb7yRLVmyhO3cuVP/Y2No2LBhLDo6mtXX1xstHzduHAsJCWHV1dX6ZampqQwA+/rrr43WXbRoEQPANmzYoF928eJFBoC98sorZtsJgHEcx/7880+j5Q8++CADwE6fPq1f9s477zAAJj8QjY2NrF+/fiw9Pd3ktZv/oEmaB/aCILCePXsyuVzO9u7da7K+uX1m6JdffjEJ4Mz5+eefGQC2bNkyk8euvfZaFhYWZvRDCYABYL/99pvJ+uYCe0de/6mnnmIAWGFhodU2W9KzZ0/Wo0cP/d9ffvkl4ziOzZw5kyUmJuqXm9s3zQN7xpqCqq1bt1p8rx06dGDl5eX65VqtlnXv3p0lJCTY1eZRo0YxAOz11183Wv7NN9+YBMkHDx5kHMcZXbxIHn30UcbzPDt//rx+WU1Njcl6tbW1rEuXLkb7iTHGRowYwcLCwlhdXZ3R8nvvvZcBYOfOndMvc+S7aYm5/S1paGgwCd4ZY+yTTz5hANi3335rtNyewF6lUrEOHTqwrKwsk3Z///33DAD74osv9MtGjhzJeJ5n+/fvN1pXOg+0NrC39xwjfQZ79+5t1O6amhoWHR3Nhg4dql+m0WhYWloaCwwMZGfPnjV67bvvvpsBMOpMuPPOOxnP8+zy5ctG67766qsMANu8ebPR9sy5/fbbmUwmYwUFBfplli7wLe2TM2fOMJ7nWZ8+fVhtba3RNnv27MlkMpnR58SR78wNN9zAQkNDzX6
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 770x440 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"fig, ax = plt.subplots()\n",
|
||
|
|
"ax.plot(x, lw=0.6)\n",
|
||
|
|
"ax.axhline(true_a / -true_b, color='red', ls='--', label='OU level a/(-b) = 2')\n",
|
||
|
|
"ax.set_xlabel('k'); ax.set_ylabel('x_k'); ax.legend(); ax.grid(alpha=0.3)\n",
|
||
|
|
"ax.set_title('Synthetic series with heavy-tailed innovations')\n",
|
||
|
|
"fig.tight_layout(); plt.show()\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 4,
|
||
|
|
"id": "9ce9d2b4",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-05-12T10:16:21.510592Z",
|
||
|
|
"iopub.status.busy": "2026-05-12T10:16:21.510324Z",
|
||
|
|
"iopub.status.idle": "2026-05-12T10:16:21.518890Z",
|
||
|
|
"shell.execute_reply": "2026-05-12T10:16:21.516422Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"a (true 1.0) -> 0.9402\n",
|
||
|
|
"b (true -0.5) -> -0.4721\n",
|
||
|
|
"IRLS iterations = 5\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"res = opt.robust_drift(x.tolist(), dt=dt)\n",
|
||
|
|
"print(f'a (true 1.0) -> {res[\"a\"]:.4f}')\n",
|
||
|
|
"print(f'b (true -0.5) -> {res[\"b\"]:.4f}')\n",
|
||
|
|
"print('IRLS iterations =', res['iterations'])\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 5,
|
||
|
|
"id": "55e42e3b",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-05-12T10:16:21.521983Z",
|
||
|
|
"iopub.status.busy": "2026-05-12T10:16:21.521704Z",
|
||
|
|
"iopub.status.idle": "2026-05-12T10:16:21.755541Z",
|
||
|
|
"shell.execute_reply": "2026-05-12T10:16:21.754070Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"OLS a, b = [ 1.00031819 -0.51377951]\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvYAAAGtCAYAAAB9QDCJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAABDrAAAQ6wFQlOh8AABE6ElEQVR4nO3deXhTZf738U+a0oWuoa1QUHaEAhUQqcDosIOoYEVkUbQVi1YGZFFQ8AcUGRZHROY3qAzKA7Ipi4VxAYSyOIOAODqOChQFlbVKKV1AaKHpef7wIQ8hbWmxSdrD+3VdXJo75z7nmybf5tPTO6cWwzAMAQAAAKjSfLxdAAAAAIDfj2APAAAAmADBHgAAADABgj0AAABgAgR7AAAAwAQI9gAAAIAJEOwBAAAAEyDYAwAAACZAsAcAAABMgGAPAAAAmADBHkCltn37dlksFi1evNjbpcBNUlJSZLFY9NNPP5V5Tv369dW5c2eX8TVr1qhVq1YKDAyUxWLR9u3bK6xOd0pMTJTFYvF2GQCqOII9gN/tUvi+/F/16tUVGxurF198UefPn/d2ib/L9u3blZKSopycHG+XUqKLFy9qwYIF6tSpkyIiIuTv76+bbrpJgwcP1meffVbsnM6dO8vX17dM+//00091//33q2HDhgoICFBUVJRatWqlJ598Uv/5z38q8qFcs++++06DBw9WWFiY5s2bp6VLlyomJkaLFy/W3LlzvV2e5s6dW+l/QP3qq6+UkpJSrh+yAFQeZfuODgBl0L9/f913332SpMzMTK1atUpTpkzRzp07tXHjRi9Xd+22b9+uqVOnKjExUeHh4d4ux8WpU6fUp08f7d69W3fccYcmTJggm82mgwcPasmSJXr33Xf14osvatKkSde0/7///e9KTk5WrVq19Oijj6px48bKycnRd999p48++khNmjRRmzZtKvhRle7AgQMuZ7i3b9+uwsJCzZ07V7feeqtjfPHixfrpp580evRoj9Z4pblz56p+/fpKTEx0ue/NN9/U/PnzPV/UFb766itNnTpVnTt3Vv369b1dDoByItgDqDCtWrXSkCFDHLeffvppxcXF6eOPP9YXX3yhtm3berE68xo4cKB2796t2bNn65lnnnG6b+LEibr33ns1efJk1a9fX4888ki59l1YWKgJEyYoODhYn3/+uW688Uan+4uKipSVlfW7H0NZXLhwQUVFRQoICJC/v7/L/T///LMkqUaNGh6ppyJVq1ZN1apV83YZAKo4luIAcBur1aouXbpIkr7//nun+44fP66kpCTVqVNHfn5+uvHGG/XEE08oIyOjxP298cYbiomJUUBAgOrXr6+UlBQVFhY6bVPamUaLxeJytnT58uXq0KGDatSoocDAQNWtW1f9+vXTvn37HPubOnWqJKlBgwaOpUYpKSkl1vnWW2/JYrFo+fLlxd7fvXt3BQUFKS8vT5J07NgxPfHEE2rQoIECAgIUGRmptm3basaMGSUe45L169dr69atuu+++1xCvSSFhITo3XffVWBgoCZMmKCLFy9edZ+XO3XqlLKzs9W0aVOXUC9JPj4+ioqKKtO+zpw5o1GjRik6OlqBgYG69dZbtXr16mK3vbTmPCsrS0888YRjzu7duyU5r7H/6aefZLFYNGXKFEn//3mqX7++LBaLPvnkEx0+fNhpqVhZ1t4fOnRIiYmJql27tuM1Onz4cJ06dcppu+zsbI0bN05NmjRRYGCgbDabYmNjHb8huFTf4cOH9cknnzjVcWnJS3Fr7C+NnT59WklJSbrhhhsUHBysnj17Ovrp/fffV7t27VS9enXVqVNHM2fOdHkce/bs0dChQ9W0aVMFBQUpKChI7dq106JFi1yO99hjj0mSunTp4qjx8p4xDENvvvmm4uLiHPvq2LGj1q1bd9WvJwD344w9ALc6dOiQJCkiIsIxdvz4cbVr104nT55UUlKSWrVqpf/+97968803tXHjRn3++eeqWbOm037mzZunY8eOKTk5WTVq1NA//vEPTZ06VYcOHdLSpUuvqbbly5dryJAh+sMf/qApU6YoODhYx48f19atW3XgwAE1b95cL7zwgmrUqKG1a9fq1VdfVWRkpCTplltuKXG/AwYM0KhRo7R48WI9/PDDTvcdPXpU27Zt00MPPaTQ0FAVFhaqR48eOnr0qJ566ik1a9ZMZ8+eVXp6urZu3aqJEyeW+hguBeOnnnqqxG2io6MVHx+vd955R7t27dIf//jHsn6JVLNmTQUHB2vv3r3auXOnOnbsWOa5lyssLFTv3r0da/W7deumI0eOaOjQobr55ptLnNe9e3dFRETo+eefV1FRkWrVquWyTVRUlJYuXarU1FSn5yk4OFhnz57V9OnTderUKb366quOOTExMaXW+9VXX6lz586qXr26hg4dqnr16un777/XG2+8oS1btmjPnj0KCwuT9NvzvW3bNj3xxBNq3bq1Lly4oEOHDiktLc2pvjFjxigyMlIvvPCCU+1Xc9ddd6lmzZqaMmWKTpw4oTlz5qhnz56aNm2annnmGSUnJ+uxxx7Tu+++q4kTJ6p+/foaPHiwY/7atWv17bffqn///qpXr55yc3O1atUqDR06VJmZmRo/frwk6cknn5S/v78WLFigiRMnOr5GjRo1cuzrscce05IlS3Tfffc5Xtupqam6//779cYbbyg5OfmqjweAGxkA8Dtt27bNkGRMmDDByMzMNDIzM419+/YZkyZNMiQZ9erVMwoKChzbP/LII4YkY/ny5U77efvttw1JxuOPP+6y7+rVqxs//fSTY9xutxvx8fGGJGPbtm2O8U6dOhn16tUrtk5JRkJCguP2/fffb4SEhBgXLlwo9fFNmTLFkGT8+OOPV/9i/D9DhgwxfHx8jCNHjjiNT5s2zZBkbNmyxTAMw/jvf/9rSDJmzZpV5n1frm3btoYkIysrq9TtZs+ebUgy/va3vznGOnXqZFit1qse49JcSUZsbKyRnJxsLFy4sFxfj4ULFxqSjFGjRjmN79y507BYLC5f34SEBEOSMWjQIKOoqMhlf/Xq1TM6derkNFbS81Taa6IkrVu3Nho0aODydf3ss88Mq9VqpKSkGIZhGDk5OYYkIzk5+ar7LK7mSy493uLGnnzySafxV1991ZBkBAcHOz3W/Px8o2bNmkaHDh2ctj979qzL8ex2u3HnnXcaYWFhTq//RYsWufTUJevWrTMkGXPmzHG5r0+fPkZoaKiRl5dX7OMD4BksxQFQYWbOnKmoqChFRUWpefPmmjZtmnr27Km0tDT5+flJ+m1N9rp169S0aVM99NBDTvMfeeQRNWrUSKmpqTIMw+m+IUOGqF69eo7bPj4+mjBhgiTpvffeu6Z6w8PDde7cOX3wwQcqKiq6pn2UJDExUUVFRVqyZInT+Ntvv6169eo5lihdOuu7bds2xxrx8sjNzXXaT0ku3X9p+/J45pln9OGHH+ree+/V4cOHNX/+fD3++ONq0KCB7rvvPmVmZl51H5eeoyt/A9GhQwd169atxHnPPfecxy8D+e233+qrr77SoEGDVFRUpFOnTjn+NWzYUI0bN9bHH38sSQoMDFRAQIA+++wz/fDDD26p58olVp06dZIk9e3b12nZmb+/v26//XZ99913TtsHBQU5/v/8+fPKysrS6dOndddddyk3N1cHDhwoUx1Lly5VYGCgBg4c6PQ1OXXqlOLj45WXl6ddu3Zd46MEUBEI9gAqTGJiojZv3qwNGzZo7ty5io6O1rFjxxQYGOjYJjMzU2fOnFHLli1d5lssFrVo0ULZ2dnKzs52uq958+Yu218aO3jw4DXV+8ILL6hhw4Z64IEHFBkZqT59+ujVV1/VL7/8ck37u1zXrl1Vr149vf32246xHTt26ODBg3r00UcdYbVevXqaMmWKNm/erNq1a6tVq1b605/+pM2bN5fpOKGhoZKuHtjL+gNASe655x598MEHys7O1r59+/T666+rRYsWev/9950+MF2SQ4cOKTIyUjfccIPLfS1atChxXmnLdNxl//79kpx/UL3834EDBxyvET8/P/3v//6v9u3
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 770x440 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Compare against a naïve OLS that is broken by outliers.\n",
|
||
|
|
"y = (x[1:] - x[:-1]) / dt\n",
|
||
|
|
"X = np.vstack([np.ones_like(x[:-1]), x[:-1]]).T\n",
|
||
|
|
"ols_ab, *_ = np.linalg.lstsq(X, y, rcond=None)\n",
|
||
|
|
"print('OLS a, b =', ols_ab)\n",
|
||
|
|
"fig, ax = plt.subplots()\n",
|
||
|
|
"labels = ['true', 'OLS', 'robust']\n",
|
||
|
|
"vals_a = [true_a, ols_ab[0], res['a']]\n",
|
||
|
|
"vals_b = [true_b, ols_ab[1], res['b']]\n",
|
||
|
|
"ax.bar(np.arange(3) - 0.2, vals_a, width=0.4, label='a')\n",
|
||
|
|
"ax.bar(np.arange(3) + 0.2, vals_b, width=0.4, label='b')\n",
|
||
|
|
"ax.set_xticks(range(3)); ax.set_xticklabels(labels)\n",
|
||
|
|
"ax.legend(); ax.grid(alpha=0.3); ax.set_title('Robust vs OLS drift estimate')\n",
|
||
|
|
"fig.tight_layout(); plt.show()\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "70fbd189",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"**Verified:** Huber IRLS recovers `(a, b)` within `0.2` even with 5 % heavy outliers."
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3 (rhftlab)",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"codemirror_mode": {
|
||
|
|
"name": "ipython",
|
||
|
|
"version": 3
|
||
|
|
},
|
||
|
|
"file_extension": ".py",
|
||
|
|
"mimetype": "text/x-python",
|
||
|
|
"name": "python",
|
||
|
|
"nbconvert_exporter": "python",
|
||
|
|
"pygments_lexer": "ipython3",
|
||
|
|
"version": "3.11.13"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|