{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# V0.1.6 - Important notes and examples of how to use Extended Least Squares\n",
"\n",
"Example created by Wilson Rocha Lacerda Junior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we import the NARMAX model, the metric for model evaluation and the methods to generate sample data for tests. Also, we import pandas for specific usage."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pip install sysidentpy"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sysidentpy.polynomial_basis import PolynomialNarmax\n",
"from sysidentpy.metrics import root_relative_squared_error\n",
"from sysidentpy.utils.generate_data import get_miso_data, get_siso_data\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating 1 input 1 output sample data \n",
"\n",
"The data is generated by simulating the following model:\n",
"$y_k = 0.2y_{k-1} + 0.1y_{k-1}x_{k-1} + 0.9x_{k-2} + e_{k}$\n",
"\n",
"If *colored_noise* is set to True:\n",
"\n",
"$e_{k} = 0.8\\nu_{k-1} + \\nu_{k}$\n",
"\n",
"where $x$ is a uniformly distributed random variable and $\\nu$ is a gaussian distributed variable with $\\mu=0$ and $\\sigma$ is defined by the user.\n",
"\n",
"In the next example we will generate a data with 3000 samples with white noise and selecting 90% of the data to train the model. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"x_train, x_valid, y_train, y_valid = get_siso_data(n=1000,\n",
" colored_noise=True,\n",
" sigma=0.2,\n",
" train_percentage=90)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build the model\n",
"\n",
"First we will train a model without the Extended Least Squares Algorithm for comparison purpose."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"model = PolynomialNarmax(non_degree=2,\n",
" order_selection=True,\n",
" n_info_values=10,\n",
" n_terms=3,\n",
" extended_least_squares=False,\n",
" ylag=2, xlag=2,\n",
" info_criteria='aic',\n",
" estimator='least_squares',\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.5626926249339153\n"
]
}
],
"source": [
"model.fit(x_train, y_train)\n",
"yhat = model.predict(x_valid, y_valid)\n",
"rrse = root_relative_squared_error(y_valid, yhat)\n",
"print(rrse)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly we have something wrong with the obtained model. See the *basic_steps* notebook to compare the results obtained using the same data but without colored noise. But let take a look in whats is wrong."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Regressors Parameters ERR\n",
"0 x1(k-2) 0.9094 0.74535677\n",
"1 y(k-1) 0.2794 0.07821222\n",
"2 x1(k-1)y(k-1) 0.1188 0.00443378\n"
]
}
],
"source": [
"results = pd.DataFrame(model.results(err_precision=8,\n",
" dtype='dec'),\n",
" columns=['Regressors', 'Parameters', 'ERR'])\n",
"\n",
"print(results)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Biased parameter estimation\n",
"\n",
"As we can observe above, the model structure is exact the same the one that generate the data. You can se that the ERR ordered the terms in the correct way. And this is an importante note regarding the Error Reduction Ratio algorithm used here: __it is very robust to colored noise!!__ \n",
"\n",
"That is a greate feature! However, although the structure is correct, the model *parameters* are not ok! Here we have a biased estimation! The real parameter for $y_{k-1}$ is $0.2$, not $0.3$.\n",
"\n",
"In this case, we are actually modeling using a NARX model, not a NARMAX. The MA part exists to allow a unbiased estimation of the parameters. To achieve a unbiased estimation of the parameters we have the Extend Least Squares algorithm. Remember, if the data have only white noise, NARX is fine. \n",
"\n",
"Before applying the Extended Least Squares Algorithm we will run several NARX models to check how different the estimated parameters are from the real ones."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD7CAYAAACPDORaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAcsElEQVR4nO3dbWwU17kH8P/MvhmD23Cd3bji+nKl3DRWSQuo0lVoKyOqYkjB0DhRG0DwoS8QpQVBJJoUXJCaEihFpYqaSI1UqepLpFAkMEWpSVsapFyii8oHEJQ2qcD0hoC9vCT4bV9m59wPuzM7Yxt2196ZOcz5/6rKu+tl9/gsefzwnOec0YQQAkREFEp60AMgIiLvMMgTEYUYgzwRUYgxyBMRhRiDPBFRiDHIExGFGIM8EVGIRYMewFi3bg3DNKfeuv/cc88CAH70o59M+bXq5bnnnkUsFsEPf/jjoIcivebmGbhxYyjoYdwTOFfVkW2e6hWjdF3DzJnT7/h96YK8aYq6BPmBgQH79WQxMDCAeDwq1ZhkxnmqHueqOjLNk18xiuUaIqIQY5AnIgoxBnkiohBjkCciCjEGeSKiEGOQJyIKMQZ5IqIAfDScQ8GHlk7p+uSJiMIu/eEort0YQWOD9yGYmTwRkc+Mgun66iUGeSIin0V0DQDgx8VXGeSJiHyma5p/7+XbOxEREQBAjPnqJQZ5IiKfCatO40O9hkGeiMhnzOSJiELMjwVXC4M8EZHPhI9RnkGeiMhnzOSJiEKMmTwRUYgxkyciCjGTmTwREdUDgzwRkc9YriEiCjHpyjU/+9nPsGzZMixbtgx79+4FAJw8eRKdnZ3o6OjA/v377edeuHABXV1dWLJkCbZv3w7DMLwZORERVVQxyJ88eRJvv/02Dh06hMOHD+P8+fM4evQotm3bhldeeQVvvPEGzp07hxMnTgAAtm7dih07duDYsWMQQuDAgQOe/xBERPcSqTL5ZDKJ559/HvF4HLFYDA8++CD6+vowe/ZstLa2IhqNorOzE729vbhy5QoymQzmzZsHAOjq6kJvb6/nPwQR0b1Eqpr8Qw89ZAftvr4+/OEPf4CmaUgmk/ZzUqkU+vv7MTAw4Ho8mUyiv7/fg2ETEd27/NwMVfUFBt977z1s2LAB3/3udxGJRNDX12d/TwgBTdNgmiY0x2H41uO1aG6eUdPz7yQeL/5oyWRTXV6vHmQck8w4T9XjXFVHlnkaGMzZt70eU1VB/vTp09i0aRO2bduGZcuW4dSpU0in0/b30+k0UqkUWlpaXI9fv34dqVSqpgHduDEEsw5XMM/ljNLYBqf8WvWSyxmIx6NSjUlWyWQT56lKnKvqyDRPt24N27enOiZd1+6aHFcs11y9ehXf/va3sW/fPixbtgwAMHfuXFy6dAmXL19GoVDA0aNH0d7ejlmzZiGRSOD06dMAgJ6eHrS3t0/pByAiChs/a/IVM/lf/OIXyGaz2LNnj/3YU089hT179mDjxo3IZrNYuHAhli5dCgDYt28furu7MTQ0hDlz5mDdunXejZ6I6B7kY4yvHOS7u7vR3d094feOHDky7rG2tjYcPHhw6iMjIgopnkJJRBRiUrVQEhFRfQkfCzYM8kREPmMmT0QUYqzJExGFWB22AlWNQZ6IyG8M8kRE4cVyDRFRiLFcQ0QUaszkiYhCiy2UREQhJtWVoYiIqL6YyRMRhRiPNSAiCjFm8kREIcY+eSKiEGMmT0QUYgzyREQhxnINEVGI+XmNVwZ5IiKfMZMnIgox1uSJiEKMmTwRUYixJq+A7OnDyF88FfQwiCgAfpZrov69FTnlTh8GAMTW/3fAIyEiv/EUSiKiMOPCKxFReHHhlYgoxHiNVyIiqgsGeSIin7FcE3J+fsBEJB+ru0bz4b0Y5IPAIE+kNB5rEHaiEPQIiChA3PEadsIMegREFCDW5MPOZJAnUhnLNWHHTJ5IaczkQ06YrMkTqcyK8X6Eegb5IDCTJ1Ka8HHptaogPzQ0hOXLl+P9998HAHzve99DR0cHVq5ciZUrV+KPf/wjAODChQvo6urCkiVLsH37dhiG4d3I72WsyRMpTaqa/JkzZ7Bq1Sr09fXZj507dw6/+c1v0NPTg56eHixevBgAsHXrVuzYsQPHjh2DEAIHDhzwbOD3NGbyREqzg7wPwb5ikD9w4AB27tyJVCoFABgdHcUHH3yAbdu2obOzEy+99BJM08SVK1eQyWQwb948AEBXVxd6e3u9Hf29ikGeSGl+LrxWvGjIrl27XPevX7+ORx99FDt37kRTUxM2bNiAgwcP4qGHHkIymbSfl0wm0d/fX/8RhwHLNURK83PhteYrQ7W2tuLll1+2769duxaHDx/Ggw8+CE0rn8QghHDdr1Zz84ya/8xE4vHij5ZMNtXl9erBGtPMmdMwXHpMpvHJhnNTPc5VdWSZp4ZpMfu212OqOcj/4x//QF9fH5YsWQKgGMyj0ShaWlqQTqft512/ft0u8dTixo0hmHU4bDmXKy76ptODU36tesnlDMTjUdy8cdt+TKbxySSZbOLcVIlzVR2Z5mlkJGffnuqYdF27a3JccwulEAIvvvgiPvroI+Tzebz++utYvHgxZs2ahUQigdOnTwMAenp60N7ePvmRhxnLNURK8/MarzVn8m1tbVi/fj1WrVoFwzDQ0dGB5cuXAwD27duH7u5uDA0NYc6cOVi3bl3dBxwKXHglUptVk/ch2Fcd5I8fP27fXrNmDdasWTPuOW1tbTh48GB9RhZmDPJESpOqT57qT7BcQ6Q0U7Ydr1RnzOSJlMZMPuwY5InUVoryPKAsrFiuIVJaHbrEq8YgH4jyJ8yLehOpR/i45ZVBPgjOTJ5Bnkg5zv/qvU70GOSD4KzJsz5PpBxnYPc6zWOQD4BgkCdSmit59zjKM8gHwc9PmIik48zkvT7igEE+CM7snZ02RMrxM7VjkA+Cq0TDTJ5INc7knQuvYeT+hIMbBxEFwl2u8fa9GOSD4MjkBRdeiZTDhdewE+yTJ1IZF17DzlWuYSZPpBoeUBZygpk8kdK44zXsTG6GIlIZd7yGHTN5IqX52WDHIB8E1uSJlOZcbGW5JoyYyRNRCTP5MOIBZURKM1mTDzkf/6lGRPLhsQYhx6OGiRTnSvS8fSsG+SAwyBMpzWQmH3I8oIxIaX6WaRnkg8BMnkhpzhjPs2vCiEGeSGnijnfqj0E+CCb75IlU5jqF0uP3YpAPhHNlnZk8kWrc58mzXBM6wuQ1XolUJsAWynDjNV6JlOZqsPP4vRjkg+BaWmcmT6QawQPKQo7dNURK41HDYWcWyrcZ5ImUw7Nrws40yrfZQkmkHF4ZKuREoRzk2UJJpB6T5ZqQY7mGSHFceA23Qh7QtOJtlmuIlCMEAJ9CQFVBfmhoCMuXL8f7778PADh58iQ6OzvR0dGB/fv328+7cOECurq6sGTJEmzfvh2GYdzpJZUmzAIQiVt3gh0MEfnOFMKK8Z6rGOTPnDmDVatWoa+vDwCQyWSwbds2vPLKK3jjjTdw7tw5nDhxAgCwdetW7NixA8eOHYMQAgcOHPB08Pcs04AWZZAnIglOoTxw4AB27tyJVCoFADh79ixmz56N1tZWRKNRdHZ2ore3F1euXEEmk8G8efMAAF1dXejt7fV08PesggFEYsXbLNcQKccUgKb5k8tHKz1h165drvsDAwNIJpP2/VQqhf7+/nGPJ5NJ9Pf313Go4SFMA4haQZ6ZPJFqXKdQepzoVQzyY5mm6foNJISApml3fLxWzc0zav4zE4nHiz9aMtlUl9erB2tMUU0AiQbkAMyYHsfHJBqjTGT67GTHuaqOLPMU0XW79+K++xo9HVfNQb6lpQXpdNq+n06nkUqlxj1+/fp1u8RTixs3hmCaU//Ndn3oFqJ6FOn04JRfq15yOQPxeBRGLgvEGgAAg4OjyEo0Rlkkk01SfXYy41xVR6Z5yhvlNupbN0eQboxN+rV0XbtrclxzC+XcuXNx6dIlXL58GYVCAUePHkV7eztmzZqFRCKB06dPAwB6enrQ3t4+6YFPhRAC/SNpXBm6Gsj7VyIKBWhWdw0PKCNSjhB2B6Xr2GEv1JzJJxIJ7NmzBxs3bkQ2m8XChQuxdOlSAMC+ffvQ3d2NoaEhzJkzB+vWrav7gKsxmB8K5H2rZhoAu2uIlCUgSntlhOe9F1UH+ePHj9u3FyxYgCNHjox7TltbGw4ePFifkU1BxsjYtye7NuCpggGN3TVEynJl8kG3UN6LMoWsfTtbyAU4komxu4ZIbe7z5L19r3AGeaMc5DOFzF2eGRDTUZNnkCdSjhCOk008fq+QBvlyYJcxk0chz0yeSGGm4/AalmsmwV2uyd7lmf4TQhRPoSzV5L3+gIlIPszkp8iVyRsSZvKAY+GVmTyRaoTjgDJm8pMgcyZv/96OlBqbGOSJlGM6jhr2OpUPZZDPF/L2belq8tYHqrOFkkhVzky+Dhv87yqUQd4Q5S3D0gX5UpTX9EjpLjN5ItUUczt/UvlQBvm86czkJSvXWJm7rgOazmMNiBRU3KRZvM1MfhIMswBdK/5o8mXyJXqkGOQ9X1snItm4AjuDfO3yZt4O8s6sXgpWJq/pgKZBMJMnUo4zk2d3zSQYpgGt9D/Zgrxw1uQ1nTV5IgU5a/Is10yCYRrQNQ26piFfkOxi4tYHqpVq8uyuIVKOs7uGC6+TYGfymi5dJm9/oHqkuOWNmTyRUoQo/nu+XK7x9v1CGeTzpgFNk7Nc4+yu0ViuIVLO2JjOYw0mwcrkdU1D3pSsXGPRIsU2SgZ5IqVYC63WdS648DoJrky+IGcmr+k6AI01eSLF2A12pfsFj1deQxnkyzV5Ccs1cLRQMpMnUk45ky/eLxQY5GtmuGrykpVr7LNrii2UgkGeSClW4m6Va0yWa2qXt2vy8nXXCFcmH+GxBkSKGZ/JexsDQhnkDZlr8o4WSk3TixcQISJliDGZvMGafO3yjpp8TrJM3k7krUye5RoipdiZfOm+ySBfO0OUWygN2Wryzs1QOjN5ItWUa/LFr+yuqZEQwlWukS+Tdx41HOHCK5FixvbJsyZfI+uCIbrVQlnIy3mxbI2ZPJGKBDP5qTFKmXsxk9chIFAQEgVS17EGrMkTqaacdGqAxiBfM6OUGVs1eQCS9cpb/1RjJk+kImdM18AgX7O8K5PXXI/JwP4lrrNPnkhFzj55TdM83/Ea9fTVA2B102ilfwoBkKxX3nnUsA4hUymJiDwnxmTyXrdQhi7IW6UZa+W6+JhEQV5wxyuRysSYKF/wOAaELsg7M3mrXCNdGyVg73hlJk+kFiukaygmo17veA1dkJ8ok5dqQ5Rdj2MmT6QiZyavwfs++dAFeVdNviQncU2e3TVEanF1UGqAwYXX2hiOTN4K83LV5EtfS9013PFKpJby2TXFDZsGM/naWEFeh3PhVZ5yTfmoYfbJE6nI7z750AV5d02+1CcvU7mGO16JlOaqyWsa8gYz+Zq4++Tl2wxVrNdopYVXZvJEqnGeXaNpXHitWd5Vk5fwWAOBYnAHiiUbdtcQKcWdyXt/0ZApBfm1a9fi5s2biEaLL/ODH/wAw8PD2L17N7LZLB577DFs2bKlLgOtliEcffKahOUaiGJwBwCdffJEqnHvhZJ44VUIgb6+PvzlL3+xg3wmk8HSpUvx61//Gp/4xCewYcMGnDhxAgsXLqzbgCsxJsjkpdoMJUQ5k2efPJFyzLGZvKwtlBcvXgQAfP3rX8eHH36Ir371q/jkJz+J2bNno7W1FQDQ2dmJ3t5eX4N8fkyffEyPyrUZCij2x6O0IYqZPJFS3DV5iTP527dvY8GCBfj+97+PfD6PdevW4Zvf/CaSyaT9nFQqhf7+/ppet7l5xmSHBACIX9UR0XTE48UfLR6JIZIAksmmKb1uPcTjUeSFgB6NIplsws0Z05AzTSnGJiPOS/U4V9WRYZ5ujBQrC9FoBJGIAVN4O65JB/n58+dj/vz59v0nn3wSL730Ej772c/ajwkhXMcLVOPGjaEpncp2e2gEET2K4Vwxe49qUdweGkE6PTjp16yXXM4AICCEhnR6ENlRAxAmBgZu1zxPYZdMNknxmd0LOFfVkWWebt0cAQAYhgnTFMgbhSmNS9e1uybHkz5P/q9//Sveeecd+74QArNmzUI6nbYfS6fTSKVSk32LSTFMAzG9/Lsrpscka6FEsRbv/MpeeSJlmPaOV3/Ok590kB8cHMTevXuRzWYxNDSEQ4cO4dlnn8WlS5dw+fJlFAoFHD16FO3t7fUcb0WGaSCqOYJ8RLIg71p4LX1lkCdSxrgDymRtoVy0aBHOnDmDr3zlKzBNE6tXr8b8+fOxZ88ebNy4EdlsFgsXLsTSpUvrOd6K8mYB0XGZvGwLr8UMXrNaKU0TiAQ4HiLyzdgDyqQ+T37z5s3YvHmz67EFCxbgyJEjUxrUVBhmfky5JipXn7wQ0MZl8uywIVKFcJRrAM3zK0OF7hqvhjAmyOQlCvLOzVDOTJ6IlOD8r714rAGDfE2MseWaSEyyzVAYV5MXPL+GSBnCUa/RUAwJXmbzoQvy+THlmrhkmbyAYHcNkcLGHlAGeFuXD12QH5vJR/Uo8gWJFl6FcO94BXgSJZFC3BfyLkZ5LztsQhjkDUT1cquKdC2UADRm8kTKGnvREIBBviZ500BUj9n3ZSvXODN5MJMnUo5wbYYqPubl4mvogvy4zVDS9cmPr8kLdtcQKcPZJ2/l8szka1A81sBRrtGjMIWJgizZ8gTdNeyTJ1LH2IuGAFx4rYkxplwTixRvy1OycS68sk+eSDWCNfmpyYsxC6+6FeQlKdkI4Vh4ZSZPpJryAWXlq9exJl+DYrnGkcmXbuckOdpAwLnwypo8kWrcNfkiboaqkilMmMJ0ZfLxUs+8NOUaAcfCK7triFQzcU2eQb4q1mX+XJuhJK7Js0+eSD32jleUg7zBhdfq5CcI8jLW5K0MnjteidTjupA3WJOviZWtJ/S4/ZhdrpGkJg8IaNaagZ3JM8gTqcLZXWPV5VmTr1K2kANQbpt03pamXCMEECn9S0OzTqFkuYZIFXZ3jcYWyppZHTTxSDmTT0QSAIBsIRvImMYSQgDWLyEuvBIpx3QeNWwfUMaafFXyZjGTjztaKBOlgG9l+YETApqdyXPhlUg1wpwgk2dNvjoTZfINpUw+I0UmL4r/L2XyGjN5IuUUJjiGkuWaKuVK2Xo84szkS+UaQ4JM3qrFWeNjJk+kHNPVQskDympiXeYv7uiuiegRRPWoHDV5qxY3pibPy/8RqcN0RPnywitr8lUpl2tirscTkbhkQb5Uk9d5QBmRapxn14DnydfGLtc4MnmgWJeXpybvLNfwgDIi1Th74u3NUIJBvio5q7sm4g7yiUhCju6aMeUajZk8kXJcffLM5GtjlWtijmMNgFKQNyTI5EsLrMzkidTlzuSLuPBapZyZQ1yP2SvWFmlq8pi4Js8dr0TqMIWwM3jwylC1yRfy40o1ANAQlaNcI+7QXcM+eSJ1mCYQ0YvR3arJ8+yaKmULOdcFQywJWRZe7T5599k17JMnUodpCuilVJ41+RplCllMizaMe1yemvyYhVdNL37KzOSJlGEKAU0vl5QjuoZ8geWaqowaGTRMEOSnx6ZhxBiFGXTGbC28xhxj1CLM5IkUYpoCEce64bREFCMZ7653EaognzFGJ8zkp8emQ0BgxBgNYFQOVjCPJcqP6Tp3vBIppCAEdEcm39gQxUiWQb4qo0ZmwiA/IzYdADCUG/Z7SG5WJh91BHktwj55IoUIU8AR4zG9IYrhjHfXuwhdkJ+oXDMjXgry+aCDfKkmH3Vn8uyTJ1KHOS6Tj2GU5ZrqZIwMpkXunMkPBxzkhWkCml4+YhilXa/M5ImUUTDHBPlEFMMM8pXlCnkYonDXcs1gbsjvYbkJs3zxbovGTJ5IJUZBIBYpx4HpDVGMsFxT2WBuEADQFJ8x7nsfizdB13Tcynzo97DchFneAGXRI9zxSqSQvGEiFi3HgWmlhVfh0SFloQnyt0tZ+sfiTeO+F9EjmJn4OK5nbvo9LDdhlnc/WDSdffJECskbJqKuTD4GoyCQM7xJ9kIT5O+WyQNA87Rm3BgNOMibBWhjDk/TojGg4N0/1YhILkbBnck3JooxwateeU+C/O9//3t8+ctfRkdHB37729968RbjWFn6vzXMnPD7yWnNuDaSDnZDlFmAFom4HtISMyCyAXf9EJFvckbBHeQbrCDvTbJX9yDf39+P/fv347XXXsPhw4fx+uuv45///Ge932acq0P9mB5rvGMm/58fa8WoMYr+kbTnY7kTYRrA2Ey+YQZEZjCgERGR30YyBhobymdsWUHeqw6baOWn1ObkyZN49NFHcd999wEAlixZgt7eXnznO9+p6s87W4uqNZgbxpWRDzAv9Yj951OplOv1Hrm/DanGZvzP1f/Fk//VCX3sAqjHzNv9SD3QgtjHk66fMZr8DxiDA9AgXK2VNLm/C6riXFUn6HkSQiAW1fHvyen4v1KMmtnUgNTMaeN2wlar0p/RRJ2XdH/+859jZGQEW7ZsAQD87ne/w9mzZ/HCCy/U822IiKgKdU8dTdN0XbRDCDHuIh5EROSPugf5lpYWpNPlunc6nbZLJ0RE5K+6B/nPfe5zeOedd3Dz5k2Mjo7izTffRHt7e73fhoiIqlD3hdcHHngAW7Zswbp165DP5/Hkk0/iM5/5TL3fhoiIqlD3hVciIpIHe/aIiEKMQZ6IKMQY5ImIQoxBnogoxBjkPVLpkLY//elPWLlyJVasWIFnnnkGH330UQCjDF61h9m99dZb+OIXv+jjyORSaZ4uXryItWvXYsWKFfjGN76h7N8noPJcnT9/Hk888QRWrFiBDRs24Pbt2wGM0keC6u7atWti0aJF4tatW2J4eFh0dnaK9957z/7+4OCg+PznPy+uXbsmhBDipz/9qXjhhReCGm5gKs2TJZ1Oi6VLl4pFixYFMMrgVZon0zRFR0eHOHHihBBCiB//+Mdi7969QQ03UNX8nVq1apV46623hBBC7N69W/zkJz8JYqi+YSbvAechbY2NjfYhbZZ8Po+dO3figQceAAA8/PDDuHr1alDDDUylebJ0d3dXfcBdGFWap/Pnz6OxsdHedPj0009jzZo1QQ03UNX8nTJNE8PDxeO9R0dH0dAw/pKhYcIg74GBgQEkk0n7fiqVQn9/v31/5syZWLx4MQAgk8ng1VdfxZe+9CXfxxm0SvMEAL/61a/wqU99CnPnzvV7eNKoNE//+te/cP/992Pbtm14/PHHsXPnTjQ2NgYx1MBV83fq+eefR3d3N77whS/g5MmTeOqpp/wepq8Y5D1Q7SFtg4ODWL9+Pdra2vD444/7OUQpVJqnd999F2+++SaeeeaZIIYnjUrzZBgGTp06hVWrVuHQoUNobW3Fnj17ghhq4CrNVSaTwfbt2/HLX/4Sb7/9NlavXo3nnnsuiKH6hkHeA9Uc0jYwMIDVq1fj4Ycfxq5du/weohQqzVNvby/S6TSeeOIJrF+/3p4z1VSap2QyidmzZ+PTn/40AGD58uU4e/as7+OUQaW5evfdd5FIJOyjVr72ta/h1KlTvo/TTwzyHqh0SFuhUMDTTz+Nxx57DNu3b1f2KOZK87Rp0yYcO3YMPT09ePXVV5FKpfDaa68FOOJgVJqn+fPn4+bNm/j73/8OADh+/DjmzJkT1HADVWmuZs+ejWvXruHixYsAgD//+c/2L8ewqvsBZXTnQ9q+9a1vYdOmTbh27Rr+9re/oVAo4NixYwCARx55RLmMvtI8hf0/vmpVM08vv/wyuru7MTo6ipaWFuzduzfoYQeimrnavXs3Nm/eDCEEmpub8eKLLwY9bE/xgDIiohBjuYaIKMQY5ImIQoxBnogoxBjkiYhCjEGeiCjEGOSJiEKMQZ6IKMQY5ImIQuz/AWanbY7F4tYxAAAAAElFTkSuQmCC",
"image/svg+xml": "\n\n\n\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import seaborn as sns\n",
"import matplotlib.pyplot as plt\n",
"\n",
"parameters = np.zeros([3, 50])\n",
"\n",
"for i in range(50):\n",
" x_train, x_valid, y_train, y_valid = get_siso_data(n=3000,\n",
" colored_noise=True,\n",
" train_percentage=90)\n",
" \n",
" model.fit(x_train, y_train)\n",
" parameters[:, i] = list(model.theta)\n",
"\n",
"sns.set()\n",
"pal = sns.cubehelix_palette(3, rot=-.5, dark=.3)\n",
"\n",
"ax = sns.kdeplot(parameters.T[:, 0])\n",
"ax = sns.kdeplot(parameters.T[:, 1])\n",
"ax = sns.kdeplot(parameters.T[:, 2])\n",
"# plotting a vertical line where the real values must lie\n",
"ax = plt.axvline(x=0.1, c='k')\n",
"ax = plt.axvline(x=0.2, c='k')\n",
"ax = plt.axvline(x=0.9, c='k')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using the Extended Least Squares algorithm\n",
"\n",
"As shown in figure above, we have a problem to estimate the parameter for $y_{k-1}$. Now we will use the Extended Least Squares Algorithm.\n",
"\n",
"In SysIdentPy, just set *extended_least_squares* to *True* and the algorithm will be applied."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD7CAYAAACPDORaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAe9UlEQVR4nO3de2wU97UH8O/srm1qIA0lu6FFlntFE5ygBuhDgpSaug3mYS9uTNoGEPSKRDhKSwSpSKhxcZVcGkpQqaIECapKqG2ixCEFE5SYRwm+QkSlpSpcEHlUYEggmDUGbC/e18zv/rG749m1g8f2zs5vZ78fKcnueNk52NHZ4/P7zRlFCCFARESO5LI7ACIisg6TPBGRgzHJExE5GJM8EZGDMckTETkYkzwRkYMxyRMROZjH7gDSXb8ehKaNbOv+s88+jYICN/7nf17MUFQjJ2NMSePHj8G1az12h5GCMZknY1yMyZxM5AWXS8G4caM/9+vSJXlNEyNO8levXkVhoWfE75NJMsZkJGNcjMk8GeNiTIPLRl5gu4aIyMGY5ImIHIxJnojIwZjkiYgcjEmeiMjBmOSJiBxMui2URET5oqc3imKLz8EkT0Rkk0uBoOXnYLuGiMjBmOSJiByMSZ6IyMGY5ImIHIxJnojIwZjkiYgcjEmeiMgGQmRn7DGTPBGRDbI12Z5JnojIBtm6gQmTPBGRDbLUrWGSJyKyg8aePBGRc3HhlYjIwdiuISJyMLZriIgcjJU8EZGDcQslEZGDceGViMjBslTIM8kTEdmBlTwRkYNxdw0RkYOxXUNE5GBs1xARORi3UBIRORgvhiIicjCpFl5ffvllVFVVoaqqCps3bwYAHDt2DH6/H5WVldi6dav+2rNnz6K2thZz587F+vXrEYvFrImciCiHSVPJHzt2DEePHsXu3buxZ88enDlzBvv27UN9fT22bduGd955B6dPn0ZraysAYO3atdiwYQP2798PIQSampos/0sQEeUaaRZevV4v1q1bh8LCQhQUFGDSpEloa2tDaWkpSkpK4PF44Pf70dLSgkuXLiEUCmHatGkAgNraWrS0tFj+lyAiyjXStGvuuecePWm3tbXh3XffhaIo8Hq9+mt8Ph/a29tx9erVlONerxft7e0WhE1ElNuy1a7xmH3hxx9/jLq6OjzzzDNwu91oa2vTvyaEgKIo0DQNiqL0Oz4U48ePGdLrB1JYGP9reb1jR/xemSJjTEYyxsWYzJMxLsZ0e9duRfXHVsZlKsmfOHECTz31FOrr61FVVYXjx48jEAjoXw8EAvD5fJgwYULK8Y6ODvh8viEFdO1az4j3j0YiMRQWehAIdI/ofTJJxpiSvN6x0sXFmMyTMS7GNLjrnbf0xyOJy+VSblscD9qu+eyzz/Czn/0MW7ZsQVVVFQBg6tSpOH/+PC5cuABVVbFv3z6Ul5dj4sSJKCoqwokTJwAAzc3NKC8vH3bwREROla2e/KCV/B//+EeEw2Fs2rRJP/boo49i06ZNWLVqFcLhMGbPno158+YBALZs2YKGhgb09PRgypQpWL58uXXRExHlqGztrhk0yTc0NKChoWHAr+3du7ffsbKyMuzatWvkkRERORgHlBEROZg0WyiJiCjzpLkYioiIMk/TsnMeJnkiIhtIs/BKGaJGEevpsTsKIpKEdFe80shonZ9CjYUhYhEonkK7wyEim3Hh1WFELBJ/oEZv/0IiygtceHWaxAwfwSRPROA+eedSeRMVImIl70CJaZxM8kQE3sjbsYTI0uZYIpKaNLf/owxjkicicHeNAyV+oNn6+CYiqTHJOxUreSJCar1n5SIsk3y2MckTEVIreSuLeib5bGOSJyKkVfJgJe8c2Ro9R0RSM26hZCXvINxCSURAah+eSd5JmOSJCFx4dS5uoSQipC28WngeJvlsYyVPREir99iucRAuvBIRUit5Ky+MYpLPNlbyRAROoXSg+A+Uu2uICEjfQslK3jmY5IkI6RdDWYdJPluSP0UmeSICxxo4FxdeiQipt/9ju8ZJWMkTEXjFq3PxYigiAi+GcjAmeSLiWAPnYiVPRGC7xoF4+z8i6mPcg8FK3kGsvDkAEeUOXvHqVKzkiQicXeNc3EJJRJBsCmVPTw+qq6vx6aefAgB++ctforKyEjU1NaipqcHBgwcBAGfPnkVtbS3mzp2L9evXIxaLWRd5rmIlT0RIbddYWfoNmuRPnjyJxYsXo62tTT92+vRp/OUvf0FzczOam5sxZ84cAMDatWuxYcMG7N+/H0IINDU1WRZ47mKSJ6K0TGBnu6apqQmNjY3w+XwAgN7eXly+fBn19fXw+/146aWXoGkaLl26hFAohGnTpgEAamtr0dLSYlngOYuVPBEhe1soPYO9YOPGjSnPOzo6MGPGDDQ2NmLs2LGoq6vDrl27cM8998Dr9eqv83q9aG9vz3zEuY5JnoiQNrvGwvMMmuTTlZSU4JVXXtGfL1u2DHv27MGkSZOgKIp+XAiR8tys8ePHDPnPpCssjP+1vN6xI36vTBs9uhB3ShiXjN8rxmSejHExptsrKupLv3feWWxZbENO8h9++CHa2towd+5cAPFk7vF4MGHCBAQCAf11HR0deotnKK5d60kZpj8ckUgMhYUeBALdI3qfTDHeKCTY04uoJHEleb1jpfleJTEm82SMizENLtQb1R93dgYxaph7HV0u5bbF8ZDfVgiB3/zmN7h58yai0SjeeOMNzJkzBxMnTkRRURFOnDgBAGhubkZ5efnwonYa47ZJdmuICBK3a8rKyrBy5UosXrwYsVgMlZWVqK6uBgBs2bIFDQ0N6OnpwZQpU7B8+fKMB5yTDNcvZ+sqNyKSm/Hqdyvzgukkf/jwYf3x0qVLsXTp0n6vKSsrw65duzITmZOkXADFi6GISLKLoWiENLXvMSt5IkLaxVAca5DjUgdH2xcHEUkjW6mAST4LhLGS58orEYHz5J0lZXcNkzwRpf2Cb2HxxySfDUzyRJQmZd2VlXyOS7kFDHfXEBHbNQ7DhVciSsUbeTuJxuqdiFIZt01aWfoxyWeBEGzXENHnYyWf60R2Ll8motyhsSfvICnVO5M8EbEn7yzcQklEabL1Wz2TfDZwrAERpTGmghHeQuO2mOSzge0aIkojslT8MclnA3fXEFEaIQAk7pBqZVZgks8GtmuIKI0QAvpdsNmuyW2CC69ElCZ1dg3bNbmNPXkiSiMEoCjxWp5XvOY6tmuIKA0HlDkJ2zVElEYT+ror2zU5j2MNiKifvizPdk2uS+nJcwslEWWvkvdY9s42uhXrxY1o2O4w+qTe54uIKJ7YFQWAsDQvODLJfxa8Ak0IqJoKt8ttdzgcNUxE/QhDJc+xBkOUHOHZE71lcyQJKb+KsZQnovQ9GFx4HZaIGrE7hDjuriGiNAICSnLhlZX88ES1qN0hxDHJE1GabC3VMclng/7TTCyyEFHei8+uUfTHVnF0ko+osiT5RCWvgJU8EQFInULJds0wSVnJM8kTEVKnULKSHyZ5kny8klcUhVsoiQhAWiVv4Xkcl+RjWkx/HDU8tlPfPnmFYw2ICEBid03yMSt584zbJqPS9OQTP0AF4MIrEQHJtKAYHlvDcUk+bEjyMaHaGImBoZJnT56IgERPnguvQ2dM8qomWZJXuIWSiOK0lH3ybNeYFkmp5OXoyXN3DRENpK8nb905TCX5np4eVFdX49NPPwUAHDt2DH6/H5WVldi6dav+urNnz6K2thZz587F+vXrEYtlP8mG1b7pk/JV8mCSJyIAiRlbyuCvG6lBk/zJkyexePFitLW1AQBCoRDq6+uxbds2vPPOOzh9+jRaW1sBAGvXrsWGDRuwf/9+CCHQ1NRkafADiRi2TcrWk1fALZREFBfP8RJc8drU1ITGxkb4fD4AwKlTp1BaWoqSkhJ4PB74/X60tLTg0qVLCIVCmDZtGgCgtrYWLS0tlgX+eYxXucpTySd312ThY5uIcoIwVPJW/oI/6Dz5jRs3pjy/evUqvF6v/tzn86G9vb3fca/Xi/b29iEHNH78mCH/GaNRwb758YWjXPB6x47o/TLhRnFB4pGCAo8cMaVjTObIGBMgZ1yM6fYURYHLFc/yY8YUWRbbkG8aomla/MrNhPg2IOVzjw/VtWs90EYwQb/jRpf+uDvYi0Cge9jvlSnhnlD8gQJEI1EpYjLyescyJhNkjAmQMy7GNDhV1SASua6rOzTs2Fwu5bbF8ZB310yYMAGBQEB/HggE4PP5+h3v6OjQWzzZFE3srlGgQJWl/23cQsmFVyJCoievSHgx1NSpU3H+/HlcuHABqqpi3759KC8vx8SJE1FUVIQTJ04AAJqbm1FeXp7xgAcTSYwycLtcUCUZa2DcQmnlflgiyh3GXCDVjbyLioqwadMmrFq1CuFwGLNnz8a8efMAAFu2bEFDQwN6enowZcoULF++POMBD0av5BWXZLtrsrDCQkQ5w7iD0sqsYDrJHz58WH88c+ZM7N27t99rysrKsGvXrsxENkwRLQqXosAFRZ7dNRrHGhBRKs6TH6aIGoECFxRFkaeST35Oc0AZESUY58lbmeWdl+QTlbwCSNOTF5oGJGLixVBEBCRm1yQWXkewoXBQjkvyUTUKRUlU8rK0a4zVO9s1RAQAhnnyVnJcko9oUbigJLZQSpLkteTCK6dQElGcSKn92K4xLaJGoCiKnJU8B5QRUYJmnCdv4Xmcl+S1aKKOV6DKMmpY07IzU5SIcgYr+WGKqlG4pOzJJ9o1TPJEeS+Z1KW84lV2qtASRbNsPfkE7q4hynvpOZ2V/BCoWiy+XVGBPJW80PStUqzjiUiv5JPPLTyX85K80Aw9eTmSvDDOk2e7hijvGe8ImvLcAo5L8jERi6d4RUFMkouhUlo0bNcQ5b30pM52zRDE59rLOGo4G7+YEVEu6NeuYSVvnipUvZKXZaxB6rg5JnmifNeXBZT4pjsLz+W4JB8fShbvyUszoCyxTsAtlEQEGNfp9COWnctxSV7VVP2KV01o0GRo2eh3hgLYriEiY61n9YXwjkrymtAgEkN/EjMf5ejLp17aZl8cRCSF9EKeSd6k5E1Ckj35+DEJ+vL6PnmFu2uICOm/0XN3jUl6Dz45u914zEYikdjj3RpW8kT5LmV+vKJw4dWs5MVPCmCo5O1P8vp9vjiFkojSxNMCK3lT1MSMGEVR9J68FKMNRN89XgXbNUR5T9MHlMWfsydvUnK0cEpPXoZxw9wnT0RG/a54te5UzkryyUoeMlbySuJjm0meKN8Z2zPxrMB2jSl61a4YevISLLymfExrbNcQ5bt+V7yykjdHFcZKPnlMhiRvnF3DJE+U7/oldSZ5cwbeJy9BUtV78hxrQETGO0PFn2ts15iT3BNv3F0jw8KrMI414O4aoryXPtaAlbxJfZW8YayBFJV8sl3DSp6I0hdaFe6TN6uv/67Iu/AKay98ICL5idQczytezVIN7Zr0Y7Yy3ONVf05EeSvt7n+s5M2KaQOMNZAhoRrv8Wp8TkR5qd8+efbkzembXWNYeJXmYijj/QEk+OAhIvvwitfh0TTD7ppERpWmXZNceI0fsDEYIrKbljZQnle8mhQzLrxCnoVXkdwnr18PxUqeiBJ4xat5KaOGpdxCqR+wKxIikkBfIa+kPLeCs5K8oV0DqbZQ9o0ajj9nkifKZ+kLr1YWfp6R/OFly5ahs7MTHk/8bZ577jkEg0G88MILCIfDmD9/PtasWZORQM2Qd3ZN+j55LaWuJ6L8kpoRFEvrvmEneSEE2tra8N577+lJPhQKYd68efjzn/+ML3/5y6irq0Nraytmz56dsYBvx5jQk1soNVnaNYbfLljJE+W39BRgZUYYdpI/d+4cAGDFihW4ceMGfvzjH+Pee+9FaWkpSkpKAAB+vx8tLS3ZS/Ja/9k1MQlm1/RbaOUWSqK8ZhxQpiiSXgzV1dWFmTNn4pVXXsHOnTvx+uuv4/Lly/B6vfprfD4f2tvbMxKoGcZ98sn/SnExFBL3eNWfspInymf9KnkZ2zXTp0/H9OnT9eePPPIIXnrpJXzzm9/UjwkhUkYMmDF+/JjhhoSiKx4oUFBYGP9ruV1uFI1yw+sdO+z3zIQgBNzuvnbN+C8Vw3OHvTGls/t7NBDGZJ6McTGmz3czHC9IPZ54nV1Y6LEstmEn+X/+85+IRqOYOXMmgHhCnzhxIgKBgP6aQCAAn883pPe9dq0Hmja8j7Xunl64FRd6IzEUFnrgUlzoDvYiEOge1vtliqZpUDWhf7OvdXTDFS6yNSYjr3es7d+jdIzJPBnjYky319kZBADEYvFOQygcHXZsLpdy2+J42O2a7u5ubN68GeFwGD09Pdi9ezeefvppnD9/HhcuXICqqti3bx/Ky8uHe4ohU4UKl8utP3crbjnaNVraPnkZYiIi2yWvkZSyXVNRUYGTJ0/ihz/8ITRNw5IlSzB9+nRs2rQJq1atQjgcxuzZszFv3rxMxntbqlDhVoxJ3iXJFkrukyeiPiljDSxeeB3RPvnVq1dj9erVKcdmzpyJvXv3jiio4VI1FR4ltZLXpBhQFl941ZcnmOSJ8pueAnjF65CoQoPb2K5xSdKuEVoiwSd/oBLERES2SU3q1l4a6bAkr8Kt9P2V5GvXJJ+zkifKZ8mpk8mevCbjPnkZqVp6T94tyTz5xD55/QOblTxRPku//Z+VnJXkhZa6u0aCdk18QSXtZl8yjFogItukL7SykjdpoHaN7WMN0m/9Fz9oSyhEJAdhqPsUwNKU4MAkn767xuaqOeU3CW6hJCLeyHvYNE1LTfIut/0Lr8YkrwxwjIjyTv92jXXnclSSV4UKt8vYrpEhyQ/QrmElT5TXUtZdFWXYo1zMcF6S77e7Rr52jZW/mhGR/LjwOkxqv3aNBPvk00capBwjorykN+XjV8IzyZskdbsGYE+eiAD09eCTKcHKhoOjknxsoHaNzUleH2FgGGvAnjxRvjPcyFtRWMmbpaVd8eqSqiffd3NxJnmi/Ja+ziq48GpOvF0jW0/e2K5JVvJs1xDltbQuLit5k1ShSdeugZq44pZbKIkoQb+Rd+Jf3CdvUnxAWerCq+1XvOoD0pLz5gAOKCPKb8Zx8grAffJmDdSusXt2jUgmeWMlb+XHNhFJT6S1cdmuMWnAi6Hs7n9ryQ+ZvlHDgpU8UV5L31nNSt6kAXvyds+T1yt5/V/syRPlOU3vyVt/FbxjkrwQAjEtBo+r77a1bpcbAgKajdW8SOnJ6wdtiYWI5GC8kbfChVdzYoldNAXGJJ+o6m1t2WiG3TUKK3ki6qvz4tsxOKDMlJgWBYDUSj6x08bWls1AlbzdLSQispWW1pTnwqsJ0UTFXJDWrgEAzc698ol98oqxkmeSJ8prycpd4RZK82KJJO9xFejH7G7XCE1F76FXEs+UvkUW1eZbEhKRrbS0+/+xJ29CVE/yxt01iXaNTZW8CHYCaryNlNKT15jkifKZcQqlAu6TNyWmt2sGqORtao+IaLjvicvVl+RZyRPlNS1t1jAHlJkQTSy8DtSTt21+TaS377Hi1pO8YCVPlNeMlTtvGmJSLFGtD7i7xq6efNSY5A2za1jJE+W1ZOXOm4YMQXTALZTxSj5mW7smBAAYNfsxAIkdNi43e/JEeU4zrruykjcnJuEWSpFo17gn3t930OXh7hqiPJc61kBJOZZpjkny4Vh8kbPIXaQfG+UeBQDojYVsiQmR+HmVglF9x9weaNcuWjqrgojkZlx4VdKPZZhjknxIjSf5UZ6+JD+6oBgAEIzesiUmvSdvTPLhINTLZ6F+9oEtMRGR/bS0m4YA1g0pc1ySN1byySR/K2ZXkg8BnkIohr37SdrNdhsiIiIZGKt2qxdfnZPk9XZNoX6sOJHk//fT9+1pj0R6U1s1BqK3K8vBEJEsYqqAxx1Pv0pia7VqUZZ3TJIPq2EUuQvhMtz+r8DlQdm4e3A5eAUXuj/Jekwi0gsUfiHl2OilW+NfC3ZmPR4ikkM0pqHAE89VrkSSD0eZ5G+rJxpEsae43/H/nrIYAPDx9XPZDgkicgtKYWpMrtHj4LqrFFoPkzxRvoqqhiTviif5UMSaXXeWJPm3334bCxYsQGVlJV599VUrTtHPzXAX7iz6Yr/jYwvH4K4vjEdb18WsxGEkwkEoo8b0O+4a/SWI4PWsx0NEcojGVBS448ndlcjC4ag1W70znuTb29uxdetWvPbaa9izZw/eeOMN/Oc//8n0afq5Hr6BLxaNHfBrX72jBG1dNrRrQt1QRvWPSRk9DhrbNUR5KxxRUVgQ35Cht2si1iR5z+AvGZpjx45hxowZuPPOOwEAc+fORUtLC37+85+b+vPJX12G4kboBgCBe780CS6XAp/Ph4ICt/5eU7xluNB9ERe6L+K/vlg65PcfDhENwe12o8Bb2i8mz/gSiEv/B0WLQvEUDv5mWTCc77vVGJN5MsbFmD6f2+1C6YSx0Hw+qAJwjfsCVG148Q32ZxSR4W0n27dvx61bt7BmzRoAwJtvvolTp07h+eefz+RpiIjIhIy3azRN07cEAfEN/sbnRESUPRlP8hMmTEAgENCfBwIB+Hy+TJ+GiIhMyHiSf/DBB/H++++js7MTvb29OHDgAMrLyzN9GiIiMiHjC69333031qxZg+XLlyMajeKRRx7BAw88kOnTEBGRCRlfeCUiInk45opXIiLqj0meiMjBmOSJiByMSZ6IyMFyPsmbHYb2zDPP4K9//asUMR06dAg1NTVYuHAhnnzySdy8edP2mA4ePAi/34+qqiqsW7cOkUjE8pjMxJV05MgRfP/735cippdffhkVFRWoqalBTU1N1obwDRbXuXPnsGzZMixcuBCPPfaY7f9fnT17Vv8e1dTU4Lvf/S6qq6ttjQkAzpw5g0WLFmHhwoWoq6tDV1d27u0wWFytra3w+/3w+/34xS9+gWAwmJkTixx25coVUVFRIa5fvy6CwaDw+/3i448/7veauro68cADD4i33nrL9pi6u7vFd77zHXHlyhUhhBC///3vxfPPP29rTMFgUMyaNUsEAgEhhBCrV68Wr7/+uqUxmYkrKRAIiHnz5omKigopYqqrqxP/+te/LI9lKHFpmiYqKytFa2urEEKIF198UWzevNnWmIxu3bolqqqqxD/+8Q/bY1q8eLE4cuSIEEKIF154Qfzud7+zNCYzcd28eVPMmDFDP7Zjx46M5YWcruSNw9CKi4v1YWhGb7/9Nn7wgx9g/vz5UsQUjUbR2NiIu+++GwAwefJkfPbZZ7bGVFxcjMOHD+Ouu+5Cb28vrl27hjvuuMPSmMzEldTQ0GB6wF02Yjp9+jS2b98Ov9+P5557DuFw2Pa4zpw5g+LiYv3CwyeeeAJLly61NSaj7du349vf/ja+9a1v2R6Tpml6ldzb24tRowa+e1s242pra8NXvvIVfO1rXwMAVFRU4NChQxk5d04n+atXr8Lr9erPfT4f2ttT7536+OOP40c/+pE0MY0bNw5z5swBAIRCIezYsQMPPfSQrTEBQEFBAVpbW/G9730P169fx6xZsyyNyWxcf/rTn3D//fdj6tSplsdjJqZgMIj77rsPa9euxe7du9HV1YVt27bZHtfFixdx1113ob6+Hg8//DAaGxtRXNz/JjrZjCmpu7sbTU1NWfmgNhPTunXr0NDQgFmzZuHYsWN49NFHbY/rq1/9Kq5cuYIPPvgAAPDuu++io6MjI+fO6SQv4zA0szF1d3dj5cqVKCsrw8MPPyxFTLNnz8bf//53VFRU4Ne//rWlMZmJ66OPPsKBAwfw5JNPWh6L2ZhGjx6NP/zhD5g0aRI8Hg9WrFiB1tZW2+OKxWI4fvw4Fi9ejN27d6OkpASbNm2yNaakvXv34qGHHsL48eMtjcdMTKFQCOvXr8fOnTtx9OhRLFmyBM8++6ztcd1xxx347W9/i1/96ldYtGhRYjR5QUbOndNJXsZhaGZiunr1KpYsWYLJkydj48aNtsd048YNHD16VH/u9/vx4Ycf2h5XS0sLAoEAFi1ahJUrV+rfNztjunz5Mnbt2qU/F0LA48n4dJAhx+X1elFaWoqvf/3rAIDq6mqcOnXK1piSDh06hAULFlgai9mYPvroIxQVFemjVn7yk5/g+PHjtselqiomTJiAN998E2+99Rbuu+8+lJSUZOTcOZ3kZRyGNlhMqqriiSeewPz587F+/fqs/OYxWExCCKxduxaXL18GEE+u3/jGN2yP66mnnsL+/fvR3NyMHTt2wOfz4bXXXrM1plGjRuHFF1/EJ598AiEEXn31Vb39Zmdc06dPR2dnp/7r/uHDhzFlyhRbYwLi/2+dOXMG06dPtzQWszGVlpbiypUrOHcufs/nv/3tb/oHo51xKYqCFStWoL29HUII7Ny5M3MfjBlZvrXR3r17RVVVlaisrBQ7duwQQgjx+OOPi1OnTqW87tlnn83K7prBYjpw4ICYPHmyWLhwof5PfX29rTEJIcTBgwdFdXW18Pv9Ys2aNaKrq8vymMzElfTJJ59kZXeNmZhaWlr0r69bt06Ew2Ep4vr3v/8tFi1aJBYsWCBWrFghOjo6bI+po6NDPPjgg5bHMZSYjhw5Ivx+v6iurhY//elPxcWLF6WI67333hPV1dWisrJSNDY2ikgkkpHzckAZEZGD5XS7hoiIbo9JnojIwZjkiYgcjEmeiMjBmOSJiByMSZ6IyMGY5ImIHIxJnojIwf4f+Z5hJ8BrCT0AAAAASUVORK5CYII=",
"image/svg+xml": "\n\n\n\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"parameters = np.zeros([3, 50])\n",
"\n",
"for i in range(50):\n",
" x_train, x_valid, y_train, y_valid = get_siso_data(n=3000,\n",
" colored_noise=True,\n",
" train_percentage=90)\n",
"\n",
" model = PolynomialNarmax(non_degree=2,\n",
" n_terms=3,\n",
" extended_least_squares=True,\n",
" ylag=2, xlag=2,\n",
" estimator='least_squares',\n",
" )\n",
" \n",
" model.fit(x_train, y_train)\n",
" parameters[:, i] = list(model.theta)\n",
"\n",
"sns.set()\n",
"pal = sns.cubehelix_palette(3, rot=-.5, dark=.3)\n",
"\n",
"ax = sns.kdeplot(parameters.T[:, 0])\n",
"ax = sns.kdeplot(parameters.T[:, 1])\n",
"ax = sns.kdeplot(parameters.T[:, 2])\n",
"# plotting a vertical line where the real values must lie\n",
"ax = plt.axvline(x=0.1, c='k')\n",
"ax = plt.axvline(x=0.2, c='k')\n",
"ax = plt.axvline(x=0.9, c='k')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great! Now we have an unbiased estimation of the parameters! \n",
"\n",
"## Note\n",
"\n",
"Note: The Extended Least Squares is an iterative algorithm. In SysIdentpy we fixed 30 iterations because it is knwon from literature that the algorithm converges quickly (about 10 or 20 iterations). Also, we create a specific set of noise regressors for estimation (order 2 and nonlinearity also equal 2). It is suposed to work in most cases, however we will add the possibility to choose the order and nonlinearity degree for noise regressors in future releases."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.7.6 64-bit ('sys_tests': conda)",
"metadata": {
"interpreter": {
"hash": "75eea971140cd20bf61a98694439aebfe611e1f965354b67c59fd475d54e3010"
}
},
"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.7.6-final"
}
},
"nbformat": 4,
"nbformat_minor": 4
}