Skip to content

Instantly share code, notes, and snippets.

@aemarkov
Created May 18, 2019 13:41
Show Gist options
  • Select an option

  • Save aemarkov/4bf31096c16b72ffca9bb78a67ea6707 to your computer and use it in GitHub Desktop.

Select an option

Save aemarkov/4bf31096c16b72ffca9bb78a67ea6707 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
" # RANSAC\n",
" ## Линейная регрессия и МНК для долбоебов\n",
" \n",
"Дан набор точек $(x_i, \\tilde y_i$). Надо аппроксимировать прямую $y = ax+b$.\n",
"\n",
"Минимизируем ошибку:\n",
"$$\n",
"F = \\sum_{i=1}^N{(ax_i+b - \\tilde y_i)^2}\n",
"$$\n",
"\n",
"Ошибка есть просто функция двух переменных $a, b$. Чтобы минимизировать найдем частные производные и приравняем к нулю. Ведь так обычно ищут экстремумы?\n",
"\n",
"Крч, немного вычислений и вот:\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"\\frac{\\partial F}{\\partial a} &= a\\sum{x_i^2} + b\\sum{x_i} - \\sum{x_i\\tilde y_i} = 0\\\\\n",
"\\frac{\\partial F}{\\partial b} &= a\\sum{x_i} + Nb - \\sum{\\tilde y_i} = 0\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Ну а это простая система двух линейных уравнений, которую можно решить на бумажке, либо, если вы тупой, как и я, нумпаем.\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"&a\\sum{x_i^2} + b\\sum{x_i} = \\sum{x_i\\tilde y_i} \\\\\n",
"&a\\sum{x_i} + bN = \\sum{\\tilde y_i}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"\\sum{x_i^2} &\\sum{x_i} \\\\\n",
"\\sum{x_i} &N\n",
"\\end{pmatrix}\n",
"\\cdot\n",
"\\begin{pmatrix}\n",
"a\\\\\n",
"b\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"\\sum{x_i\\tilde y_i} \\\\\n",
"\\sum{\\tilde y_i}\n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import sys"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def make_data(x, coefs, n_data, n_noize):\n",
" a,b = coefs\n",
" xl = np.random.uniform(x[0], x[1], size=(n_data))\n",
" yl = a*xl+b\n",
" \n",
" _min = min(x[0], a*x[0]+b)\n",
" _max = max(x[-1], a*x[-1]+b)\n",
" \n",
" xn = np.random.uniform(x[0], x[-1], size=(n_noize))\n",
" yn = np.random.uniform(a*x[0]+b, a*x[-1]+b, size=(n_noize))\n",
" \n",
" data = np.concatenate((np.vstack((xl, yl)).T, np.vstack((xn, yn)).T))\n",
" np.random.shuffle(data)\n",
" return data\n",
"\n",
"\n",
"def regression_line(data):\n",
" x, y = data[:,0], data[:,1]\n",
" N = len(x)\n",
" A = np.array([\n",
" [np.sum(np.square(x)), np.sum(x)],\n",
" [np.sum(x), N]\n",
" ])\n",
" B = np.array([np.sum(x*y), np.sum(y)]).reshape(2,1)\n",
" a, b = np.linalg.solve(A, B)[:,0]\n",
" return a, b\n",
"\n",
"def draw_line(x, coefs):\n",
" a, b = coefs\n",
" y = a*x+b\n",
" plt.plot(x,y, '-b')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Линейная регрессия хорошо работает без выбросов"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4.999999999999998 2.000000000000012\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFYVJREFUeJzt3Xtw3HW9xvHnYy8WhBFKYq2kNYxAe1CHywkMF6dcSoWm5eLIcERhekKxDNgDeFlb/YMWZ9SycYqco0etQBuPUC4VoTQo1ApSPRw0hXqhl7FAC4WUbooV6SC29HP+2MVms7vJJtn9Xd+vGSbZX9PuZ0d4fPrd7/e35u4CAMTfu8IeAABQGwQ6ACQEgQ4ACUGgA0BCEOgAkBAEOgAkBIEOAAlBoANAQhDoAJAQI4N8soaGBm9ubg7yKQEg9tatW9fj7o0D/Vyggd7c3Kyurq4gnxIAYs/MtlXzcyy5AEBCEOgAkBAEOgAkBIEOAAlBoANAQhDoAJAQBDoAJASBDgB1tGrVbn3iE6uVy/XU/bkIdACog337pA9/WLrggsP0wAPTtGTJj+v+nAQ6ANTYQw9Jo0ZJGzbkH1999XJdffXldX/eQI/+A0CS7d4tHX74gccf+9g/9MQTo2V2WSDPT0MHgBo4++ziMJdO1oUX3iqz4GaoqqGb2VZJf5P0tqR97t5iZmMl3SOpWdJWSZe6+1/qMyYARNOLL0of/GDxtVyuR0uXXqq2trZAZxlMQz/b3U9w95bC4/mS1rj7MZLWFB4DQGqMG1cc5o88IrlLDQ0NymQyamhoCHSe4ayhXyTprML3HZIelzRvmPMAQOT9/vfSCScUX3MPZ5beqm3oLulRM1tnZnMK18a5e3fh+x2SxpX7jWY2x8y6zKwrl8sNc1wACJdZcZg/80w0wlyqPtA/5u4nSZou6XNmNqX3L7q7Kx/6Jdx9ibu3uHtLY+OAH7gBAJG0erWK3uBsbMwHed+mHqaqllzc/eXC151m9lNJp0h61czGu3u3mY2XtLOOcwJAaPruVNm2TZo4MZxZ+jNgQzez95jZoe98L+njkv4kaaWkWYUfmyXpwXoNCQBhWLasOMzPPDPfyqMY5lJ1DX2cpJ9a/lWNlHSXu//czH4n6V4zmy1pm6RL6zcmAARn/35pxIjia6+91nefefQMGOju/ryk48tc3yVpaj2GAoCwzJ4t3XHHgcdXXSX98IfhzTMYHP0HAElvvCEdemjxtb//XXr3u8OZZyg4+g8g9d773r5hfr+y2fZYhblEQweQYtu3SxMmFF979dUedXQ8F/ix/VqgoQNIJbPiMJ8/P7+D5X3vC+fYfi0Q6ABS5f77S/eVZ7Pt+uY3w5mnllhyAZAafYN8ypRfaebM38ZyeaUcAh1A4n3609Ly5cXXWltnqKOjQw0NZ4YzVB0Q6AASrW8rP/XU/9V55z2quXM7YrlO3h8CHUAilfukoPxdEU8v/JM8vCkKIFHcS8P8kkt+rlyuJ5yBAkRDB5AYlVv5+UGPEgoaOoDYe+ON0jC/9to7U9HKe6OhA4i1yq38M0GPEjoaOoBY2rChNMx37ozOx8GFgYYOIHYqt/J0o6EDiI1Fi0rDfPr0CwjzAho6gFgo18pbW2do8eLFwQ8TUTR0AJF2/PHlb6aVy/Wos7NTkyZNCmewCKKhA4isymvlmaBHiQUaOoDIMavcylEZDR1ApPQN8qYm6aWXJFr5wAh0AJHAVsThY8kFQKj27y8N87PP/j/CfAgIdAChMZNGjCi+ls226957jw5noJgj0AEErru7tJXPnn2f3BXbD2iOAtbQAQSq3Fp5NtuemM/1DBOBDiAQDz8szZhRfK27W3r/+yV2sNQGgQ6g7tjBEgzW0AHUzbXXlob5okXthHmd0NAB1AVr5cEj0AHU1KhR0r59xdcONHLWyuup6iUXMxthZs+Y2arC46PM7Ckz22Jm95jZ6PqNCSAOzPoLc9TbYNbQr5e0sdfjmyXd4u5HS/qLpNm1HAxAfJS7mZY7YR60qgLdzJokzZB0W+GxSTpH0orCj3RIurgeAwKItr5B3txMkIel2jX0b0v6sqRDC4+PkLTb3d/5y9V2SUeW+41mNkfSHEmaOHHi0CcFEClsRYyeARu6mc2UtNPd1w3lCdx9ibu3uHtLY2PjUP4IABGyd29pmB999J2EeQRUs+RyhqQLzWyrpLuVX2q5VdJhZvZOw2+S9HJdJgQQGWbS6D7bH1pbZ2jVqpZwBkKRAQPd3b/i7k3u3izpU5J+6e6fkfSYpEsKPzZL0oN1mxJAqDZvLm3lc+bcw+d6RsxwTorOk/QFM9ui/Jr67bUZCUCUmEmTJxdfy2bb9YMf/Bt3RYyYQR0scvfHJT1e+P55SafUfiQAUXDrrdINNxRf27Rpl1auvIPTnhHFSVEAJSrvYDlCmQynPaOKm3MB+KemptIwv/nmduVyPeEMhEGhoQOQ1F8rp5HHBQ0dSLlyx/azWVp5HNHQgRSjlScLgQ6kEMf2k4klFyBlCPPkoqEDKUGQJx8NHUi4PXtKw/zUU9cT5glEoAMJZiYdckjxtWy2XQ891BTOQKgrllyABHrySen004uv/epX0pQpEjtYkotABxKm3Fp5LtfDjbRSgCUXICGuu640zHfvzr/xSZinAw0dSAB2sECioQOxxrF99EZDB2KKY/voi4YOxAytHJXQ0IEYoZWjPzR0IAZo5agGDR2IuMr7ymnlKEZDByKqXCvP5XrYV46KCHQgYnp6SoP8pJOe5bQnBsSSCxAhld/0/HDQoyCGaOhABNx1V2mYjxnTymlPDAoNHQhZuVY+atRoLVv2P8EPg1gj0IGQHHectHFj8bUXXujRffctVVvbK6yXY9AIdCAE5Vr5ggUL1dy8UJkM2xExNAQ6EKBKe8qXLl2qtra5wQ+EROFNUSAglVp5Q0ODMpkMSywYNho6UGe0cgSFQAfqqPK+8gbWylFzAy65mNkYM/utmf3ezJ41s5sK148ys6fMbIuZ3WNmo+s/LhAP3EwLYaimob8l6Rx3f8PMRkn6tZn9TNIXJN3i7neb2fclzZb0vTrOCkSeu/SuMjWJm2khCAM2dM97o/BwVOEfl3SOpBWF6x2SLq7LhEBMmJWGuTsf0ozgVLXLxcxGmNl6STslrZb0nKTd7r6v8CPbJR1Z4ffOMbMuM+vK5XK1mBmIlOeeK11emTDhZyyvIHBVBbq7v+3uJ0hqknSKpMnVPoG7L3H3FndvaWxsHOKYQDSZSUcfXXxtwYKFevrpk2nlCNyg9qG7+25Jj0k6TdJhZvbOGnyTpJdrPBsQWV/7WmkrP+20LymX69HChQsJc4RiwDdFzaxR0l53321mB0maJulm5YP9Ekl3S5ol6cF6DgpERbmtiNlsu9ra5hPkCFU1DX28pMfM7A+SfidptbuvkjRP0hfMbIukIyTdXr8xgfCV24r4jW8slrs46YlIGLChu/sfJJ1Y5vrzyq+nA4lXuZW3BT8MUAH3cgH6UemAEK0cUUSgAxXQyhE33MsF6KNckLe2zlBnZ6ckTnsiumjoQC+Vwnzx4sXBDwMMEg0dUH93RZSkziBHAYaMho5U27t3oDAH4oOGjtQiyJE0NHSkzpNPlob5NdcQ5og/GjpShVaOJKOhIxUuv7w0zNevJ8yRLDR0JB6tHGlBQ0dilTu2P336BYQ5EouGjkQq18rPPXeabrnlO8EPAwSEQEei9L+8sjrIUYDAseSCxGCtHGlHQ0fsEeRAHg0dsUaYAwfQ0BFLBDlQioaOWHnzzdIwHz2aMAckAh0xYiYdfHDxtWy2XW+9Fc48QNSw5ILIW7tWmjKl+Nqxx35Hl13Wo7a2ueEMBUQQgY5I6+9zPfmAZqAYSy6IpCuvLA3zD33oTG3atFmZTIYwB8qgoSNyyrXyM888SytW/IQgB/pBoCMyygX51KnTdN55H1db2wrCHBgAgY5IKBfmra0ztHjxdzRp0qTgBwJiiDV0hKrcLW6z2Xa5S52dnYQ5MAgEOkLT3w4WAIPHkgsCVy7IJ0/+F61du1YNDZngBwISgoaOQFVaK3/ggQd40xMYJgIdgWCtHKi/AQPdzCaY2WNmtsHMnjWz6wvXx5rZajP7c+Hr4fUfF3Hz+uulQT5jxluslQN1UM0a+j5JX3T3p83sUEnrzGy1pH+XtMbdF5nZfEnzJc2r36iIm8q3uH23JNbKgVobsKG7e7e7P134/m+SNko6UtJFkjoKP9Yh6eJ6DYl4+cUvSsP8qqvuUy7XE85AQEoMapeLmTVLOlHSU5LGuXt34Zd2SBpX4ffMkTRHkiZOnDjUORETlbYiZjI0cqDeqn5T1MwOkfQTSTe4++u9f83dXVLZjxhw9yXu3uLuLY2NjcMaFtE1a1ZpmG/Zsou1ciBAVTV0MxulfJjf6e73Fy6/ambj3b3bzMZL2lmvIRFtldfKj6CZAwGqZpeLSbpd0kZ3X9zrl1ZKmlX4fpakB2s/HqKs3FZEdz4ODghLNQ39DElXSPqjma0vXPuqpEWS7jWz2ZK2Sbq0PiMiiviQZiB6Bgx0d/+1pDL/+UqSptZ2HERd5bsiLpbE4SAgTNzLBVWr3Mo7gx4FQBkEOgbE8goQD9zLBRW5l4b5yJF7OCAERBSBjrLMpHf1+bcjm21Xd/eb3BURiCiWXFBk1y6pb15Pm/Yb3XXXJO5VDkQcgY5/qrxWfkbQowAYApZcoIcfLg3zdet44xOIGxp6yrGDBUgOGnpKXXVVaZh//eu3EOZAjNHQU6jSLW65KyIQbwR6ivS/vMIOFiDuWHJJCdbKgeSjoSccQQ6kBw09wQhzIF1o6AlEkAPpRENPkHI30zr5ZMIcSAsaekKUa+W5XA830gJShIYeczt2lIb5smX5Vk6YA+lCoMeYmTR+fPG1bLZds2aV/3kAyUagx1BnZ7mbab3GaU8g5VhDj5nKO1jG6qSTOO0JpBkNPSY++9nSMN+7lx0sAA6goccA+8oBVIOGHmEHH1wa5tlsO2EOoCwaekRxi1sAg0VDjxizyq08k8mwtxxARQR6hPQN8jFjnFYOoGosuURA5Tc9TXzwBIBq0dBD9PbbpWF+ww3sYAEwNDT0kLAVEUCt0dAD9tJLpWH+6KOEOYDhG7Chm9kdkmZK2unuHylcGyvpHknNkrZKutTd/1K/MZOBVg6gnqpp6Msknd/n2nxJa9z9GElrCo9RwSOPlIb5jTd+lzAHUFMDNnR3f8LMmvtcvkjSWYXvOyQ9LmleDedKDA4IAQjKUNfQx7l7d+H7HZLGVfpBM5tjZl1m1pXL5Yb4dPGTyZSG+fTpM7Rp02YOCAGoi2HvcnF3N7OKiwfuvkTSEklqaWlJxSJDuVa+YMFCLVzYGfwwAFJjqA39VTMbL0mFrztrN1J8HX98aZjncj3KZts1d+7ccIYCkBpDbegrJc2StKjw9cGaTRRTlXewNCiT4bQngPobsKGb2XJJT0qaZGbbzWy28kE+zcz+LOncwuNUKnczrdbWGcrlesIZCEBqVbPL5bIKvzS1xrPETt8gv+KKv+ujH/0vtbV18KYngMBx9H8IKi+vjBE30wIQFo7+D8K+faVhvnw5pz0BRAMNvUoc2wcQdTT0AezYURrm8+b9kDc9AUQODb0flVv5Z4MeBQAGREMvY+3a0jDfs4clFgDRRkPvg7VyAHFFQy+47bbSMN+/nzAHEB80dJVv5blcj8w4HAQgPlLd0K+55s2Kx/Y56QkgblLb0PNBftA/H8+c+ZamTPlPju0DiK3UBfrJJ0tdXcXXDjRyju0DiK9ULbmYFYf5j36Uf9OTRg4gCVLR0BsbpZ4+BzvZvQIgaRLd0Pfuzbfy3mF+3XU/IswBJFJiG3q5rYjZbLva2tqCHwYAApC4QP/rX6XDDiu+tmuXNHasxJueAJIsUYHOsX0AaZaINfRXXikN80WLvsUtbgGkSuwbet8gP/bYF7R581GSvhTKPAAQltg29K1bS8M8m23Xb35zaCjzAEDYYtnQDz9c2r37wONPfvIRrVhxnnjTE0Caxaqhd3fnW3nvMM9m2/X97/9reEMBQETEJtA//3npAx/o/XiZ3KVMJsPRfQBQTJZcrrzyTS1dmr8z4k037dFBB/03B4QAoI9YBLpZp6RG3XTTet144/VirRwASsUi0G+++SxNnryUVg4A/YhFoDc0NCiToZUDQH9i86YoAKB/BDoAJASBDgAJMaxAN7PzzWyzmW0xs/m1GgoAMHhDDnQzGyHpu5KmSzpO0mVmdlytBgMADM5wGvopkra4+/Pu/g9Jd0u6qDZjAQAGaziBfqSkl3o93l64VsTM5phZl5l15XK5YTwdAKA/dd+H7u5LJC2RJDPLmdm2Qf4RDZLS9kkVaXzNEq87TdL4mqWhv+4PVvNDwwn0lyVN6PW4qXCtIndvHOyTmFmXu7cM9vfFWRpfs8TrDnuOIKXxNUv1f93DWXL5naRjzOwoMxst6VOSVtZmLADAYA25obv7PjObK+kRSSMk3eHuz9ZsMgDAoAxrDd3dH5b0cI1mqWRJnf/8KErja5Z43WmSxtcs1fl1m7vX888HAASEo/8AkBCRDfQ03lbAzCaY2WNmtsHMnjWz68OeKShmNsLMnjGzVWHPEhQzO8zMVpjZJjPbaGanhT1TEMzs84V/v/9kZsvNbEzYM9Wamd1hZjvN7E+9ro01s9Vm9ufC18Nr/byRDPQU31Zgn6Qvuvtxkk6V9LmUvG5Jul7SxrCHCNitkn7u7pMlHa8UvH4zO1LSdZJa3P0jym+o+FS4U9XFMknn97k2X9Iadz9G0prC45qKZKArpbcVcPdud3+68P3flP8PvOT0bdKYWZOkGZJuC3uWoJjZeyVNkXS7JLn7P9x9d7hTBWakpIPMbKSkgyW9EvI8NefuT0h6rc/liyR1FL7vkHRxrZ83qoFe1W0FkszMmiWdKOmpcCcJxLclfVnS/rAHCdBRknKSlhaWmm4zs/eEPVS9ufvLkr4l6UVJ3ZL+6u6PhjtVYMa5e3fh+x2SxtX6CaIa6KlmZodI+omkG9z99bDnqSczmylpp7uvC3uWgI2UdJKk77n7iZL2qA5/BY+awrrxRcr/H9oHJL3HzC4Pd6rgeX57Yc23GEY10Ad9W4GkMLNRyof5ne5+f9jzBOAMSRea2Vbll9bOMbMfhztSILZL2u7u7/wNbIXyAZ9050p6wd1z7r5X0v2STg95pqC8ambjJanwdWetnyCqgZ7K2wqYmSm/prrR3ReHPU8Q3P0r7t7k7s3K/+/8S3dPfGNz9x2SXjKzSYVLUyVtCHGkoLwo6VQzO7jw7/tUpeDN4IKVkmYVvp8l6cFaP0Hd77Y4FCm+rcAZkq6Q9EczW1+49tXCiVwkz39IurNQWp6X1BbyPHXn7k+Z2QpJTyu/q+sZJfDUqJktl3SWpAYz2y5pgaRFku41s9mStkm6tObPy0lRAEiGqC65AAAGiUAHgIQg0AEgIQh0AEgIAh0AEoJAB4CEINABICEIdABIiP8HmwLzAlF6N/0AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"data = make_data((0, 10), (5, 2), 100, 0)\n",
"plt.plot(data[:,0], data[:,1], 'ok', markersize=1)\n",
"\n",
"a, b = regression_line(data)\n",
"draw_line(data[:,0], (a, b))\n",
"print(a,b)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Но с выбросами полный пиздец"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.409041313229257 14.872409701991671\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt0nHW97/H3txSEhnKRVuwBKnhZiagLe9qtQN3SDa1gi2333sr2BjWiKBJQg4Qe3a62btkbJsegEtSNQijLC6BoQRoRiiCyFKSlXtA2G1CQegqd0QIShLb0e/74TUiapJ2ZZOa5fl5rzUpmnuk835k++c7v+f4uj7k7IiKSfhPiDkBEROpDCV1EJCOU0EVEMkIJXUQkI5TQRUQyQgldRCQjlNBFRDJCCV1EJCOU0EVEMmJilDubMmWKH3nkkVHuUkQk9datW1dy96mVnhdpQj/yyCNZu3ZtlLsUEUk9M3u0muep5CIikhFK6CIiGaGELiKSEUroIiIZoYQuIpIRSugiIhmhhC4ikhFK6CIiu1Eqlejs7KRUKsUdSlWU0DMkLQdfWuIU6enpoaOjg56enjG/RpTHe6QzRaWxBg4+gAsuuCDmaHYvLXGKtLa27vJzLKI83pXQM6QeB18U0hKnyJQpU8adhKM83s3dG76TAbNmzXKt5RKPUqlET08Pra2tTJkyJe5wRKQGZrbO3WdVel5VNXQze8TMfmtmvzKzteXHXmpmt5nZg+WfB4836KRIS423ljjrUQsUkd1LQt6opVP0n9z9jUO+JZYCt7v7a4Dby/czIS3Jr5Y4W1tbKRQKKnOINEgi8oa7V7wBjwBThj3WB0wr/z4N6Kv0OjNnzvQ0KBaLXigUvFgsJnpfUcaZVfoMpV5GO5bqdXwBa72aXF3Vk+CPwP3AOuCs8mNPDtluQ+/v7paWhB6lQqHggBcKhbhDySV9/tJI9Tq+qk3o1Y5yeYu7/9nMXgbcZmYbh7Xy3cxG7V01s7OAswCmT59e5e7yQyM+4qXPP5+iGiQQ+fFVTdYfegOWA58iwyWXRtHpfWX6jCQKaTszo14tdDNrAia4+9/Kv78N+BxwE7AEuLj888a6f9tkjCbUVKbPSKKQ1TOzakouhwI/MLOB53/b3W8xs/uA683sTOBR4LTGhZkNWT2I6kmfkUShHhOGkkgTi0REhknaRLy6TiwSkSAJk0ek8RIxpnwMtJaLSA1U48+HtJb+lNBFapDWP3SpTVpr7Cq5iNRg4A89CXVVqY8sldGU0EUaKEvJIqvSWi8fTapLLknriRYZTjX35MtSGS3VCV1/LJJ0WUoWWTBaIzCt9fLRpDqh649Fki5LySILst4ITHVC1x+LiNQi641AdYqOgzq8RNIl66OUlNDHIUu94yJZ0dfXxwknnMCcOXPo6+uLO5xIpbrkEresn76JpFF7ezt33XXXi7+vXr065oiio4Q+DqrhiyRPV1cXzzzzDGZGV1dX3OFESgldJIU0B2P3mpub+elPfxp3GLFQDb3B1HEqjaD+GxmNWugNlvVxrxIP9d/IaJTQG0x/eNIIee2/Ualpz1RyabCsj3sViVK9Sk1ZLYWqhS4iqVGvM96slkKV0EUkkUqlEt3d3Tz77LNMmjSJtra2upWasloKVUIXkUTq6elhxYoVL95vamqqW2s6q30QSugikkitra309/e/2ELPWmu6EczdI9vZrFmzfO3atZHtT0SSTyNXKjOzde4+q9LzNMpFRGI10EG5ZMmSzI06iZoSekJkdRiVpEccx2CpVKK/v5+5c+fS29urma/jpISeEJrKLXGL4xgc6PicPXs2hUJBdfJxUqdoQtQ6jEp1R6m3OIbyDd2njuPxU6doSnV2dtLR0UGhUMjk8CsRGVRtp6ha6CmV1YkRIjJ2SugpldWJESIyduoUlVTSqKBkK5VKdHR0MG/evNxd1zNOVSd0M9vLzNab2c3l+0eZ2b1m9pCZXWdm+zQuTMm74Qlco4KSraenh87OTtasWUN7e3vc4eRGLSWXjwMbgAPK9y8BLnX3a83sa8CZwFfrHJ8IMHJ1PPUhJFtrayvFYpH169fn7rqecapqlIuZHQ6sBC4C2oF3AEXg5e6+w8yOA5a7+8l7eh2NcpGx0jBNybN6T/3/ItAB7CzfPwR40t13lO9vAg6rOcqYqQ6bHnm4UIiORxmvigndzE4Ftrj7urHswMzOMrO1Zra2WCyO5SUaRnVYSZK0HY99fX3MmzePjo4OfQklRDU19NnAQjObD+xLqKF/CTjIzCaWW+mHA38e7R+7+xXAFRBKLnWJuk5Uh5UkSdvx2N7ezpo1a1izZg1Tp07VMNoEqGmmqJnNAT7l7qea2XeBG4Z0iv7G3b+yp3+vGrpIdvT19dHW1saMGTPo6OjIdDksblHMFL0QuNbMPg+sB64cx2uJSMo0Nzdz2223xR1Gov3qV/ChD8ERR8ANN8CEBs/8qenl3f1Odz+1/Psf3P1N7v5qd3+Xuz/fmBBFJC7qqK3dD38IhxwCZjBjBqxbB6tWwXPPNX7fmikqIqMqlUosWbIkVR21cXjhBfjyl0MCN4OFC+Gvfw3bXvIS+M53YOdOmDSp8bEooYvICAPJvLe3l/nz56emozYqzzwD7e0hgU+cCB//+OC25ma4+25wD63yd787PC8KSugiNchLCaKnp+fFZL5y5Up1eAKPPQaLF4fkPHkyXHrp4LaTT4aHHw5JfONGmD07nhiV0EVqUO1Y8bQn/tbWVgqFQiKTeZSf7bp18MY3hiQ+fTrceOPgtrPPhq1bQxK/5RZ45Sujj28Ed4/sNnPmTE+7YrHohULBi8Vi3KFkQto+z2rjLRQKDnihUIgosvxo5Ge7c6f797/vfsAB7iFV73q75BL3bduijw9Y61XkWCX0GukPtb6y+nmm5YsqLXEOVe+Yt293/8IXRk/gTU3u3/1uSPRxxeeuhN4wafwDSDJ9nvHK6hdqJU895d7WNnoSf93r3O+5J+4Id1VtQtc1RUVyLE+rWD76KJxzDqxePXLbO94Bl10Gr3hF9HFVQ9cUFZGKsn4pw3vvhTPPhN/9buS2c8+Fz38eDjhg5La00igXSYy0jwyR+LmH5H344WFkyrHH7prMu7pg+/bwvC9/OVvJHNRClwQZflUiGZ9SqUR3dzcAbW1tmS2pvPAC3HNPmF6/ahU89NDgtgMPhKuvDuPH80AJPaHyVNsckLblY5Oup6eHFStWANDU1JSpL8nnnoM1a0IC/+EPYcsW2HtvOPFEOP98OPXU0ErPGyX0hMpjazXr9dyoDDQGFi5cSH9/P5CNL8mtW0OH5qpVYSJPf3+YsblgQWiBn3JKaJHnmRJ6Qqm1KmM1tDGwfPnyeIMZp8ceC7MzV62CO+8M5ZVp0+D000MSnzMnLIAlgYYtJlyl0kseSzMy0tDjAEjtMeEODzwwWA+///7w+GtfGxL44sUwa1bj1xVPGg1bzIhKpZc8lmZkpOHHQZqOhRdegJ//fDCJ/+EPgyNULrkEFi0KKxhKZUroCVep9KLSjED6joO//x1uu22wU7NUgn32gblz4cILw0SfadPijjJ9VHIRSaE0ltr+8pfBTs0f/xiefTZ0Yg7t1Jw8Oe4ok0klF9lFGhOA7F5aSm2PPBI6NW+8Ee66K5RXDjsMPvCBkMRPOCG0zKU+lNBzIi0JQEbX19dHe3s7XV1dNDc3J7bE4g6/+c1gPfxXvwqPv+51sHRpSOIzZ0Z3BZ+8UULPiaQmAKlOe3s7vb29AKxevTpRY/Z37AiXXBsYXvjIIyFhH388dHaGTs3XvCbuKPNBCT0nkpQApLLhLfKuri6AF3/G7dln4dZbBzs1//rXMB583jz4938PMzUPPTTuKPNHCV1yJQ19CaVSicWLF7Nx40YgtMibm5tZPdq6r5HGBTffHJL4rbeGkSoHHRSS9+LF4bqa++8fa4i5p4QuuZLkvoSBL5v+/n42btxIS0tL7C3yP/5xsB5+992wcycccQR86EOhlPLWt4Y1VCQZlNAlV6LqS6j1TKCvr+/FVvmyZcsoFAqxnEW4h47MgST+m9+Ex9/wBvjMZ0JLfMYMdWomVjWXNarXLQuXoJPsaOTl72q9tNv8+fMd8JaWFi8Wi5Femm/7dvef/MT9vPPcp08Pl2GbMMH9H/8xXGvzoYcaHoJUQJWXoFMLXXKrkeWXas4Ehrbih3Z6Tpkyhc7OzoaWhvr7w+SeVatCXXzrVth3X3jb22D58lAXnzq17ruVBlNCl9xqZPmlmlFFw79QhnZ6NiK2YjGMSFm1Kky7f+45OPjgMM1+8eKQzJua6rY7iYGm/otEKOpVER9+eLAe/vOfh07NV7xicOXCt7wFJqpZl3ia+i+SQI1eFdE9LDk7kMQfeCA8fswx8NnPhiR+zDHq1MwqJXSRCDWilLJ9O/z0pyGB33gjbNoU1gt/61vh0kvD8MKjjqrb7iTBKi4Tb2b7mtkvzezXZvY7M1tRfvwoM7vXzB4ys+vMTEvsiAzT19fHggUL6OvrAwZr6+Mtsfztb/C978H73w8ve1mYoXnVVeHiD1dfDU88AXfcAZ/4RLqSealUorOzk1KpFHcoqVRNC/154ER3f8bM9gbuNrMfAe3Ape5+rZl9DTgT+GoDYxVJneFrsIzHE08MdmquWQPPPw+HHDJYD583DyZNqkfU8enu7mbFihX09/en/vJ5cajYQi8Pg3ymfHfv8s2BE4HvlR9fCSxuSIQiKTK8Rd7V1cX8+fOrmvE5Wuv0wQfDAlezZ4cLPnz4w/D738PHPhbKLI8/Dj09oayS9mQu41dVDd3M9gLWAa8GLgceBp509x3lp2wCDtvNvz0LOAtg+vTp441XJLHGuwZL6DC9kMceezmTJ5/OqlUheUOYnbl8eWiJv+EN2e3UbGtro6mpSauCjlFVCd3dXwDeaGYHAT8AWqrdgbtfAVwBYdjiWIIUSbLxrsGybVu4ov2GDW0ceOBHueyyyey1V7j4w0c/CgsXhqGGeaBVQcenpmtnu/uTwB3AccBBZjbwhXA48Oc6xyYCJL+jbOhQxEKhwM9+9jOaK1zV+Omn4frr4b3vDTMyTz4ZrrtuP046aTLXXANbtsDtt8O558aXzJP+uctIFVvoZjYV2O7uT5rZfsA84BJCYn8ncC2wBLixkYFKfiV1hcSBlvnChQsBKk4Q2rwZbropDC28/fbQMp86Fd71rlADnzsX9tsvqugrS+rnLrtXTcllGrCyXEefAFzv7jeb2e+Ba83s88B64MoGxik5ltSrLVWT8Pr6Bif53HNPeOxVrwot78WL4bjjYK+9ooq4Nkn93GX3NPVfMieqi1iMtp+dO+G++waTeLl/lFmzQgJftChcXzOrnZrSGJr6L7nViFLBaMl7oAPv+efhlltCAr/pplBamTgR5syBtrbQqXnEEXUJQ2SPlNAlcxpRKhj+JfHUU/CjH4Uk3tsbZm42NcHb3w4nnfQ3nniih3POeW9iL3Mn2aSELpnTiKFvra2tPPXU/sDpnHxymFa/fXuYdv/ud4dyyoknhjXFOzu/xvLlHUya9Lw6EyVSSugiQ5RKJQqFAuvXr+eyy7pxby7Xw6fwy1+eDcCrXx3WSFm8GN785pGdmupMlLgooYsMceWVPXR23g0sYtasJvr7w+NvehNcdFFI4q997Z47NTU5RuKihC6599xz8JOfDIxMOR+4ALMdvPGNz/O+94VOzcNGXdhCJFmU0CWXHnkEzjkndGgOmDwZ5s+fwKJF8Pa3T+Sgg/TnIemiI1Zy4xe/gA9+cHBs+FBnnvldLr/8XbzkJdHHJVIvNa3lIpIm7nDddWE6vRkcf/yuyfw//uMZNm8uUSh0cvHF/6RkLqmnhC6Zsn07XHxxSOATJoQhhc89F7YdfHCokxcKnYDxkpd8lZe/vD5XEJLRaYGvaKnkIqm3dSssXQpXXDFy24wZ8I1vwPTpYabn7NmtzJ6tYYVR0QJf0VJCl0Tb3bosDz8MZ58Nt9028t/867/Cl76068iUzs5dE4uSSzQ0Jj9amUvoUS3MJNEY2sI77rgL+OAHw2XZhjv//HBFn/33H/11lFjioTH50cpcQtcpXna4w+TJH2XChHY6OkauMdvdHa7oU83ys0oskgeZS+h5b4ml/Qxl2zYoFOCznx14ZPKL26ZODRdEXrAgltBEEi9zCT3vLbE0nqH85S/Q0QFXXTVy2z/8A3z963DMMdHHJZI2mUvoeZeWM5T/+Z9QLrnjjpHb/u3f4NJLYdq06ONKkrSfbUn0NA49YwbOUJKYAO68E446KowRb27eNZlfeCH094e6+bXXKpnD4NlWT09PrHFoLHl6qIUuDeMOK1fC7k4W/vu/4UMfChOAZKSknG2lsYyXV0roUlfPPQf/9V/wuc+N3DZtWujUPPnk8e0jL6WIpPQHJeWLRSpTQpdxKxbhU5+Ca64Zue3440NL/PWvr9/+1GKMVlK+WKQyJXQZkw0b4CMfgZ/9bOS2970PvvAFOPTQxuxbLUaR0al6KVUplUqcddb1HHbYC5jB0Ufvmsw/8xl49tlQN//mN+uTzHfXGZfkjl+ROKmFLru1c2cYG/7hDwNMAU7bZfuVV8IHPtC4Tk2VVkRqo4Quu/j738O1My+6aOS2gw56miuvdP7lXw6MJBaVVkRqo5KL8MQToe5tBpMm7ZrMTzgh1MvdYevWAxqazIeXWFRaEamNEnpOPfBAGIFiBi9/OXz724PbliyBLVtCEr/zTmhpiSampEykEUkrlVxy5Mc/DjXvxx8fuW358jBbc999o45qkEosIuOjhJ5hO3eGq/icffbIbWZw9dVw+unh9yTQeGeR8VHJJSWqXU/j2WdDS9ssrBM+NJm/8pVh/RT3kOzPOCM5yVxExq9iQjezI8zsDjP7vZn9zsw+Xn78pWZ2m5k9WP55cOPDza891Zc3b4bTTgvJuakprCc+4KSTwsqG7uGybXPmRBfzcFrkKVr6vPOnmpLLDuB8d7/fzCYD68zsNuADwO3ufrGZLQWWAhc2LtR8G15f/vWvw/jw++4b+dwzz4RLLoFDDokywso0rjxa+rzzp2JCd/fNwOby738zsw3AYcAiYE75aSuBO1FCb5gpU6Zw9NEX8NrXwmgNrosuCuup7LNP9LFVS52e0dLnHY1ELRbn7lXfgCOBPwEHAE8OedyG3t/dbebMmS7V27HD/bLL3EPBZNfb3nu7f+tb7jt3xh2lSL4VCgUHvFAoNGwfwFqvIkdX3SlqZvsDNwCfcPenh30pOOC7+XdnmdlaM1tbLBZr/8bJmWeeCS1tM5g4Ec49d3BbczPcfXdI6du2wXvfm8xOzb6+PhYsWEBfX1/coYg0XGtrK4VCIRFnQhZycYUnme0N3Az82N27yo/1AXPcfbOZTQPudPfmPb3OrFmzfO3atXUIO1s2bYLzzoMf/GDktpNPhssvh1e9Kvq4xmrBggX09vYyf/58Vq9eHXc4IqlnZuvcfVal51WsoZuZAVcCGwaSedlNwBLg4vLPG8cYay7df3+4Ws/69SO3ffSj8J//CQendNxQV1fXLj9FJBoVW+hm9hbgZ8BvgZ3lhz8N3AtcD0wHHgVOc/e/7um18t5Cv/HGMFPzySdHbrvkEvjkJ2HvvSMPS0QSrm4tdHe/m9DpOZqTag0sT3bsgO7ukKiH22+/MFPzXe9KZh28klKpRHd3NwBtbW3x9+5LzRI1OkPqQlP/6+zpp+Gzn4Uvf3nktte9Dr7xDTj22OjjqqdSqcSSJUvo7e0FoKmpSeOcU6ie49T15ZAMmvpfB3/6EyxcGFraBx64azJfsAD++McwMuWBB6JP5rXMFqz03IHRK4VCgd7eXubOncuyZctobW3VrMQUqufoDK2UmRDVjG2s1y1L49B/+Uv3N7xh9DHi557r/uSTcUcY1DJGdk/PLRaL3tLS4oDPnTvXC4WCF4vFMe1HsqdYLI44JqR+qHIcukouVXKHG24InZr9/SO3f+ELYejhxIR9orXMFtzTc3t6eti4cSMtLS10d3fT3Nxc9b+V7NNKmclQ1Tj0eknbKJcdO+CLX4TRjtMDDoCeHvjnf05np2atVCMViU/dRrnkzVNPwac/DV/5yshtxxwTOjVnVfxYs0ctMJHkU0IndFqecw786Ecjty1aFDo5p0+PPi4RkVrkdpTLL34BRx8dyiWvfOWuyfyTnwzDD91h1SolcxFJh9wkdHe49tpwzUyzcIHkDRsGt3/pS6Fm7g5dXTB5cnyxRk1DDmWsdOwkS6YT+vbtcPHFIYFPmADveQ88/3zYdsghYSr+wGDD884Ll2zLI40hlrHSsZMsmauhb90KS5eGiyMPN3MmfP3rMGNG9HElyfARKxpyGI8sjBzSsZMw1QxWr9etUROLHnrIfd680Sf5vPOd7ps2NWS3qVQsFn3u3LkO+LJly+IOJ9c0GUuqRdYnFt19N3zwg/DggyO3XXABLFsWLpgsu+rp6WHNmjVxhyGodSv1l5qE7g7f+laYqfnCCyO3X345fOQj+a2DV6u1tZX+8lTXtra2mKPJN43tl3pLRafoypWhU/P00weT+aGHwurVg8WVj31MyXw0t9xyCy972cu45ZZbgJBEli9fzvLly1NbtxWR0aUioQ9485vh178OCfzxx2H+/LgjSqahQ8nOOOMMisUiZ5xxRtxhRUZD6SSvUlFyWbIk3KQ63d3drFixgv7+fq655hrOOOMMrrnmmrjDikw91/kWSZNUJPS4ZGFY2SmnnMKWLVviDiNS6myUvFJC34O0tPSGf/G0tbXR1NSU24SmzkbJKyX0PUh6S28gkff397NixQogfPHUktCycBYiIkGqOkWjNpAYk5rohp5BjPVSYpq6vWfqYJU0UQs9hQZa1QsXLgQYV+s66WchcUtL2U0ElNBTZXcllj09t1Kyr1SeyUpJZqzvIw1feFn5P5LxU8klBQZO+7u7u6susdSrlJKVksxY30fSy26Qnf+jJElrqU0t9IQrlUosWbKE3t5eli1b9mIir5Rg6tWyTEMLtRpZeR+jyfJ7i0taS226SHTCdXZ20tHRwfz581m5cmWiW4oiWZG0MpYuEp1iQw+moa2vJBxYInmQ1rkMqqEn0NCaaBpquI2U1lpmNZLy3pISh4yfWugJsbtWed6ltZZZydC+EYj3vWX1M84jJfSIVKrJDf+j0h9WkNUvt56eHnp7e5k/f37s7y2rn3EeKaFHpFIrSH9Uo0trLbOSJPWNZPUzzqOKo1zM7CrgVGCLu7++/NhLgeuAI4FHgNPcfWulneV5lMvQFjqQqB70LEjKqISkxCHZUu0ol2o6Ra8GThn22FLgdnd/DXB7+b7swdDOzYEJQt3d3XGHlRlJmVyTlDgknyqWXNz9LjM7ctjDi4A55d9XAncCF9Yxrsy55557aG1tHfMfulp+e5aUklVS4pCccveKN0Jp5YEh958c8rsNvb+n28yZMz2vWlpaHPCWlhYvFoteKBS8WCxW/e+XLVvmgC9btmzEtoHX27hxY82vKyLJB6z1KnLsuMehl3e220K8mZ1lZmvNbG2xWBzv7lKrp6eHlpaWhowtHzjNb29v1+m+SIJEPcZ/rKNcnjCzae6+2cymAbu9xpm7XwFcAaFTdIz7S71jjz2WDRs2jPnf7+kqRAOPLVy4kDlz5uh0XyQhoh7jX9VaLuUa+s0+OMqlE/iLu19sZkuBl7p7R6XXycMoF9W6RWRAvfJB3Ua5mNl3gF8AzWa2yczOBC4G5pnZg8Dc8n1BoxxEZFDUS3dUM8rlPbvZdFKdY0mtvr4+2tvb6erq0igHEXlR1Gfsmik6DgP/Wbfeeitr1qwBYPXq1Zp1JyJA9DV0JfRxGPjPuuCCC9hnn33o6uqKOyQRSZCoz9iV0MegnhdpFpHsinqdHCX0MdByoyKSRLrARZWGThBobW2teJFmkTzTRTPioYReQV9fHwsWLKBQKOgqQjmjpDR2Gr4bD5VcKmhvb6e3t5dt27apVZ4zKq2NnYbvxkMJvYKBkStdXV00NzfHHI1ESUlp7HTRjHhUNfW/XvIw9T8PtLyBSLTqeYGLzFOttDaqj4okk0ouqFZaK5UiJC46O9wzJXSUoGqV1fqokkXyqfG1Zyq5EP2KaEmgMtPIz0ClpOTTHJA9Uws9p9TSGfkZ6Ewt+bJ6dlgvSug5peQ18jNQspC007BFEZGE07BFEZGcyWRCV4efiORRJhO6RivUTl+CIumXyU5RdfjVTqNeRtK4dEmbTCZ0jVaonb4ER4r7S05fKFKrTCZ0qZ2+BEeK+0su7i8USZ/U1tDzUPPNw3tMsrhnEGtWpNQqtQk9Dx2feXiPsntxf6FI+qS25BL36XAU8vAeRaR+UjFTVJ1DIpJnmZopqtKDiEhlqSi5qPQgIlJZKhK6htSJiFSWipKLiIhUpoQuIpIR40roZnaKmfWZ2UNmtrReQYmISO3GnNDNbC/gcuDtwNHAe8zs6HoFJiIitRlPC/1NwEPu/gd33wZcCyyqT1giIlKr8ST0w4DHhtzfVH5sF2Z2lpmtNbO1xWJxHLsTEZE9aXinqLtf4e6z3H3W1KlTG707EZHcGs849D8DRwy5f3j5sd1at25dycweHcO+pgB5XHIwj+9b7zkf8vieYezv+xXVPGnMa7mY2UTgf4CTCIn8PuC97v67Mb3gnve1tpp1DLImj+9b7zkf8vieofHve8wtdHffYWZtwI+BvYCrGpHMRUSkOuOa+u/uvUBvnWIREZFxSMtM0SviDiAmeXzfes/5kMf3DA1+35Guhy4iIo2Tlha6iIhUkPiEnrf1YszsCDO7w8x+b2a/M7OPxx1TVMxsLzNbb2Y3xx1LVMzsIDP7npltNLMNZnZc3DE1mpl9snxsP2Bm3zGzfeOOqRHM7Coz22JmDwx57KVmdpuZPVj+eXA995nohJ7T9WJ2AOe7+9HAscA5OXjPAz4ObIg7iIh9CbjF3VuAY8j4+zezw4DzgFnu/nrCCLl3xxtVw1wNnDLssaXA7e7+GuD28v26SXRCJ4frxbj7Zne/v/z73wh/4COWVMgaMzscWAB8I+5YomJmBwJvBa4EcPdt7v5kvFGa2Hu3AAAB4UlEQVRFYiKwX3kuyyTg/8UcT0O4+13AX4c9vAhYWf59JbC4nvtMekKvar2YrDKzI4EZwL3xRhKJLwIdwM64A4nQUUAR6CmXmr5hZk1xB9VI7v5n4P8CfwI2A0+5+63xRhWpQ919c/n3x4FD6/niSU/ouWVm+wM3AJ9w96fjjqeRzOxUYIu7r4s7lohNBP438FV3nwH0U+dT8KQp14wXEb7M/hfQZGbvjzeqeHgYYljXYYZJT+g1rxeTBWa2NyGZf8vdvx93PBGYDSw0s0cIZbUTzeyb8YYUiU3AJncfOAP7HiHBZ9lc4I/uXnT37cD3geNjjilKT5jZNIDyzy31fPGkJ/T7gNeY2VFmtg+h8+SmmGNqKDMzQk11g7t3xR1PFNz9/7j74e5+JOH/+CfunvlWm7s/DjxmZs3lh04Cfh9jSFH4E3CsmU0qH+snkfGO4GFuApaUf18C3FjPFx/X1P9Gy+l6MbOB04Hfmtmvyo99urzMgmTPucC3yg2WPwCtMcfTUO5+r5l9D7ifMKJrPRmdNWpm3wHmAFPMbBOwDLgYuN7MzgQeBU6r6z41U1REJBuSXnIREZEqKaGLiGSEErqISEYooYuIZIQSuohIRiihi4hkhBK6iEhGKKGLiGTE/wcJkba7B1SHcQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"data = make_data((0, 10), (5, 2), 100, 100)\n",
"plt.plot(data[:,0], data[:,1], 'ok', markersize=1)\n",
"\n",
"a, b = regression_line(data)\n",
"draw_line(data[:,0], (a, b))\n",
"print(a,b)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## RANSAC\n",
"Поэтому запилили алгоритм RANSAC (Random Sample Consensus), который так прост, что даже программист на HTML сможет его реализовать:\n",
"\n",
"1. Взять $M$ случайных точек из выборки\n",
"2. Сделать линейную регрессию\n",
"3. Посчитать, сколько точек из всей выборки отклоняется не более, чем на $p$ от нашей модели\n",
"4. Если это число больше текущего максимального, обновить максимальное (ну и сохранить хорошие коэффициенты)\n",
"5. PROFIT"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4lOW5x/HvzWLZlG0QF+Rg1Qu0btTU2uPSqsW6oHjaHluxFUdaOAoWBRLQLqKtFkKN1hNUqDjiqYJVa0WhRUEFz1EREFEUEaSlqCwZBBfUIvCcP96JE7bMJJmZd/t9ritX8s6SuSeTued572cz5xwiIhJ+zfwOQERECkMJXUQkIpTQRUQiQgldRCQilNBFRCJCCV1EJCKU0EVEIkIJXUQkIpTQRUQiokUpHyyRSLgePXqU8iFFREJv0aJFaedcl1y3K2lC79GjBwsXLizlQ4qIhJ6Zrc7ndiq5iIhEhBK6iEhEKKGLiESEErqISEQooYuIRIQSuohIRCihi4hEhBK6iEgR/e//wh13QCl2+yzpxCIRkbjYtg2OOw7eeMM7HjAA2rYt7mOqhR4l6TSMH+99F5Gma+R76vHHoWVLL5l3Js3K/xpP20+L/75UQo+SVAoqKrzvQaYPHgmLBr6nPvsMOnaECy7wjk8/HWrGpTjsrtK8L1VyiZJkcufvQVX7JgEoL/c3FpH6NOA9de+9O99s8WI4/nggnQTL73c0lblSVOozysrKnBbn8kk67SXSZBISCcUiUiAffAAdOmSP+/eH++8v7GOY2SLnXFmu2+VVcjGzf5jZa2b2ipktzFzWycyeMrMVme8dmxp0YESxJBCkckwi4bXMlcwl5Cors8m8M2nSo8Zz/+/9yxsNKbmc7pyrG+loYI5zbqyZjc4cjypodH4JS0mgIS3dsJRjREJg3To48MDs8YgR8LuumbzRGd/yRlNq6P2Ab2V+ngI8S1QSeimTX1PKDw354KltFYtIk4wcCbfckj1euxYOOACvVg47540SlxfzTegOeNLMHDDROTcJ6OqcW5u5fh3QtRgB+qKUya8pZwNqdTed6vmSp7ffhsMPzx5XVu7ylt1T3ijx2X6+Cf0U59y7ZrY/8JSZvVn3SuecyyT73ZjZIGAQQPfu3ZsUbCQ1JSmr1d10YSmvSWE18IO8f3+YOjV7vHkztG+fx+OUuNGVV6eoc+7dzPcNwKPAicB6MzsQIPN9w17uO8k5V+acK+vSJeeWeNG2p85WdRDurNQd0smk19TSWU685DlIYPFiMMsm81TKm8KfVzKHkr+/c7bQzawt0Mw591Hm57OAG4HpwABgbOb7Y8UMNBLUGsyt1H8jneXEU46Ws3PepKC5c73jjh3hvfegVasSxddI+ZRcugKPmlnt7R9wzv3NzBYAfzKzgcBq4KLihRkRqnnnpr+RlEI9H+TPPusl81rTp8P555cmrKbSxCIREbzFtI46Clas8MaUX9s1xbBXkrQ4wP9yaEEnFolIRhQnnQmPPuotprVihXe84IoUI9ZX0OJ/AjARrwG0lotIQ6gfJFI++QT23x+2bPGO+/SBWbPANibhUEJX+lNCF2kI1fgj4+674ac/zR4vWQLHHps5CGlnuRK6SEOE9I0uWZs2QadO2eMrL0ozoSwFByUB/+vlTaEaukgxqeYeKL/97c7JfNUqvGQelIXrmijcLXRN25agU809EN57Dw4+OHs8ahSMHZs5iFAZLdwJXW8WCboIJYuwGjYMbr/d+7kzaf7xqxTtrqpTXolQGS3cCV1vFgm6CCWLsHnrLejZM3tcVQXXbMs0AtsRydcl3AldbxYR2YVzcNFF8PDD2cs++AD22489L3EbIeoUbQp1eIkEynvvwXe/C0seXs7jnMdfKpfjXCaZQ+QXw1NCb4ogbesmEmPOweTJ3tT9lTOX83KbU+nLTPo9O9zv0Eoq3CUXv6mGL+K7Vatg0CCYMwdOOw1m2HDaza2BLl28wnmMqIXeFBE/fRMJsu3b4bbb4Jhj4KWX4K674JlnoN3EKjj3XHjuuZ17RWNACb3YVGeXYoj5/9Xrr8PJJ8M113hL3b7xBgweDM2a4SXxGTNil8xBCb34VGeXYojp/9XWrfDrX0Pv3rByJdx/Pzz+OHTr5ndkwaAaerGpzi7FEMP/qwULYOBAeO01+OEPvclCcd/VcldqoReb6uxSDDH6v/rkE+9k5KSTYMeGNEsHjGfqf6eblswjWrJSQheRwJo7F447zsu9AwfCgiEpvjKlAKWmiJasVHIRkcD54AO48Wdp9r2vmpH7fUKfS9vw5ZuHAkloRdNLTREtWWlPUREJlBkzvBEr/d8bT6WryF5RWRnbpT7y3VNULXQRCYSaGrj6anjgATj6aPjhPUl4fotXRG/TJnKt6WJQQhcRX7maNEuuTvGDvyX5+0cJxoyBa6+FffZJwFlj/A4vVJTQRcQ3774LT/ZJkVxWweT2z9J5zhSOPDX6I3eKRaNcgiKiw6gkREr4P+gc/OEP3mJav/p7klVHnsspH8zkyBejNeqk1JTQgyKiw6gkREr0P7hyJZx5preg1unHpHllWIovP1rldXqqTt4kKrkERUOHUWk/VSm0Ig/lq11M65e/hJYtvRb6wPdT2KgK6ExsR7AUkhJ6UDR09yXtpyqFVsQdwJYuhcsv96bvn38+3HlnZtPmdBIMtcwLRAk9rCI6MUKiZetWuPlm76tDB5g2zdsezixzA20jWVBK6GGlN4IE3Esvea3y11+HSy7xyi2qDhaXErqIFNQnn8DYkWna3VnJHV9azPa7qjl9cPzWJvdD3qNczKy5mS02sycyx4ea2XwzW2lmD5rZPsULU2Jv1yF1GuYZSE8/7e0g9OmdKSoYz2n/ms3p0+O1r6efGjJscRiwrM7xOOBW59zhwCZgYCEDE9nJrkPqNMwzUDZvhp/+1BuO2KwZXPiXpFcS/Pa3Y7evp5/yKrmYWTfgPOAmYLiZGXAG0D9zkynAGODOIsQosnsnsDqFA2P6dLjiCli3zvuMHTMGWrdOQL9Kv0OLnXxb6LcBFcCOzHFnYLNzblvm+B3g4ALHVnw6bQ+PXTd0iOIGDyH7f9ywwds5qF8/6NwZ5s+HceOgdWu/I4uvnAndzPoCG5xzixrzAGY2yMwWmtnCmpqaxvyK4tFpuwRJSP4fnfP28jzw18cwu7sxps8BLPxbmrKci7tKseVTcjkZuMDMzsVbWn4/4PdABzNrkWmldwPe3dOdnXOTgEngrYdekKgLRaftEiQh+H9cs8Yrr8yYAVy/lOTzcP1T6+H+lIbRBkCDNrgws28BI51zfc3sIeAR59w0M7sLeNU5d0d999cGFyLhtGMHTJrknUBs3+5NFPpDy2NY94+lXLuiKyP+sDRa5a+AKcUGF6OAaWb2G2AxMLkJv0tEAmrFCvjJT2DePG/QyqRJcOihMIzX/A5NdtGghO6cexZ4NvPzKuDEwockIkGwbRucOnYIL34+kZadBjN58gSSyTrT9iVwtHyuiOxmyRI46CdDeHHbHdBsOzt6T+Tyy5XMg05T/0XkC//6F/zmN3DTK0NwJ9zhrYQIDD5hsL+BSV7UQhdpiJCNFW+IF16A3r29hM4JE79I5leWXcmE8yb4GttOIvwaNJUSuvgrbG/OfMeKh+h5ffwxXH01nHwybNkCf/0rXPG1wTS35sFL5hD88fp+vvbOuZJ9nXDCCS70amqcq6z0vkvTVVY6B973MMj39Q/J8zr/zisdv2ruOPdKN2SIcx9+6HdEeQj6e7AIrz2w0OWRYxs0Dr2pIjEOffx4r3VQWamJFIUQ1a30Av68Nm2CkSPhnm4toNl2mtGc7ddvy31Hya0Ir32+49CV0Bsq4G9UkVwefRSuvBJqauC4a4ewpOVEBp8wOHilFfmCErqI7GT9erjqKnjoITj+eJg8Gb76Vb+jknzkm9DVKSoScc7BfffBkUfCY4/BTTd528MpmUePEroER4hGhoTF6tVwzjkwYICX0Jcsgeuug5Yt/Y5MikEJXYIj6MPRQmTHDjht3BB6pJox6yTj1N8O4bnnoFcvvyOTYtJM0aCKY+drCJaPDYP+9w9h6lsTwXaAeX1kz2+dSLNm6vSMOiX0oKptrUJ8hkfW7kIkjfL553DLLTD104nQbDsAhuFwmrofE0roQaXWqjTARVOG8NCqibBwMF8+bDCrExqKGEcathh0uUovcSzNyBc++wx+/Wu4uYUmCEWZhi1GRa6OQnUkxtb//Z83nvzmm6HXx97aK/9VptJKnKnkEnS5Si8qzcTORx/Bv980hKVfmki74wcz6/YJnHXWBEDllbhTyUUkRGbNgkGD4J+XeSWW5tacbb9SiSXqVHKRnWnSTqhtnr+cV7qdx1VnL6dNG/iP7l6JRaNXpC6VXOIijsMgI+KRR2C//sPps3Um0w+DHotn0KqVSiyyOyX0uFCtPVzSaT78fYqfLU4yZUaCfkdWcWJn6HV3FbTyOzgJKiX0uNCkHU8Ihnm6mjTvfetiDn5jNoc128LYsWMYMaInLVrM8Ds0CTgldImXIJee0mnevyXFE3/awqWrZgNwxRWQGOVzXBIaSugSL6UqPTXwTGD7G8v54PQL6bThTd5peT0vnXM9ZV+DxFVDixunRIpGuUi81JaeEonijvxpwISvZctg/snD6bThTf7ZtheXvjSUE+8bSrN2bQsfl0SaWugSX8Usv+RxJvD52jTzkil+/HSSbq2r+PNxcMi0KqxXIrt3bTFik8hSQpf4Kmb5JUcn9KJFMO/8FNesraDyWOjzZDldu9bp9NSoJGkEzRQVKaFP16R55tIUyXlJunSBad9JcfQtwR1xI8GQ70xRtdBFSmTePHjheylGpSuo/prXKu/QQeUUKZycnaJm1srMXjKzJWb2upndkLn8UDObb2YrzexBM9un+OGKhM+HH8KQIfDNb8Kf2iRZ8dNK/nNmkg4d/I5MoiafUS7/As5wzh0HHA+cbWYnAeOAW51zhwObgIHFC1MknOZOWs6Crucx+47lXH01zHsjwRGTylVi2RutOdQkORO683ycOWyZ+XLAGcDDmcunABcWJUKRENq4ES69FD4aPJwzP5vJi98Yzq23QluNRKxfdbU3uqe62u9IQimvcehm1tzMXgE2AE8BbwObnXO163a+AxxcnBBFQiKdxpVXsO7YPpx3xHKmToW3r6hi+9nn0jFVldf91TqVpsirU9Q5tx043sw6AI8CvfJ9ADMbBAwC6N69e2NiFAmFD36fov3vxnMAML79cNovmsGxx/YE8lyDJcjLEpTK0KHeaYyGazZKg0a5OOc2m9kzwDeADmbWItNK7wa8u5f7TAImgTdssYnxigSOq0mz4MoUV/z1Avo3r+G7hy7mG3+posVXGviLNPZci8g1UT6jXLpkWuaYWWugD7AMeAb4fuZmA4DHihWkxFxQSxHpNBtHj+e+E6s58eEKftJ1Ov3erOTQFU/R4is9G/776i5LEARB/bvLXuVTQz8QeMbMXgUWAE85554ARgHDzWwl0BmYXLwwJdYCuBH29u3w3OUpOo+r4N334MXvVjL4hSSHH+53ZAUUwL+71C9nycU59yrQew+XrwJOLEZQIjsJWCni9ddh4EBYOT9JZS+47MEkBx0bkFZ1IQXs7y65aeq/RE8xNrFIp9n2hxS3bk7y81sT7Lcf3H47XHwxmBXmIUT2RptES3wVoVSw5sYULa6rYENliu9/31vytn//vSRz1Z7FJ1rLRaKngKWCTz6B66+HKdVJfrYffPuOJN+5JMedNPxQfKKELtHT1KFvy5fD8OHM/2EVP7qxJytXwqBBCa6qLKd9+zzur9qz+EQJXaSudJrt519I8xVvUjMT3GEzePppOP30BvwOjaUWnyihi9SxrCLFkSveZBm9WHpZFa9OgDZt/I5KJD9K6CJATQ1cfTXMeiDJtV3h9PuSjD4rgkMRJdI0ykVizdWkWdx/PKf0SvPQQ3DVmARX/bOcryqZSwiphS7xkxmnvujYJNPOTjGeCoYeAqfPLefoo/0OTqTxlNAldnZMTtFsdAXTgBRJOneC8gVJmnf1O7IIKsYkL9krJXSJh0xiefrfklw0OkkSL5k/MCvBWWdpRErRaEx+SSmhSyxs+4M30/OvwEbKmXdiORtegGbqRSoujckvqej9O2vadbQU4PWcOhUOuC5JOZWkSDJ/Psyfr2ReEkFbEjjiotdC1yletDTh9fzoI9hvv9qjBG//Rzk1j2gxLYmu6CX0uJ/iRa0TqpGv5+23w7Bh2eNly6BX3hsnioRT9BJ63KddR+0MpYGvZ00N7L9/9vjKK2HChCLEJRJA0UvocRfjM5Rf/AJuuil7vGYNdOvmXzxNFrWzLSk6dQtFTQw7oVav9uritcn8hhvAuZAncwjOFnAaaBAaaqFLqA0cCPfckz3euBE6dfIvnoIKytlW1Mp4EaaELuGTTrPutymOrkqyEe9M5K67YPBgn+MqtKD0BwXlg0VyUkKXUHEOJp+S4ifLK0gC//2lcjZuhLZt/Y4swoLywSI5KaFLaDz/PJx8MnQmyXLg1MlJxl/ud1QiwaFOUcmPHx1jmcfcvj7Nccd5yRxgv0MT3Ly1nAsuj0/Hr0g+lNAlP36MuMg85ugDUrz6qnfRnDmwahW0bFm6METCQiUXyU8pOsbqjLv+rF2CY3+bpB/eqoinnAJz52r9FZH6KKFLfkrRMZZpkS9aBGUPlgMJfkc5CxfCCScU96FFokDtHQmMj77vrYj4nQe9s4CLLoIdO5TMRfKlFroEQlUVjBiRALyzgLfegiOO8DcmkbBRCz0sIjr9ev16b9r+iBHe8bBh3lhzJXORhlNCD4ugrOtRQKNGwQEHZI/ffRduu82/eETCLmdCN7NDzOwZM3vDzF43s2GZyzuZ2VNmtiLzvWPxw42xZBIqK8M9/TpzlvGPhWnMvKcDcPPNXqv8oIP8DS9yInpWJ3uXTw19GzDCOfeyme0LLDKzp4DLgDnOubFmNhoYDYwqXqgxF4Xp15mzDG95cu+5bNoEHTr4GVSEaVGt2MnZQnfOrXXOvZz5+SNgGXAw0A+YkrnZFODCYgUp4bdkCSQqsvt63n231ypXMi+iKJzVhUGAzoQaNMrFzHoAvYH5QFfn3NrMVeuArnu5zyBgEED37t0bG6eESZ0JQq5zgj59vBmekOCuduVs2ACtW/sdZAxE4awuDAJ0JpR3QjezdsAjwNXOuQ+tzk67zjlnZm5P93POTQImAZSVle3xNhIx1dVwww2sfmMLPe4d88XFjz4KF+o8TqImQMsL55XQzawlXjK/3zn358zF683sQOfcWjM7ENhQrCAlXLbvgObAPfd6xz17wtKl0EKzHiSKAnQmlM8oFwMmA8ucc1V1rpoODMj8PAB4rPDhSdhMnw5dfz2UciqZwFCefRbefFPJXKQU8nmbnQz8GHjNzF7JXHYdMBb4k5kNBFYDFxUnRAmsOrXyT9sm6NoVPvoIIMGi08upmeNNGhKR0siZ0J1z/wvs7W15ZmHDkdBIp2HAAJg5kwUL4MSHsqecixfD8cf7GJvkp84HcpM3FS/k75JG00zRqGvIkKpct629fvnyL5L5E5zLOQ8l6Uya+48fj6tJK5mHRSFnH0dwJnMYqbIZdQ0ZUlXfbdNpuPhimD0bnnwSZs/mCc7lMqawkQTpUePpPK4CUnk8jgRDIUdnBGikR5wpoUddQ95o9d02lfKSOTBudm/SnEWKJMmRCcaPB9JJ6Jzn40gwFHJ0RoBGesSZOVe6oeFlZWVu4cKFJXs8KaB0mlnnV/PCi1DNUDaSYO3anRfXEpHiMLNFzrmyXLdTC11yWrkSjjgiAYwBvDL6yJG+hiQie6CELvW6+GKYNi17vHkztG/vXzwisnca5SJ7tHixN4a8Npnfe6+3mJaSuUhwKaHLTsMVd+yAb34TvvpV76qOHeHTT71RiiK7CdBKg6KELvDFcMWVP0/RvDnMm+dd/Pjj8P770KqVv+FJgGn8eaAoocfRLq2qz3+UZGynSk6a5A05PPpo+Pxz6NvXzyBjIAqtW625HijqFI2buhOEtmzhz8eO4XvfS1C7g9Bzz8Epp/gbYmwEaB3tRtP480BRQo+bOhOEfnMT/HKbd/FZZ8Hf/qbFtEpKsyulwJTQo27XRZOSSRbO3cITM6B621AAXn0VjjnG5zjjSK1bKTAl9KgbMwYmTIDVq3n/xmo6d8lOELrsMvVliUSJOkWjqG5n26xZAGycOovOnbM3+fvfI5zMo9DZKNIISuj1CWtiqK72Otuqq1l/y//wBr049/3/AeDaa70JQj16+BtiUWkoncSUSi71CfkohJkz4bwbTgKWAbB+Pey/v78xlYQ6GyWm1EKvT1jG2O5yJrHybG9Pz0sXeJ2et97qtcpjkcwh29monXMkZtRCr0/QRyHUjmDZsgVuuAHn4D9fKueRR7Ljyj/8EPbdN4/foa3DREJPLfQwqm2R19bKgTU/q6TLqCSPPOLd5I9/9Frl9SZzUL05l7D2o0gsqYUeRrVJ+Prr2TGuknMeTPLky17rumtXWL0avvSlPH+X6s31C3k/isSLWuhhUttavOACqKxk7jFDaT6q/Itk/te/wrp1mWSeb8syV705Ki3Uxj6PMPSjROU1CpKQ/k3VQg+6dNorrdS64Qa2bYMv31nOmjXeRb17w4IF0Lx5nfsVqmUZlRZqY59H0PtRIDqvUZCE9G+qhB50qRTccIP38/XXs+SSSs68LsnGzNUvvAAnnbSH+xWqlBKVkkxUnseeRPm5+SWkf1NtEh1EdUeeAFRXs3UrHDJuKBt2eKWRvn1h+nQtpiUSB9okOsx2Od27Y/8xDBmSvfr11+Goo/wJTQpIQ0alwJTQgyjTMn+/X5LOdVrggwbBxIk+xeSXqCa9dNrb12/mTO/YzzptVP/GMaRRLqWSq9e87vWJBGO2lNO5Z/bNtXp1DJM5RHecfCrlJfNzz/W/ThvVv3EMqYVeKrl6zTPXb94MHW/OXv+rX2X7RGMppJ1TOdV9Xn63iqP6N46hnJ2iZnYP0BfY4Jw7OnNZJ+BBoAfwD+Ai59ymXA8W607RXKe16TSP9E0xeH6SjXjX19T4/14PDZUNJMLy7RTNp+RyL3D2LpeNBuY4544A5mSOpT51J/AsXw7nned9B5YtA+uS4Pvzy9lIgupqb9q+8lIDBKVsENIJKRINOUsuzrl5ZtZjl4v7Ad/K/DwFeBYYVcC4omf5chg+HKqqvO8zZ+KAfs1n8Pjj3k2aN4fNm6Fduz3cXy3Q+gWlbBDSCSkSDY2toXd1zq3N/LwO6FqgeKJr6FBvc+atW6G6mk2b4Oszq1iRuXraNPjBD+q5f3W1V0zfssXbVm5P4pz0gzKjMygfLBJLTe4Udc45M9trId7MBgGDALp3797Uhwuv3r1h9mx2HNebr/XvycsvzwDgkENg5UrYZ58m/O5dltEFgpHc4igoHywSDCVuZDU2oa83swOdc2vN7EBgw95u6JybBEwCr1O0kY8XfhUVvPl+F065JTtt/8knoU+fPO8/dCi0bbvnll+d1RcDv5CUSJyUuATX2IQ+HRgAjM18f6xgEYXdHj6R//Uv6HFMgnXrvBf061+H55+HZg2ZBVBfyy9IQ+BEJKvEJbicKcXMpgIvAD3N7B0zG4iXyPuY2Qrg25ljgd1GWzzwALRq5S1rC/DSS/Diiw1M5rloyzWRYCrxezOfUS4X7+WqMwscS3jVHcGS+ST++D+T7Ftn2v53vwsPP6zFtERiJSQ1dIHsi/Xkk94IFoAZM/j9PuVcfWj2Zm++CT17+hOiiPgoJDV0geyLVV4O++zDxl9UkajTAr/ySpgwwb/wRMRnJa6hK6E3Rm3L/IILvONkkp/fmuDmf8/eZM0a6NbNn/BEJCBKPIxVCT1fdWthdU6jVl9UTo8u2ZvdeCP88pc+xSgisaaEnq+6tbDM6dOwV5LcXpG9ycaN0KmTD7GJBE2cZy37SOuh51K7kNapp34xaWfpugRWUc7tD3j/qHfd5S2mpWQeMVpoq/GCslhazKiFnktmIS0A98QMzjkHZs3yrmrVymuVt2njY3xSPFpoq/G0po0vlNBzqaoCYNElVZTVOZ95+GH43vd8islvcTmdVlJqPK1p4wsl9By2H96T3mtm8Nol3vGXv+yNK2/Z0t+4fBWXlquSkoSMauiw11rpjBnQogW89pp3PGcOvP12zJM5eC1WLQImEjhK6LBbB85nn3mNs759vatPPRW2b4czzvAxxiDR2jHiF3VU10sJHXZqcd53H7Ru7XV2AixaBPPmFXgxLQkmJYvg0+iZeqmGDpBI8MGgcjp0yF70gx/A1KkRXkwrLh2b9dn1bxCXvoEwU0d1vZTQgVtugZEjs8dvvQVHHOFfPCWh5LX730DJIvjUUV2vWCf09evhgAOyx8OGwW23+RdPSSl57f43ULKQkDPnSrcrXFlZmVu4cGHJHq8+o0Z5ZfNa770HBx7oXzwiIntjZoucc2W5bhe7rr5Vq7y6eG0yHzvWm7avZC4iYRfNhL6X0Qo//jEcdlj2eNMmr6UuaISHSAREM6HvMrRpyRKvVf7HP3pX33231yqvO6ol9jQcbHf6kJOQiWanaKaTy12W5Mwz4JlnvIv33dfrCG3d2sfYgkqdpLvzeySQhpZKA0UzoScSzPt6Od/cP3vRo4/ChRf6F1LgaYTH7vz+kPP7A0VCJ7wJfS+tl23b4OijvWXMwduceelSb02W0FELzV9+f8j5/YEioRPeGvoear5/+Yu3cFZtMp8711sZMZTJHFTXjjutmSMNFNZUt1Pr5dNPYf/94eOPvYvOOANmz47AtH210ESkAcLRQt/TaINM6+We6QnatMkm81de8Za5DX0yB7XQRKRBwtFC30Pn0ObN0LFj9iaXXJIdligiEkfhSOi7lB7GjYPRo7NXv/22t5OQiEichSOhZ0oPa9fCQXVKKSNHepUYEREJS0IHrrlm55UQ162Drl39i0dEJGhC0Sk6ZEg2mf/ud960fSVzEZGdNSmhm9nZZrbczFaa2ejc92icvn3htNO8jtARI4r1KCIi4dbohG5mzYE4c5U9AAAD4ElEQVQJwDnAUcDFZnZUoQKr65xzvElC7dsX47eLiERDU1roJwIrnXOrnHNbgWlAv8KEJSIiDdWUhH4wsKbO8TuZy3ZiZoPMbKGZLaypqWnCw4mISH2K3inqnJvknCtzzpV16dKl2A8nIhJbTUno7wKH1DnulrlMRER80JSEvgA4wswONbN9gB8C0wsTloiINFSjJxY557aZ2VBgFtAcuMc593rBIhMRkQZp0kxR59xMYGaBYhERkSYIxUxRERHJzZxzpXswsxpgdSPumgDiuPV6HJ+3nnM8xPE5Q+Of978553IOEyxpQm8sM1vonCvzO45Si+Pz1nOOhzg+Zyj+81bJRUQkIpTQRUQiIiwJfZLfAfgkjs9bzzke4vicocjPOxQ1dBERyS0sLXQREckh8Am9VJtoBIWZHWJmz5jZG2b2upkN8zumUjGz5ma22Mye8DuWUjGzDmb2sJm9aWbLzOwbfsdUbGZ2TeZ/e6mZTTWzVn7HVAxmdo+ZbTCzpXUu62RmT5nZisz3joV8zEAn9FJuohEg24ARzrmjgJOAITF4zrWGAcv8DqLEfg/8zTnXCziOiD9/MzsY+BlQ5pw7Gm/ZkB/6G1XR3Aucvctlo4E5zrkjgDmZ44IJdEInhptoOOfWOudezvz8Ed4bfLd15qPGzLoB5wF3+x1LqZhZe+A0YDKAc26rc26zv1GVRAugtZm1ANoA7/kcT1E45+YB7+9ycT9gSubnKcCFhXzMoCf0vDbRiCoz6wH0Bub7G0lJ3AZUADv8DqSEDgVqgFSm1HS3mbX1O6hics69C/wO+CewFvjAOfekv1GVVFfn3NrMz+uAgm53H/SEHltm1g54BLjaOfeh3/EUk5n1BTY45xb5HUuJtQC+CtzpnOsNbKHAp+BBk6kZ98P7MDsIaGtmP/I3Kn84b4hhQYcZBj2hx3ITDTNriZfM73fO/dnveErgZOACM/sHXlntDDP7o78hlcQ7wDvOudozsIfxEnyUfRv4u3Ouxjn3OfBn4N99jqmU1pvZgQCZ7xsK+cuDntBjt4mGmRleTXWZc67K73hKwTl3rXOum3OuB95r/LRzLvKtNufcOmCNmfXMXHQm8IaPIZXCP4GTzKxN5n/9TCLeEbyL6cCAzM8DgMcK+cubtB56scV0E42TgR8Dr5nZK5nLrsusPS/RcxVwf6bBsgpI+hxPUTnn5pvZw8DLeCO6FhPRWaNmNhX4FpAws3eA64GxwJ/MbCDeyrMXFfQxNVNURCQagl5yERGRPCmhi4hEhBK6iEhEKKGLiESEErqISEQooYuIRIQSuohIRCihi4hExP8D9nfxQKX36fEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def separate(data, coefs, p):\n",
" x, y = data[:,0], data[:,1]\n",
" a, b = coefs\n",
" y1 = a*x+b\n",
" dist = abs(y-y1)\n",
" return data[dist <= p], data[dist > p]\n",
"\n",
"def ransac(data, iters, sample_size, p):\n",
" best_cnt = 0\n",
" best_coefs = None\n",
" best_inliers = None\n",
" best_outliers = None\n",
"\n",
" for i in range(iters):\n",
" idx = np.random.randint(len(data), size=(sample_size))\n",
" sample = data[idx,:]\n",
" coefs = regression_line(sample)\n",
" inliers, outliers = separate(data, coefs, p)\n",
" if len(inliers) > best_cnt:\n",
" best_cnt = len(inliers)\n",
" best_coefs = coefs\n",
" best_inliers = inliers\n",
" best_outliers = outliers\n",
" \n",
" return best_coefs, best_inliers, best_outliers\n",
"\n",
"\n",
"coefs, inliers, outliers = ransac(data, 50, 5, 0.1)\n",
"#plt.plot(data[:,0], data[:,1], 'ob', markersize=1)\n",
"draw_line(data[:,0], coefs)\n",
"plt.plot(inliers[:,0], inliers[:,1], 'og', markersize=2)\n",
"plt.plot(outliers[:,0], outliers[:,1], 'or', markersize=1)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment