Skip to content

Instantly share code, notes, and snippets.

@jclevesque
Last active July 1, 2022 02:26
Show Gist options
  • Save jclevesque/5273ad9077d9ea93994f6d96c20b0ddd to your computer and use it in GitHub Desktop.
Save jclevesque/5273ad9077d9ea93994f6d96c20b0ddd to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import scipy.odr as spodr\n",
"from sklearn.linear_model import LinearRegression"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def odr_fit(x_o, y_o):\n",
" def linear_func(coefs, x):\n",
" w, b = coefs\n",
" return w*x + b\n",
"\n",
" # Create a model for fitting.\n",
" linear_model = spodr.Model(linear_func)\n",
"\n",
" # Create a RealData object using our initiated data from above.\n",
" data = spodr.RealData(x_o, y_o)\n",
"\n",
" # Set up ODR with the model and data.\n",
" odr = spodr.ODR(data, linear_model, beta0=[0., 1.])\n",
"\n",
" # Run the regression.\n",
" out = odr.run()\n",
"\n",
" # Use the in-built pprint method to give us results.\n",
" # out.pprint()\n",
"\n",
" w, b = out.beta\n",
" return w, b"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1072c00b8>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4k1X7wPHvabpLoS2ljEIBWbL3kiGoDF9Z4h6vIjIFFV7hhygUZAiKggtUQMCFKAqIggIKiKBM2RuZLS2FLqC7yfn9kRaa0pHStEna+3NdXjTPyklqn/s56z5Ka40QQgiRycXeBRBCCOFYJDAIIYSwIIFBCCGEBQkMQgghLEhgEEIIYUECgxBCCAsSGIQQQliQwCCEEMKCBAYhhBAWXO1dgNsRGBioa9SoYe9iCCGEU9mzZ88VrXWF/I5zysBQo0YNdu/ebe9iCCGEU1FKnbPmOGlKEkIIYUECgxBCCAsSGIQQQlhwyj6GnKSlpREWFkZycrK9i2JXnp6eVK1aFTc3N3sXRQjhpEpMYAgLC8PX15caNWqglLJ3cexCa010dDRhYWHUrFnT3sURQjipEtOUlJycTPny5UttUABQSlG+fPlSX2sSQhROiQkMQKkOCpnkOxBCFFaJaUoSorRatTecWeuOczEuiSp+XoztUY9+zYPtXSzhxEpUjcEZ/fnnnzRs2JBmzZqRlJREREQEvXr1yvOcn3/+mdDQ0GIqoXBkq/aGM37FQcLjktBAeFwS41ccZNXecHsXTTgxCQx29vXXXzN+/Hj27duHl5cXs2fPZvDgwXme88ADD/DTTz+RmJhYTKUUjmrWuuMkpRkttiWlGZm17ridSiRKAgkMNhIaGsp777134/Xrr7/O+++/n+c5Cxcu5LvvvmPixIk89dRTAPzwww/07NkTgDlz5jBw4EAADh48SKNGjUhMTEQpRZcuXfj555+L6NMIZ3ExLqlA24WwRsnsY/jlVYg8aNtrVmoM98/MdffAgQPp378/o0aNwmQysWzZMjZu3EizZs1yPH7p0qUMGjSIrVu30qtXLx5++GHOnDmDv78/Hh4eALz88st06dKFlStXMn36dD799FO8vb0BaNWqFX/++SePPvqobT+ncCpV/LwIzyEIVPHzskNpRElRMgODHdSoUYPy5cuzd+9eLl26RPPmzalevTr79u2z+hoRERFUqHAz8aGLiwtLliyhSZMmDB06lA4dOtzYFxQUxMWLF236GYTzGdujHuNXHLRoTvJyMzC2Rz07lko4u5IZGPJ4si9KgwYNYsmSJURGRjJw4ECuXbtGp06dcjx26dKlNGjQwGKbl5fXLXMQTp48SZkyZW4JAsnJyXh5yVNhaZc5+khGJd0ko7QKr2QGBjt58MEHCQ0NJS0tjaVLl2IwGApUY6hbty5nz5698To+Pp6XXnqJLVu2MHLkSL7//nsefvhhAE6cOEGjRo1s/RGEE+rXPFhufBkyR2ll1qAyR2kB8h0VgHQ+25C7uztdu3bl0UcfxWAwFPh8Hx8fatWqxalTpwAYPXo0I0aMoG7dunz22We8+uqrREVFAbBp0yYeeOABm5ZfCGcno7RswyY1BqXUIqAXEKW1vuUxVin1FDAOUMA1YLjWen/GvrMZ24xAuta6lS3KZA8mk4nt27ezfPlyq89ZsmSJxeuRI0eyZMkSpk2bxqJFi25sr1at2o2AcenSJZKSkmjcuLFNyi1ESSGjtGzDVjWGJUDPPPafAe7WWjcGpgLzs+3vqrVu5sxB4ciRI9SuXZt7772XOnXq3PZ1HnzwQfJbtvT8+fO8++67t/0eQpRUuY3GklFaBWOTGoPWeotSqkYe+//K8nI7UNUW7+tIGjRowOnTp21yrUGDBuW5v3Xr1jZ5HyFKGhmlZRv26GN4Hvgly2sNrFdK7VFKDbFDeYQQJUS/5sHM6N+YYD8vFBDs58WM/o2l47mAinVUklKqK+bA0DHL5o5a63ClVBCwQSl1TGu9JYdzhwBDAEJCQoqlvEII5yOjtAqv2GoMSqkmwEKgr9Y6OnO71jo8498oYCXQJqfztdbztdattNatsk4CE0IIYVvFEhiUUiHACuC/WusTWbb7KKV8M38GugOHiqNMQgghcmaTwKCU+gb4G6inlApTSj2vlBqmlBqWcUgoUB6Yp5Tap5TanbG9IrBVKbUf2Ams0Vr/aosyOQtJuy2E7azaG06HmRup+eoaOszcKOnHb5PSWtu7DAXWqlUrvXv3bottR48epX79+nYq0e0bNmwYHTt25OmnnwZg7NixdOzYkb59++Z6jtaaFi1asG3bthtJ9bJy1u9CiMLIPusZzCOSpPP5JqXUHmumBcjMZxvZtWsXTZo0ITk5mYSEBBo2bMihQ3m3iuWXdrtz584WKTU6duzI/v37Je22EDmQWc+2UyJzJb218y2OxRyz6TXvDLiTcW3G5bq/devW9OnThwkTJpCUlMTTTz9N9erVC5V2+/nnn2fJkiW89957nDhxguTkZJo2bQpI2m0hsrPnrOeSlrivRAYGewkNDaV169Z4enrywQcfFDiJXva024888ghTp05l1qxZLFq0iAEDBtzYJ2m3hbBkr7UpSmLivhIZGPJ6si9K0dHRXL9+nbS0NJKTkzGZTIVKu+3t7U23bt348ccf+e6779izZ8+NfZJ2WwhL9pr1nFcTlgQGwdChQ5k6dSpnzpxh3LhxfPTRR4VKuw3m9Bi9e/emU6dO+Pv739guabeFsGSvtSlKYuI+CQw28sUXX+Dm5saTTz6J0WjkrrvuYuPGjdxzzz1WXyNr2u3atWsD0LJlS8qWLctzzz1nceymTZuYMWOGTT+DEIVl77Z2e8x6LonLq8qoJBt55pln+OGHHwAwGAzs2LHDqqCwZMmSG4vvwM2025kuXryIyWSie/fuN7ZJ2m3hiDLb2sPjktDcbGsv6XMJxvaoh5eb5forzp64TwKDg8madvuLL76gbdu2TJ8+HReXm78qSbstHJEzDRe15US4kpi4Tya4lUDyXQh7qPnqGnK6myjgzEzHWW2wNE+EkwluQohi5SyL5DhTzcZeJDAIIWzCWdraS+IoIluTwCCEsAlnaWt3lpqNPclwVSGEzTjDIjmy/Gf+JDAUo8mTJ1OmTBnGjBlj76IIUWrZayKcM5HA4ADS09NxdZVfhRDFxRlqNvYkfQw2NHv2bBo1akSjRo147733AJg+fTp169alY8eOHD9+c9RDly5dGDVqFK1ateL999+3V5GFEOIWJfIxNfLNN0k5atu02x7176TSa6/lun/Pnj0sXryYHTt2oLWmbdu2dOrUiWXLlrFv3z7S09Np0aIFLVu2vHFOamoq2edjCCGEvZXIwGAPW7du5cEHH8THxweA/v37s2bNGh588MEbq6z16dPH4pzHHnus2MsphKOxd36l/KRHR5N67jzeLZrbuyjFpkQGhrye7B1JZhARorRy5LUMtNZcXb2aS2/OQHl7U3v9OpSbm13LVFykj8FGOnXqxKpVq0hMTCQhIYGVK1fywAMPsGrVKpKSkrh27Ro//fSTvYsphEOx5yzkvPIlpYWHc2HwEC6OexX3mjUJWTC/1AQFKKE1Bnto0aIFAwYMoE2bNoB5HYWWLVvy2GOP0bRpU4KCgmjdurWdSylE4VjT7FOQpiF7zULOtaZiMtH5yB9EzZkDQMXXX8f/ySdQBkNelytxbJJETym1COgFRGmtb1k9RimlgPeB/wCJwACt9T8Z+54FJmQcOk1r/Xl+7ydJ9PIm34UoCtYknytogroOMzfmuJZBsJ8X2161fi2Tgsrpfatdu8T/HfyB2lGn8enUicqTJ+EW7Dh9HbZgbRI9W9UYlgAfAV/ksv9+oE7Gf22Bj4G2SqkAYBLQCtDAHqXUaq11rI3KJYSwEWuWsCzoMpf2moWctUbiakrn4ZObefL4BpIMHlR5+y02VW3BrC+PczFuH1X8vOh6ZwU2HbvssB3ktmaTwKC13qKUqpHHIX2BL7S5erJdKeWnlKoMdAE2aK1jAJRSG4CewDe2KJcQwnasafYpaNNQ5s317XVHiEwMp7J3SLHcdDNXXasbe55Re5dT82oEm4Ob8WPHxxlWreUtzUxfbT9/41xH6iAvKsXVxxAMXMjyOixjW27bb4vWGnOrVenljOtrCOdgzRKWt7PM5QNNgvgzfjZ/R/zNj/1WE+gVaJsC5+H/7g7h2Ix3eeDEH8R6lmVy2+c4ENKEGf0a51jryS6vWlBJ4DSjkpRSQ5RSu5VSuy9fvnzLfk9PT6Kjo0v1jVFrTXR0NJ6envYuiiiBrEmrXdDU2ynGFEZtHsX6c+sZ2mRosQSFhL//puHE4fQ5sZk/63Zg2L1jCKvf6kY/iLUd3yU5TXdx1RjCgWpZXlfN2BaOuTkp6/bNOV1Aaz0fmA/mzufs+6tWrUpYWBg5BY3SxNPTk6pVq9q7GCWCo0+8Km7WJJ8rSIK6xLREXtr0EjsjdjKx3UQerfdokZbfGB/PpbfeJn7FCtyrVyfki88Z3qYNw7Mco7XG39udmMTUfK9XktN0F1dgWA2MVEotw9z5HK+1jlBKrQPeVEr5ZxzXHRh/O2/g5uZGzZo1bVNaUeo58sQre7Im+Zw1x1xLvcYLv73AgSsHmN5xOr1r9bZlMW9xdf16IqdOxRgTS/nBgwgcMQKXbDXriPgkQn88TExiKq4uinRT7q0PJT1Nt00Cg1LqG8xP/oFKqTDMI43cALTWnwBrMQ9VPYV5uOpzGftilFJTgV0Zl5qS2REthD0VdHSNsF5cchxDfxvKidgTvHP3O3Sr3q3I3istKopLU6dxbcMGPBrUp9onn+DVsKHFMUaT5qvt55i17jjpJhOv/edOyvt4MHvDCcLjkjAohVHrG/8GF7D26Iw1T1uNSnoin/0aGJHLvkXAIluUQwhbkeUfi8aVpCsMXj+YC9cu8H7X9+lctXORvI/WmvgffuDS27PQKSlUeOV/lB8wIMfZy2/9eoz5W07TqU4g0/s1JqS8ObeZwUVZ1BqNWt+oKRQkKDhjzVNmPguRg9sZXSPyFnE9gsEbBhOVGMW8e+fRpnKbInmf1AsXiAgNJfHv7Xi3akWlqVPwyNbMnJxmJCElnfJlPPhvu+rUr+xLv2bBFqMabVFrdNaap9OMShKiODnLwvbO4vzV8wz4dQAxSTHM7za/SIKCNhqJXrSY0737kHzwEJUmTyLki89vCQp//XuFnu9tYczy/QBUC/DmweZVbxnqbotao7PWPKXGIEQOZPlH2/k37l8Grx9MmimNhT0W0qB8A5u/R/LxE0RMmEDywYOU6dqVSpNCcatUyeKY2IRU3lx7lOV7wqhe3ptBne7I85q2qDU6a81TAoMQuZDlHwvvaPRRhm4YisHFwOIei6ntX9um1zelpnLl44+JXrAQQ9myBM+ZjW/Pnrc8/e89H8ugz3cTn5TG8C61ePneOni65Z0YzxbpOuyV8qOwJDAIIYrE/sv7Gb5hOGXcy7Cw+0JCyobY9PqJ/+wlYsIEUk+fplzfvgS9Og5Xf3+LY0wmjYuL4o4KZWge4s8r3etSv3JZq65vi1qjs9Y8bZJdtbjllF1VCOE4dkXuYsTvI6jgVYEF3RdQpUwVm13beD2By3PmELt0Ka6VK1H5jSmU6dTR4ph0o4lF287w66FIvh3aHjeDdKdC8WdXFUIIALaGb2XUplFULVOVBd0XUMG7gs2ufX3LFiImTSY9MhL/p54iaPQoXLKthHgwLJ5XVxzg8MWr3Fc/iISUdPy83W1WBltw9LkNEhhEqeXof5yOxprv6/dzvzNmyxjq+NXh026f4u/pn8vVCiY9NpZLM2ZwdfVPuNeqRfWlX+Pd3HIN5uQ0I++uP85nW89QvowH855qwf2NKjlcYk1nmNsggUGUSs7wx1lcrF2VLb/v66d/f2Litok0CmzEvPvmUdbdurb8vGitubp2LZemTcd47RqBLwyn/LBhuLjfWgMwuCi2norm8TYhjOt5J+W8bLcU56p/wpi1/oRNHiKcYW6DBAZRKjnDH2dRygwG4XFJKMyrZEHuATK/72v5ieVM/XsqrSu15sN7PsTbzbvQZUyLiCDyjSlc37wZz8aNCZk2Dc96dS2OuXwthQ9+P8nYnvUo6+nGyhfuyne0UYFEHeXM2tl4nD5PeOpLQOEfIpxhboP0yIhSyRn+OItK5tN/5vj67MNPMm/4WeX1fX155Eum/D2FjsEdmXvv3EIHBW0yEfvNN5zu1ZuE7dsJGjeOGsu+sQgKWmu+23WB+2b/wbe7LrD7rDnFmk2CgskIx3+Bz/vAvHZUObuKOJMXLphuHJLTd2St3OYwONLcBqkxiFLJWSce2YI1C9FkDwS5fV/lq/7J27vW0K16N97q9BZuhsI136ScOUPExIkk7d6Dd/t2VJ4yBfdq1SyOOX35Oq+vPMTfp6NpXcOfGf0bUzvIt1DvC0ByPOz9CnbOh9izUDYY7g2l/ZoqxHBrs9jtPkQ4w9wGCQyiVHKGP86iYs0NLXuAvPX70vhU3ECK70aa+t3L9u33U/fX9bfd/q7T0ohevIQrH32E8vSk8vRplOvfP8eO4xm/HOPQxXimP9iIJ1qH4OKSe+eyVQMMrpyEHZ/CvqWQlgDV2sF9k+HOXmBww2vbRrDhQ4QzzG2QwCBKJWf44ywquT39Z8opQFp+XwkEhPxKqs8WWgf8h7933k1Smnlhm9tpf086fJiICRNJOXoU3+7dqTjhddyCgiyO2XMuliBfD6oFeDOlb0MMShFUNu+VCvPsMG9aGf79HXZ8Aqd+A4M7NHoI2g6FKpajnYriIcLRZ9XLBDchSpnsN0zgRgd0fmsNGE1GpmyfwoqTK3imwTOs/L0FF+OSbzku2M+Lba/ek2c5TMnJXPnoI6IXL8EQ4E+l0FDKdrNcm+Fachqz1h3ny+3n6N+8Ku8+2tTqz9lh5sZbAqAPSTxfZjv/K7cJok9BmYrQehC0fA7K5D7foqQMbZYJbkKIHN1ubSnNlMbrW1/nlzO/MLTJUEY0G8G8H9bmeGx+zVUJO3cSOTGU1HPn8HvkYYLGjMFQrpzFMesORzLpx8NcupbMs+1rMKaAT+hZyxCiLvGsYT2PGDZTNj0JPFtC/4XQoC+45j/5zdGf8G1NAoMQpVBBb3SpxlTG/DGGTRc2MbrlaAY2GggUvBPfeO0aUbPeIe6773CrVo2QJYvxadfuluO+3nGO11ce4s5Kvnzy35Y0q+ZndVlvlKGcJzWu7eI5w6/c47IPIy6sMbVljVdfFgweVuDrlSYSGIQQeUpKT2LUplH8dfEvxrcZz5P1n7yxryDt79c2biRy8hukX7lCwHPPUeGlF3HxuhlATCbNlYQUgnw96dWkCkmpRp69q0bB8xylJsD+ZfziPpey7v9yRZflQ2M/vk6/j2tugcy4v3GJaRoqKhIYhBC5up56nZEbR7I3ai9T7prCg3UetNhvTbNU+pUrRE6fzrVffsWjbl2qzv0Ir8aNLa5z4tI1xq84SEJKOj+92JFyXm75rpdwi9hzsGsB/PMFJMdTtnJT9tR+kzFHanE23kgVPy9mZAQsmfWeN+l8FkLkKD4lnuG/Dedo9FFmdJpBz5o9C3S+1pr4H3/k0oyZ6MREczqL559HZUlnkZxmZN6mU3z8x7/4eLgy4YEGPNQi2Pr8RlrD2a3m0UXH1wIK6veCtsMhpB3kcJ2cOqXBug5zZyedz0I4CUds1ohOimbIhiGciT/D7C6z6RrStUDnp4aFEzlpEgnbtuHVvDmVp03Fo1Yti2MuxCTy7KKdnL6SQL9mVZjYqwHly3hY9wZpSXBwuXn+waVD4BUAHV42jzAqVzXPU0vzrHdr2SQwKKV6Au8DBmCh1npmtv1zgMz/s7yBIK21X8Y+I3AwY995rXUfW5RJCGfgiMn8LiVcYvCGwURcj+Cjez/irip3WX2uNhqJ/fprot57HwVUnDAB/yefQLnc7CfQWqOUolI5T2oFlWFyn4Z0rmtlau74MNj1GexZAkkxULER9P4AmjwKbtZNOCvNs96tVejAoJQyAHOBbkAYsEsptVprfSTzGK316CzHvwhknUGSpLVuVthyCOGMHC2ZX9i1MAatH0RcShyfdPuElhVbWn1uyqlTREyYSNK+ffh06kTlyZNwC775GbTW/Hwggk/++Jelg9tRzsuNBc/k26phbi66sMPcXHRkNaCh3n+g7TCo0THH5qK8lOZZ79ayRY2hDXBKa30aQCm1DOgLHMnl+CeASTZ4XyGcniM1a5yJP8Pg9YNJSk9iYfeFNApsZNV5OjWVKwsWcOWTTzH4+FDl7bco27u3RT9BWGwiE1cdYtPxyzQOLkd8Ylr+abHTU+DQD+aAELEfPMtB+xHm5iL/6rf9OUvzrHdr2SIwBAMXsrwOA9rmdKBSqjpQE9iYZbOnUmo3kA7M1FqvskGZhHAKjtKscSL2BIPXDwZgUY9F1Auw7uk5af9+czqLkycp26sXFV8bj2tAwI39JpNm0bYzvLv+BAATHqjPgLtq4JrXENRrkbB7kfm/hMsQWA96zYEmj7HqcByzPj3OxbhDhbqhl7YJawVV3J3PjwPfa62z1p2ra63DlVJ3ABuVUge11v9mP1EpNQQYAhASYttFxYUoKvl1LDtCs8ahK4cY9tswPAweLOy+kJrlauZ7jikxkcvvf0DMl1/iGhRE1U8+xrdLl1uOUwq2nrpC2zsCmNavEVX980jJHbbbXDs4vNKc+rpuD3Puoju6glIO2R9THOwxOMEWgSEcyJoXt2rGtpw8DozIukFrHZ7x72ml1GbM/Q+3BAat9XxgPpiHqxa61EIUMWtuZPZu1thzaQ8jfh+Bn4cfC7svpKpv3iN6ABL++ouI0EmkhYXh98TjBL3yCoYyZW7sT0o18sHGkzzZJoRqAd7Me6oFXm6GnIegpqfCkR9hx8cQvgc8ykKbIebmovKWo5gcrT+mONgrGNoiMOwC6iilamIOCI8DT2Y/SCl1J+AP/J1lmz+QqLVOUUoFAh2At21QJiHsztobmb2aNf66+Bcvb3yZSj6VWNh9IRV9KuZ5vDEujktvvU38ypW416hB9a++xLuVZefxlhOXeX3VQS7EJFGprCfP3lUDb/ccbjPXo8wji3Z9BtcjoXxt+M870PRx8Mh5bQVH6o8pLvYKhoUODFrrdKXUSGAd5uGqi7TWh5VSU4DdWuvVGYc+DizTljPq6gOfKqVMmFeTm5l1NJMQziy31NaOcCPbdH4Tr/zxCjXL1WR+t/mU9yqf67Faa66tW0/ktGkYY2MpP2QIgSNewMXj5pyD6OspTFtzlJV7w7mjgg/LhrSj3R05XPPiXvPcg0M/gDEVat0LfT8y/+uSd+oLR+mPKU72CoY26WPQWq8F1mbbFprt9eQczvsLaJx9uxDObtXecIu1lLOy943slzO/MP7P8TQo34CP7/uYch7lcj027VIUkVOncP233/Fs0ICQBfPxrF//luPmbf6Xnw9c5KV7avNC19qWS2wa0+HYT7D9E7iwHdx8oMWz5v6DwDpWl9sR+mOKm72Cocx8FuI25dUpOGvd8RyDggK73shWnlzJpL8m0TyoOXPvnUsZ9zI5Hqe1Jm75cqJmvYNOTSXqycFMcG1M2OenqeIXwdge9WhWzY+kNCP1K5fl5fvq8FjratStmKUZKCEa/vkcdi2Eq+HgXwN6vAnNnzYPPS0ge/fH2IO9gqHkShLiNuS02I2Xm4EZ/RvTr3kwNV9dk2NgADg784HiKWQ2S48uZcbOGdxV5S7e6/oeXq45P3WmnjtHROgkEnfswLt1aw49MYIxf8dafFZXF4VS0KyaH8uHZZsZHXnQ3Fx0cDmkJ8MdXcyT0ep0BxcDomBsOSpJciUJm3DEPD6OIL9OQT9vN2IT0245z987n0ldReR/6+awIXIR6dcacHDPQ6wLjLnl96jT04n5/Asuf/ghytWVSpMn4/foIzz59uZbPmu6SePp5sJHT7YwbzAZzUnstn8C57aCqxc0fcLcXBR0a9OTsJ49BidIYBC5Kq3jxq2RX6dgbhXx4q6ga6158dcZ/BH1DWnxTUm++ChJpN3ye0w+fpyI114n+fBhytxzD5UmheJW0TxKKbfPmpJmoqJbEmybDzsXQPx5KBcC3aaam4u8A3I8Tzi+Aq6AIUqTvJ6KS7vcOv8yt8cn3VpbyGt7UdBaM2v3LP6I+obUuFYkX3wM88DBm79HU0oKUe+9x+n+DxH173lmtH6ax6v1Z83F9BvXqVjW85Zr11ZhzPFZArMbwIaJ4BcCj30FL+2FDi9JUHByUmMQuSop48aLojksv05Bew+tNGkT07ZPY/mJ5aTGtCflUm+yPwf6/XuEM/3eIfXMGTbVaMPHDXpx3d0b4pMZv+Ig8Ulp7DgTTWJqOp6uLqSmp9PVZS8DDOvoZDiEUbtDo0fN/QeVZHBhSSKBQeTK3jc3Wyiq5rD8RsjYc2hluimdidsm8vPpn3m+0fN8t6ExF0m+sd87LZlnj6yl15m/0cFVmN1tJBt8alhcIynNyOTVh3FzdWFM50q0ifuFCkeXEKwvEUV5DtcfRcNeL4FP7vMfspP+KuchgUHkqiSMGy/KmaN5dQoWdGilrW6aacY0xv05jg3nNvBi8xcZ0mQINVxuBsfWkUcZuf8HApPiiev5IO2nv8ZvUzfneK2a6iIrmx+m3O7lkHodqrWFtjMJqt+bIEPBOtGlv8q5SGAQuSoJ48bt2Rxm7WgSW900k9OTGb15NFvDtzK21VieafjMjWsYrsZx5a2ZtDu9m3C/Slx4dQo9H+sGWNYMFSY6uxzkOcOvdDHshyPu0LC/eXRRcIsCff6srA3QUqtwDBIYRJ6cPT2xMzSH2aJWk5iWyIsbX2RX5C5C24fySN1HAHMH9NWf11DvzTepff06gSNGUG/oEFyyrLvcv0UwizcdpL/LFp41rKeWSwSXtR9H7xxJ/V4vQ5mgQn9GawK01CochwQGUaI5Q3NYYWs1V1Ov8sJvL3DoyiGmd5xO71q9AUiLiCBi8mQS/tiCZ9MmhEydimfdujfOi09KY/6q3wk4vIS/PTbjSxL7TLWY7DYWAHehAAAgAElEQVSa5j2fpW/L/NNvW8uaAF0as6c6KgkMokRzhuawwtRqYpNjGbphKCfjTvLO3e9wX/X70CYTscuWcfmdd9FaU3H8q/g//TTKkDHrWGt2/r6C5G3zeMW0B+3mgq7fD+56gWZVW1EU6+xaE6BLyii4kkACgyjxHL057HZrNZcTLzN4/WDCrofxQdcP6FS1EymnzxAxcSJJe/bgc1d7Kk2ZgnvVjDUWUhPgwLewYz5tLh8lTpXjSouXCOo6HMpWLsqPaFWAdoZmv9JCciUJ4QAK2ul68fpFBq0fxJWkK8y9dy6tyjcj+rNFXJk3D+XlRcX/+z/K9X/QvDhO3HlMOxaQtmsJHulXoVITUlsNwaXxQ7h65LGi2m2W7Xblln/qoZbBbDp22WFrfM5EciWJEseZRqwUtKwFqdWcu3qOQesHkZCawILuC6h7yZUzwx4l5dgxfHv2pNLrr+EaGAhnt8KOT9DH16I1bDC2Jqr+AAY+8QTuOa2mlsvnKK4O4ZxqFV3vrMAPe8KlQ7qYSY1BOIX8splmP9aeAaQgZS2oU7GnGLxhMEaTkU87fUDA0g3ELF6Ca2AglSaF4tv5LnNW0x2fwqVDJLmWZUlKF352u58hfe6mT9MqOS+xmYsOMzfm2LwT7OfFtlfvKdRncYb3L2mkxiBKlIKMg7f3kMeiGl1zJPoIQzcMxc3Fjc8qjEY9N46Y8+fxe+QRgoY+heHYMpg9CJJiIKghfzWcxHN7atKn5R18/UB9/Lzd83+TbOzdIWzv9y+tJDAIp2DtDcIRhjwWxc1sX9Q+XvjtBSoYvXnncBPSV72GW0g1Qt4eh0/KZvisPaBJrdWTs3X+S90299PapPm6dRytatx+Qjt7dwjb+/1LKwkMwilYe4NwhCdMW9/MdkbsZOTGkXQ5483AX5NJj1lHQO/2VKh+Cpd/XgbPcuh2w1lfpg/jN17F84ILf7TSrDkQUeAmtezNcNnb+KF454E4wzyUkkjSbgunMLZHPbzcLFf/yukGkV867OJgbVmtsSVsC+N/HMaYVfDcV5dw9zBRo1cqFX1+wIVkeGA2Fwb8wzNhfRj60xVCArxZ9Fxr1hyIYPyKg4THJaG52aS2am94ru+V2QyX9Zwf9oTzUMtggv28UJjb9m3RV2Ktfs2DmdG/sd3ev7SSGoNwCtZOVHOEJ0xbTapbf2Yda+aOYdbvJrzSkqnQNIHy9SJQd/Yw5y66oyunryTwnw/+xNXFhSl9G/JU2+oYXBTPL9ld4Ca13JrhNh27bNeOXkefh1IS2SQwKKV6Au9jXgVkodZ6Zrb9A4BZQObjykda64UZ+54FJmRsn6a1/twWZRKOxRYjhay5QTjKTOfC3sx+2fIZV998l2FnNe6BqVTtkIZHlyehzWAoX4vYhFT8laJmoA8v3lOH/i2CqVzuZq3odprUHKEZTjiGQgcGpZQBmAt0A8KAXUqp1VrrI9kO/VZrPTLbuQHAJKAVoIE9GefGFrZcwnEU90ghZ37C1PGRbJw+kEq/nKEKmoDO7gQ9NwrV4knw8CUhJZ3ZPx/h210X+OXlTlQL8GZE19q3XOd2+jmko1dkskUfQxvglNb6tNY6FVgG9LXy3B7ABq11TEYw2AD0tEGZhAORJUKtcHEfKZ88w+7enamy+gyRIa7UmjeRip/sQ7UfCh6+bDoWRfc5W/hs6xn6NqtCOe/c10TIqZ9DYQ7KHWZuzLGvwZZ9I8K52aIpKRi4kOV1GNA2h+MeUkp1Bk4Ao7XWF3I51zkf9USupIkiF8Z0OPYTpq0fc2X9ES4fLYPJw8CmZ+5k0P99i7ured6ByaQZ9e0+Vu+/SJ2gMnw/rH2+Q1CzNqmFxyWhMFfJIfcam6M0wwn7K67O55+Ab7TWKUqpocDnQIF6s5RSQ4AhACEhIbYvoSgy0kSRTWIM7FkCuxaSdDqKi3sqkBrry7aGishB/+G1HjNxdbn5p+nioqjg68Er3eoy9O5auLvmX9HP2qdjUApjtgwHuXVEO3MznLAdWwSGcKBaltdVudnJDIDWOjrLy4XA21nO7ZLt3M05vYnWej4wH8wpMQpTYFG8HGGkkEOIPAQ7PoGDyzElpRB17k5i90CivxfvP5JM3QeeYELb13BRLvx7+ToTVh5iTI96tKzuz8ReDax+m+x9OtmDQiZramz2Ti8i7MMWgWEXUEcpVRPzjf5x4MmsByilKmutIzJe9gGOZvy8DnhTKeWf8bo7MN4GZRIOpFQ3UZiMcHytOXfR2T/B1YvrnvcR+VsYaZFRnLinFtObneXR5s/xSqtXSDNqPv3jJB9uOoWnqwuXr6UU+C1z6tPJSX41NkdILyLso9CBQWudrpQaifkmbwAWaa0PK6WmALu11quBl5RSfYB0IAYYkHFujFJqKubgAjBFax1T2DIJx1PqmiiSYuGfL2HnAog/D+Wqkd5uPFG/XSL+519xq1mDNWPv4nPXnQxv+gLDmw7nn/OxvPrDQU5GXadXk8qE9m5AkK9ngd/ampqANTU2R0gvIuzDJn0MWuu1wNps20Kz/DyeXGoCWutFwCJblEOI22WzJpOoY+bmogPfQloiVO+I7j6Na2cNRE6diTE+Hr8hz/NWw9P8Hvkno1uOZmCjgQDsOhtLQko6iwa04p47K972Z8mtT8egFCatrf58Mmig9JKZz6LUK3STickEJ9fDjo/h9GYweECTR6DtMNJUEJFvTOH6xo14NmxI0Kcf8X8X5/F3xN+81vY1/NO78tuRS9zXoCLPd6zJf9tVx8ejcH+WufXpFDSVhAwaKL0kMIhSL795FrnWJJLjYe/XsHM+xJ4B3ypwbyi0eBbtFUDc8u+JmjUQnZ5O0NixuD/Rn5F/vMy+y/sY0zyUTTtrsv7IHu6uW4H7GlTEzeCCm6HwU4ts1acjgwZKL1moR5R6NV9dQ25/BV5uhltujB9096XbtR9h39eQeh2qtYW2w6B+bzC4kXr2LBETQ0nctQvvtm2pPOUNkir5MXTDUI7HHOf+Sv9j9bYg0k0mRt9Xl4Eda9okIBQFGZVUsshCPUJYKbcmE6W4ERQUJjq7HGQAv9L19/3g4gaN+psDQnALAHR6OjELF3L5w49Q7u5EDhrNhLRaRMz/G9+ai1DuVxhYdxKzV7nRqY4f0/s1JqR8/msu25M9Bw1IULIfCQyi1MupycTNRZFm0niTzEOGLQwwrKOWSwRR2o85aQ8x+tUZ4Huzgzj56FEiXp9A8pEjlLn3Xvb1H8z/bY4gWV/Cu/oCjIZ4TGEDqNaoFV8970mH2uVzXWLTWW6IRVlOGSprXxIYRKmXU5u8f0o4/dLW8KhhM2VVEvtMtRiV+gJrTO0I8vNldEZQMKWkcGXuPKI/+wyDnx/B772Hb4/uPP7WJpKJwrv6QpQhkaTzz2NMqsE760/kmcLaWW6IRV1OGSprXxIYHJyzPD06u37Ng+nXrAqc+QN2fIrp2C8YDS78YmrD4vSe7NV1bhyb2fmauHs3ERMmknr2LOX696fi/43lpzMJzHprExGJ58xBQaWTeH4wpuSqQP5DPZ3lhljU5cxrqKz8TRQ9CQwOzFmeHp1eaiIcWGaenXz5GHgH8oXrw8xL6EIU/haH+nm50btOOSLeeIO4b5bhFhxMtc8WUqZDB1btDWfs8v0Y3cLxqv4ZoEg8PwRTSqUb5+c31NNZ5g4UdTlz6/cp5+UmfxPFwDGHQgig5KWrXrU3nA4zN1Lz1TW5pn4uVnHnYf1EmF0ffh4NBnfoOw9GH8av1xtccwu0ONzLzcCsqtc43as3ccu+JeDZZ7jjp9WU6dABgNAfD2F0P4939QWgXUk8N9QiKFgz1NMRlia1RlGXM7cU4FkHBGRy5r8JRyU1BgdmzVOZs1SrHab2ozWc22aenXxsDaDMw0zbDoOQduahSNza71DXM50p59dQdvkmDHVqU/X99/Bq2hSAdKOJJX+dJUGdwLvaErSxDInnB6HTbqbGDrbyd+MscweKupy5zcUY/e2+HI93tBqVs5PA4MDym3nqMDdbKxRFm3SBgmJaEhz83txcdOkgePlDh5eh1fPgVy33a3avS9ewf7j05kyMCQkEvjiSwMGDUe7mtRIOhcczfsVBjsTtwivkS0xpfiSdH4ROL2fx9taumWzPhIOr9oYzefVh4pLSAPD3dmNS74Y5vndxlDOnobKZ60tk52g1KmcngcGB5fdU5iwdlWD7Nmmrg2J8OOz+DHYvhqQYCGoIvT+AJo+Cm1ee10y7GM71/33ExcijeDVtSuVpU/Goc7MTevb648zd/C9lA47jW/1LjCkVSDo3EG0sY3Fd/zxWWsuJPeYOZPaPpJluTvWLTUxj7Pf7b5Qp67FZA8Kcx5oVW3mdpUbl7CQwOLD8nsqcpaMSbJ93J8+g2KwKXNhpzl10ZDVoE9T7D7QbBjU63Wguyu2aSpvodeYvBhz5BRdt4ps2D/PG4skog2Wbt4+HKx2aXuBAyuc0KN+A3kGhhJ49TVqWedRuBsWk3g1v6zMWp1nrjlsEhUxpRm3xoGHvWmqpTuFejCQwOLi8nh4dPclZ1ifLcl5uuBkUacabN5/CPOnlFPzcSaPt1T9h/hsQsQ88ykG74dBmMPjXsOqaVa9FMWrvdzSMOcueoLp80PRhLvsEMMVg4Mr1FKb+fIRuDSrSq0kVKlTez76zc2lRsQVz752Lj5sPXq6+TnnTyuthIus+R6illroU7nYggcGJOXK1OvuTZVxSGm4uCn9vN+IS0wp908waFCsQy9Ouv/Ok4TcqqKuQVg8eeBeaPA4eZfK5kplOS2Pwuc302v8LyQZ33m3xGL9VawVK4QKMX3GAXw5FkpCSTpOqfnx99Gtm7pxJhyodmNN1Dl6u5mBclDetohxokNtDRua+TM5USxW3TwKDEyvuanVBbkw5PVmmmTTe7q7sDe1e6LKM7VGPpStW8gRrecBlO66Y+IPmuLd/gQ7dH861uSgnSQcPETFhAg8eP862qs34qFFf4jx9b+w3Ad/svEDNQB++H9aezZe+ZebO97mn2j3MunsW7gb3Qn+eTLl9x0XdhDO2R71b+hjA3BSW9UHD0WupwjYkMDi54qpWF/TGVGRPlumpcHQ1/XZ/TD/DbhLw4qv0bvzq3Ycn7+9SoO/ClJTE5Q8+JObzz3ENDKTqvLkE+tfj2nf7zcNas0lOS+eX8CXMPzCf+2vez/SO03FzKVjHcl7y+o6Lugkn8xr5jUpy5FqqsB0JDMIqBb0x2fzJ8vpl2LPEPMLoWgQE1IL738an2ZMM9PBlYAEvl7B9OxETQ0m7cAG/xx4jaMwrGHx96QeMynGsvCbG43vmH9jGQ3UeYmK7iRhcDDkcd/vy+o6LownHmocM6fwtHSQwCKsU9MZksyfLi/vMcw8OfQ/GVKh9n3m4ae37wKXgE/eNV68SNWsWccu/x616CCGff45P2zYAXEtO450cZ9Ca8Ki0Cnf/naTG3MX6LZ1o5hVp85thXt+xIzXhSOdvySeBQViloDemQj1ZGtPh2E/mgHD+b3DzgRbPQpshUKHubX+Gqxs2cGnKVNJjYig/6HkCR47ExdMTgPWHIwn98TCXriXTqU4gu87EkJxuAox4Vvket3J7SbnSldTL3blIMuNXHGT3uRg2HbtssyfnvL7jktKE4ywz9Us7CQzCKrdzYyrwk2VijLm5aNdCuBoOftWhx5vQ/GnwLJfv6blJv3yZyKnTuLZ+PR7161P1k4/xanhzbsGpqOsM/WoP9Sr68vHTLWge4s+qveG8ve4wsT5f4lr2MClRPUiN7nrjnKQ0I19vP39jxoItOoPz+o5LQhOOvedACOvZZGlPpVRP4H3AACzUWs/Mtv9/wCAgHbgMDNRan8vYZwQOZhx6XmvdJ7/3k6U97aPInvYiD8HOT+HAd5CeDDU7Q9vhULcHFKIdX2tN/IoVhL05E1NSMl/d2Z3trXryv/sb0KdpFXafi6VNTXM+oz9OXOauWuVvLLGZnJ7M6M2j2Rq+lZTIXqTGdrTqPYP9vKxOf5GTkvxE3WHmxhxrRIX9zoT1rF3as9CBQSllAE4A3YAwYBfwhNb6SJZjugI7tNaJSqnhQBet9WMZ+65rra0bbJ5BAkMJYDLC8V/MyezO/gmuXuY0FW2HQcUGhb586oULRISGkvj3dg4H3sGcpg8T7hsEgIerC1XKeXE2JoG1L3WifuWyFucmpCXw4sYX2R25m9D2ocxZEZDrGP/sFHBm5gOFLn9JlNva2vKdFZ/iXPO5DXBKa306442XAX2BG4FBa70py/Hbgadt8L7CGSXFwt6vYOd8c9rrctXgvjegxTPgHZD/+fnQRiMxX3zJ5Q8+QLm48FW7x1hasSVa3eyoTkk3cTYmgVkPN+XOSr4W519Nvcrw34Zz+MphHqo29kZQUGBxU8v+OlNpG89fkBqOI3Wgi7zZYj2GYOBCltdhGdty8zzwS5bXnkqp3Uqp7UqpfjYoj3BEUcfMax7MbgDrJ5gDwqNfwEv7oOMomwSF5OMnOPvEk0S99RY+bdpwx88/sbRSa4ugkElrcHVRFusuxyTHMGjdII5EH+HRkNdYtqnCjRuZxhwMwNz08VS7kBzXC3C2zuDCyOwzCI9LQnOzzyC3dTZyW2OhNH1nzqJYO5+VUk8DrYC7s2yurrUOV0rdAWxUSh3UWv+bw7lDgCEAISEhxVJeUUgmE5xcb24uOr0JDB7Q5BFoMxQqN7Hd26SmEv3JJ1yZvwCDry9V3n2Hsv/5D0opKpXzJCI+OcfzsnZ8RiVGMXj9YMKvh/PhPR8y7stUktIsn241lu3hraoHlNj+AGsUdG5LSehALy1sERjCgWpZXlfN2GZBKXUf8Dpwt9Y6JXO71jo849/TSqnNQHPglsCgtZ4PzAdzH4MNyl2qFWknZ3I87Ftqbi6KOQ2+VeCeidDyOfApb5v3yJC4d6953eV//6Vs795UfG08rv7+aK35+cBFrqek35K8L1PmTaxNHcWg9YOITorm4/s+pnWl1lyMW5Pj+2Wda1Dax/PfzqS70v6dOQtbBIZdQB2lVE3MAeFx4MmsByilmgOfAj211lFZtvsDiVrrFKVUINABeNsGZRJ5mLDqoM2HWgJw5ZR5dNG+pZB6Haq1hXsmQP0+YLBd6ggAU0ICUXPeI/brr3GtXIlq8z+lTOfOgPnzTFx1iI3HomgcXI4eDSvyzvoTOV4nIuECz/46g4S0BOZ3n0/TCuZV2aQ9PH/yHZVchQ4MWut0pdRIYB3m4aqLtNaHlVJTgN1a69XALKAMsDyjTTdzWGp94FOllAlzf8fMrKOZhO2t2htuERQy3XbeHZMJ/t1obi46tQFc3KDRQ9B2KAS3sFm5s7r+51YiJ00iLSIC/6eeosKoURjK+ADw3e4LTF59GK1hwgP1GXBXDVwNLnyz88ItNzEX90uUqfEZKekGFvVYxJ0Bd97YV1ImlBUl+Y5KLpv0MWit1wJrs20LzfLzfbmc9xfQ2BZlELnL2mzkolSOo2mggHl3Uq7D/m/Ms5OjT4JPEHQZb24u8q1ok3Jnlx4bS9TMt4j/8Ufc77iD6l9/jXeL5hbHeLsbaFMzgKl9G1EtwPvG9uw3MRfPMLxDFlHGw5MlPRdzh98dFtexdXt4SZyfIH0GJZdNJrgVN5nHYL3ss03zYtVEo5gz5pnJ/3wJKfFQpbl5MlrDB8HVdumns9Jac3XtWi5NfxPj1auUHzyIwOHDcXF3JynVyHu/n6BCGQ8GdbqDzP+fVQ5ptzNvzpEpR/EOWYKfR1mW9lpCtbLVbjnWlnL6HXi5GZjRv7HcREWxKs55DMKB5TRyJCcK8m8C2P4x/DrePBu5QV9zQKjaqkBrHxRUWmQkkW9M4fqmTXg2bkzI4kV41jOX88+Tl3l95SHOxyTybPvq5s+RR1n6NQ+mcqUwXtz4OUHeFVnQbQGVy1QusrJnym30zivf3bqeshCOQAJDCWdN85ACnmoXkv8Nqvpd0HkMtBoIZavYpoC50CYTcd99R9Ssd9BGI0HjxhHwzH9RBgMxCalM+/kIK/aGc0egD98Mbkf7WvmPdtoStoXRm0YTUjaEBd0XEOgVWKSfIVNuvwOj1pIrSDgkCQwlXG4jRwxKYdK6YO3ClZua/ytCq/aG88V3W3h8y1c0jj5NQsPmNJkzE/csc1fORifw88EIXrynNiO61sbTLf98SuvPrmfclnHUDajLp/d9ip+nX67vb+s287yWzSzu9ZKFsIYEhhIut5Ejjti+vWrXOXa+/SGTDq8j1cWVOc0e4c/a7ZkRbaC5TwJbTlzmv+1r0CLEn23j7qGCr4dV113972ombptI0wpNmXvvXHzdfW85ZtXecIvVy8B2w3hz+h1kZW2nf0nswBaOSQJDCecsI0eSjxzB8+VR/DfmAtsqN2Ju0/7EepaFdBOhPx4iJd2Eh6sLvZpUwd/H3eqg8O2xb5m2YxptK7flg64f4O3mfcsxeXXQ2+KJPvPcV77bjzGHwR7WjPuXlNWiOElgKAUcebapKTmZK3PnEb1oEb6u3kxr/Qzbgi3TZVxNTqdnw0pM7tMQfx/rRz59fvhz3tn9Dp2rdmZ2l9l4GHIOJvl10Nti+czM7/92x/0X9ZrPQmQlgUHYTeKuXeZ0FufOUe6h/rzk3p6TSbeOKgrwdueT/7a0+rpaaz7Z/wnz9s+je/XuzOw0E7c8Zl7nd+O31UzewtTeimPNZyEySWAQxc54/TpR77xD3LJvcatWjZDFi/Bp354Re8P5v+8PkGo03TjW09WF0N7Wr8+gtWbOnjksPryYPrX68MZdb+Dqkvf/5nl1Dtt6Ju/t1t4k/YQoTrZIuy2E1a5t3MTpB3oR991yAgYM4I4fV+HTvj1R15LZcOQSqUYTAT7uKMwT7mY+1MTqG6lJm5i+YzqLDy/msXqPMbXD1HyDAuScDhrA39vNYTrpJWW1KE6lrsYgIzvsIz06mkvT3+Tq2rV41KlD1Y8+xKtxY0wmzTc7z/Pm2qOkpJsY070uQzrXwt21YM8sRpORSX9N4sd/f2RAwwH8r+X/8pzslpUzdNA7QxlFyVGqUmJIaoLip7Um/scfiZoxE1NiIuWHDyNw0CCUu7kTecTSf1hzIIJ2dwTw5oONuaNCgVZ5BSDNmMb4reNZd3YdLzR9gWFNh1kdFKwhDxOipJCUGDmQkR3FKy08nIhJk0nYuhWvZs2oPH0aHrVqkZJuxMVows3gQt+mVbi7TgUeaVX1tm7mKcYUXtn8Cn+E/cErLV9hQKMBNv0MMkxUlEalKjDIyI7ioY1G/nj7Y8p9vQCtYWXbR2j50iBq1Aph9vrjzNv8L+kmTXDG0/ejrW8viV1iWiIvb3qZ7RHbeb3t6zx+5+M2/iTyMCFKp1IVGGRkR9FLOXWKw6PHUfHkEXYH1ePDZg8R5R3A8pWHmfvHGU5GXb9xbGGevq+lXmPk7yPZd3kf0zpMo2/tvjb9HJnkYUKURqVqVJKM7DBbtTecDjM3UvPVNXSYuTHXxdsLQqemcnnePM482B/j+XPMavEEE9sPIso7AIDkdJNFUMiU+fRdEHHJcQxeP5gDlw/wVue3iiwoQO4PDfIwIUqyUlVjkJEdRdNmnnTgABGvTyDl5EnK/ud+HtdtiPO4NR9Rbgry9H0l6QqD1w/m/NXzvNf1Pe6udvftFNlqskqZKI1KVWAAx04PURxs2WZuSkzk8gcfEvPFF7hWqEDVefPwvacr3jN+Jy4++ZbjDUrddq4ggMiESAavH8ylxEvMvW8u7Sq3K1B5b4c8TIjSqNQFhtLOVm3mCX//TcTEUNLCwvB7/DGCXnkFg68vxyKv4ma4tYXSy83AQy2D+WFP+G09fV+4eoHBGwYTnxLPp90+pXlQ83zPsZXS/jAhSh8JDKVMYTvgjfHxXHr7beJ/WIF79epU//ILvFu3JjnNyLu/HmP+ltOU9XLj6bYhbDwWRUR8ssVTdqvqAbk+fec2X+B03GkGrx9MiimFhd0X0jCwoU2/E0cl8yeEvUhgKGUK02Z+/Y8/uDhhAsaYWMoPHkzgiBdw8fQE4PK1FBZvO0u/5sG8/p/6uWZBze3pO7e+j4ik0yy7MAGFYlGPRdT1r3s7H9vpyPwJYU82GZWklOqplDqulDqllHo1h/0eSqlvM/bvUErVyLJvfMb240qpHrYoj8hdv+bBzOjfmGA/rxv5iKyd+a1NJlwrVKDm8u8IeuV/xBtdWPjnabTWVAvwZtOYLrzzSNMCpcbOlFPfR4rhDPOOv4KbixtLei4pNUEB8u4LEqKoFbrGoJQyAHOBbkAYsEsptVprfSTLYc8DsVrr2kqpx4G3gMeUUg2Ax4GGQBXgN6VUXa11/qvXi9t2u23mvl27UqZzZ3BxYdXecKb8fISrSWl0rBPInZXKUqmc522XKXsfh8H7NF5Vl2BKL8MX939BlTJFu8a0o5H5E8KebFFjaAOc0lqf1lqnAsuA7APL+wKfZ/z8PXCvMuc/6Ass01qnaK3PAKcyricc1IW4FJ5ZtJNR3+4jJMCbn17syJ2Vyhb6uln7OAw+J/CqthhTuh9lY0eVuqAAMn9C2JctAkMwcCHL67CMbTkeo7VOB+KB8laeKxyE0aT576Id/HMulsm9G/DD8LuoX7nwQQFuTj40lDmCV7XPMaUGosOHM65b6XxOkMmYwp6cpvNZKTUEGAIQEhJi59KUTgYXxTuPNCXYz8vmT66ZTVszf79EfEIt/K8/x//1a15qO1pl/oSwJ1sEhnAgaxa0qhnbcjomTCnlCpQDoq08FwCt9XxgPpjTbtug3OI2tK4RcFvnWTP00tz38SzwrA1K6vxk/oSwF1s0Je0C6iilaiql3DF3Jq/Odsxqbv61Pwxs1OaFIFYDj2eMWqoJ1AF22qBMwoFkDr0Mj0tCc3PopS1yNAkhbK/QNQatdbpSaiSwDjAAi7TWh/IyBjMAAAi2SURBVJVSU4DdWuvVwGfAl0qpU0AM5uBBxnHfAUeAdGCEjEgqeQqShkMmdQlhfzbpY9BarwXWZtsWmuXnZOCRXM6dDky3RTmEY7J26KVM6hLCMZSqtNvCPqwdeimTuoRwDBIYRJGzduilTOoSwjFIYBBFzto0HDKpSwjH4DTzGIRzs2bopSyKI4RjkMAgHIZM6hLCMUhgEA5FJnUJYX/SxyCEEMKCBAYhhBAWJDAIIYSwIIFBCCGEBQkMQgghLEhgEEIIYUECgxBCCAsSGIQQQliQwCCEEMKCBAYhhBAWJDAIIYSwIIFBCCGEBQkMQgghLEhgEEIIYUECgxBCCAuFCgxKqQCl1Aal1MmMf/1zOKaZUupvpdRhpdQBpdRjWfYtUUqdUUrty/ivWWHKI4QQovAKW2N4Ffhda10H+D3jdXaJwDNa64ZAT+A9pZRflv1jtdbNMv7bV8jyCCGEKKTCBoa+wOcZP38O9Mt+gNb6hNb6ZMbPF4EooEIh31cIIUQRKWxgqKi1jsj4ORKomNfBSqk2gDvwb5bN0zOamOYopTwKWR4hhBCFlO+az0qp34BKOex6PesLrbVWSuk8rlMZ+BJ4Vmttytg8HnNAcQfmA+OAKbmcPwQYAhASEpJfsUu1VXvDmbXuOBfjkqji58XYHvVkHWUhhNXyDQxa6/ty26eUuqSUqqy1jsi48UflclxZYA3wutZ6e5ZrZ9Y2UpRSi4ExeZRjPubgQatWrXINQKXdqr3hjF9xkKQ0IwDhcUmMX3EQoNDBQQKOEKVDYZuSVgPPZvz8LPBj9gOUUu7ASuALrfX32fZVzvhXYe6fOFTI8pR6s9YdvxEUMiWlGZm17nihrpsZcMLjktDcDDir9oYX6rpCCMdT2MAwE+imlDoJ3JfxGqVUK6XUwoxjHgU6AwNyGJb6tVLqIHAQCASmFbI8pd7FuKQCbbdWUQUcIYTjybcpKS9a62jg3hy27wYGZfz8FfBVLuffU5j3F7eq4udFeA5BoIqfV6GuW1QBRwjheGTmcwkztkc9vNwMFtu83AyM7VGvUNfNLbAUNuAIIRyPBIYSpl/zYGb0b0ywnxcKCPbzYkb/xoXuJC6qgCOEcDyFakoSjqlf82CbjxbKvJ6MShKi5JPAIKxWFAFHCOF4pClJCCGEBQkMQgghLEhgEEIIYUECgxBCCAsSGIQQQliQwCCEEMKCBAYhhBAWlNbOl8FaKXUZOGejywUCV2x0LVuSchWMlKtgpFwFU1LKVV1rne8Kmk4ZGGxJKbVba93K3uXITspVMFKugpFyFUxpK5c0JQkhhLAggUEIIYQFCQwZy4U6IClXwUi5CkbKVTClqlylvo9BCCGEJakxCCGEsFDqAoNSKkAptUEpdTLjX/88ji2rlApTSn3kCOVSSjVTSv2tlDqslDqglHqsCMvTUyl1XCl1Sin1ag77PZRS32bs36GUqlFUZSlguf6nlDqS8f38rpSq7gjlynLcQ0oprZQqlhEu1pRLKfVoxnd2WP1/e2cbmlUZxvHf34ZIpFlGulIpYULL3mxISmWSxFygfTQaNBiBLoI+9EHYl6hPGgZFQsSIlmCvnwQVYqIMtGmC2qQP86WgpTkoCyKqRVcfzi2dM55nz1l7zv0ctusHB65zzs1zfrvPvXPdL+fZpH1l8JK0XNIRSafDveyI4PS+pDFJ56qcl6S3g/PXklYX7ZTT67ngMyzpuKQHpn1RM5tVG7AL2BHiHcDOScq+BewD3imDF7ASaAnxHcAVYGEBLjcAF4EVwFzgLNA6oUwP8G6ItwKfRKijPF4bgBtDvL0sXqHcfGAQGALayuAFtACngVvC/u0l8XoP2B7iVuC7CF6PA6uBc1XOdwCHAAGPACeKdsrptS51/zbVw2vWjRiALUB/iPuBZyoVkvQwsBj4oixeZjZiZudDfBkYA2p+WeV/sAa4YGaXzOwv4OPgV833c+BJSSrAZUpeZnbEzH4Pu0PA0oKdcnkFXgd2An9EcMrr9QKwx8yuAZjZWEm8DFgQ4puBy0VLmdkg8PMkRbYAH1rCELBQUnOjvczs+PX7R53a/GxMDIvN7EqIfyR5+GeQNAfYDbxSJq80ktaQ9LYuFuByJ/B9an80HKtYxsz+Bn4FFhXgMlWvNN0kPbyiqekVph2WmdmBCD65vUhGoSslHZM0JKm9JF6vAp2SRoGDwEsRvGox1fbXCOrS5mfkv/aUNAAsqXCqN71jZiap0mtZPcBBMxutZye4Dl7XP6cZ2As8b2b/1E1wBiGpE2gD1pfAZQ7wJtDVYJVKNJFMJz1B0tMclHSfmf3SUCt4FvjAzHZLWgvslbTK23t1JG0gSQyPTvezZmRiMLON1c5Juiqp2cyuhAdspaHzWuAxST3ATcBcSb+ZWdVFxUheSFoAHAB6w3C2CH4AlqX2l4ZjlcqMSmoiGe7/VJDPVLyQtJEk2a43sz8LdsrjNR9YBRwNHY0lwH5Jm83sVAO9IOn1njCzceBbSSMkieKrBnt1A+0AZvalpHkkfxcoxlRXNXK1v0Yg6X6gD9hkZtP/PYyxeFKmDXiD7CLvrhrlu4iz+FzTi2Tq6DDwcsEuTcAl4G7+Wxy8d0KZF8kuPn8aoY7yeD1EMr3WErFN1fSaUP4ocRaf89RXO9Af4ttIpkoWlcDrENAV4ntI1hgUoc7uovoi79NkF59PRmxjk3ktBy4A6+p2vVg/WFk2knnww8B5YAC4NRxvA/oqlI+VGGp6AZ3AOHAmtT1YkE8HMBIesr3h2GvA5hDPAz4LDfIksCLS/avlNQBcTdXP/jJ4TSgbJTHkrC+RTHN9AwwDW0vi1QocC0njDPBUBKePSN70GycZSXUD24BtqbraE5yHI97DWl59wLVUmz813Wv6N58dx3GcDLPxrSTHcRxnEjwxOI7jOBk8MTiO4zgZPDE4juM4GTwxOI7jOBk8MTiO4zgZPDE4juM4GTwxOI7jOBn+BRgNQ+XJrdMgAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10724bdd8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# If the noise on input and output is the same, least rectangle or orthogonal \n",
"# regression works best\n",
"n = 100\n",
"x = np.linspace(0, 1, n)\n",
"y = x\n",
"\n",
"x_o = x + np.random.normal(0, 0.2, n)\n",
"y_o = y + np.random.normal(0, 0.2, n)\n",
"\n",
"plt.plot(x, y, '--')\n",
"plt.scatter(x_o, y_o)\n",
"\n",
"lin_x = LinearRegression(fit_intercept=True)\n",
"lin_x.fit(x_o[:, None], y_o)\n",
"plt.plot(x, lin_x.predict(x[:, None]), label='y=f(x)')\n",
"\n",
"lin_y = LinearRegression(fit_intercept=True)\n",
"lin_y.fit(y_o[:, None], x_o)\n",
"plt.plot(lin_y.predict(y[:, None]), y, label='x=f(y)')\n",
"\n",
"w, b = odr_fit(x_o, y_o)\n",
"plt.plot(x, x*w + b, label='odr')\n",
"\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1073fa908>"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VMX6wPHvZNNDQnoCCZDQew0gvSlgQ8Du9Sr2exURUBQsKAiCF6WIV6/1B/YuoihYKFIEpUrvNSEVEtKTzc7vjw0hQMom2Zq8n+fhIbt79pw5ye57Zt6ZM6O01gghhKhb3BxdACGEEPYnwV8IIeogCf5CCFEHSfAXQog6SIK/EELUQRL8hRCiDpLgL4QQdZAEfyGEqIMk+AshRB3k7ugClCc0NFTHxMQ4uhhCCOFStmzZkqq1DqtsO6cN/jExMWzevNnRxRBCCJeilDpuyXaS9hFCiDpIgr8QQtRBEvyFEKIOctqcf1kKCws5deoUeXl5ji6Kw3l7exMdHY2Hh4ejiyKEcEEuFfxPnTqFv78/MTExKKUcXRyH0VqTlpbGqVOniI2NdXRxhBAuyKXSPnl5eYSEhNTpwA+glCIkJERaQEKIanOpmj9Q5wP/efJ7EKL2WLItnjkr9pOQnkvDQB8mDWvFyC5RNj2mywV/IYSoTZZsi2fKNzvJLSwCID49lynf7ASw6QXApdI+rmzt2rW0a9eOzp07k5uby+nTp7nuuusqfM8PP/zA1KlT7VRCIYQjzFmxvyTwn5dbWMScFfttelwJ/nby8ccfM2XKFLZv346Pjw9z587lgQceqPA91157Ld9//z05OTl2KqUQwt4S0nOr9Ly1SPCvoqlTpzJ//vySx8888wwLFiyo8D3vvvsuX3zxBc899xz/+Mc/APj6668ZPnw4APPmzePee+8FYOfOnbRv356cnByUUgwcOJAffvjBRmcjhHC0hoE+VXreWlw35//TZEjcad19RnaAq2dXuMm9997L6NGjGT9+PCaTic8++4yVK1fSuXPnMrf/5JNPuP/++1m3bh3XXXcdN910E0ePHiUoKAgvLy8AHnvsMQYOHMi3337LzJkzeeutt/D19QUgLi6OtWvXcsstt1j3XIUQTmHSsFYX5fwBfDwMTBrWyqbHdd3g7yAxMTGEhISwbds2kpKS6NKlC02aNGH79u0W7+P06dOEhV2YdM/NzY1FixbRsWNHHnroIfr06VPyWnh4OAkJCVY9ByGE8zjfqSujfSxVSQ3dlu6//34WLVpEYmIi9957L5mZmfTr16/MbT/55BPatm170XM+Pj6XjdE/ePAg9erVuyzQ5+Xl4eNj2+afEMKxRnaJsnmwv5TrBn8HGjVqFFOnTqWwsJBPPvkEg8FQpZp/y5YtOXbsWMnjjIwMxo0bx++//87YsWP56quvuOmmmwA4cOAA7du3t/YpCCHqOOnwrQZPT08GDRrELbfcgsFgqPL7/fz8aNasGYcOHQJgwoQJPPLII7Rs2ZL33nuPyZMnk5ycDMCqVau49tprrVp+IYSQmn81mEwmNm7cyJdffmnxexYtWnTR47Fjx7Jo0SJmzJjB+++/X/J8o0aNSi4KSUlJ5Obm0qFDB6uUWwghzpOafxXt2bOH5s2bM2TIEFq0aFHt/YwaNYrKlqk8ceIEr776arWPIYQQ5ZGafxW1bduWI0eOWGVf999/f4Wvd+/e3SrHEUKIS0nNXwgh6iAJ/kIIUQdJ8BdCiDpIgr8QQtRBEvztRKZ0FkI4Ewn+diJTOgshnIkE/yr666+/6NixI3l5eWRnZ9OuXTt27dpV4Xsqm9K5f//+F00P0bdvX3bs2CFTOgshbMZlx/m//OfL7Duzz6r7bB3cmqd6PFXhNt27d2fEiBE8++yz5Obmcuedd9KkSZMaTel83333sWjRIubPn8+BAwfIy8ujU6dOgEzpLISwDZcN/o40depUunfvjre3N6+99lqVJ3a7dErnm2++mRdffJE5c+bw/vvvM2bMmJLXZEpnIYQtuGzwr6yGbktpaWlkZWVRWFhIXl4eJpOpRlM6+/r6ctVVV/Hdd9/xxRdfsGXLlpLXZEpnIYQtuGzwd6SHHnqIF198kaNHj/LUU0/x+uuv12hKZzBP9XD99dfTr18/goKCSp6XKZ2FELYgHb5V9MEHH+Dh4cEdd9zB5MmT+euvv1i5cmWV9nHplM4A3bp1IyAggHvuueeibWVKZyGELVgl+Cul3ldKJSulyhz2opQaqJTKUEptL/7nsoPX77rrLr7++msADAYDmzZtYvDgwZW+b9GiRSULtMCFKZ3PS0hIwGQyMXTo0JLnZEpnIYStWKvmvwgYXsk2a7XWnYv/TbfScV1W6SmdP/jgA3r27MnMmTNxc7vwJ5EpnYUQtmKVnL/W+nelVIw19lWXnJ/S+a677uKuu+667HWZ0lkIYSv2zPn3UkrtUEr9pJRqV9YGSqkHlVKblVKbU1JS7Fg0IYQAbTQ6ugh2Y6/gvxVoorXuBCwElpS1kdb6ba11nNY6rvQ4eCGEsCVjairxTz5J/ONPOLoodmOXoZ5a63Olfv5RKfWGUipUa51qj+MLIcSllmyL55Wf9tJ5+0ru2bscH1MhoQ/ejzaZUG61fyCkXYK/UioSSNJaa6VUD8wtjjR7HFsIIS61ZFs87727jElbvqRl+im2hbXgvW438Wj/IYysA4EfrBT8lVKfAgOBUKXUKeB5wANAa/0/4Cbg30opI5AL3Ka11tY4tjN64YUXqFevHk88UXeakEK4iqKMDBJemMbs/etI96rH7Lh/sCaqMyjFnBX7GdklytFFtAtrjfa5vZLXXwdet8axXJnRaMTdXW6qFsIRtNacW7qUpP/MoX/aGZY27cOHbYaR43Fh+pSE9FwHltC+6kb7xsrmzp1L+/btad++PfPnzwdg5syZtGzZkr59+7J///6SbQcOHMj48eOJi4tjwYIFjiqyEC5nybZ4+sxeSezkZfSZvZIl2+Krva/8gwc58c+7SHhqMp7R0cy4/ine6jjyosAP0DCw7syj5bLV0MSXXiJ/r3WndPZq05rIp5+ucJstW7bwf//3f2zatAmtNT179qRfv3589tlnbN++HaPRSNeuXenWrVvJewoKCti8ebNVyypEbbZkWzxTvtlJbmERAPHpuUz5ZidASVpmybZ45qzYT0J6Lg0DfZg0rNVlKRtTTg6pb7xB2qLFuPn5ETl9GoE33cRtO06zo9T+AXw8DEwa1spOZ+h4Lhv8HWXdunWMGjUKPz8/AEaPHs2yZcsYNWoUvr6+AIwYMeKi99x66612L6cQrmzOiv0XBWaA3MKikpx8ZRcHrTVZv/1G4syXMJ4+Tf3Rowl/4nHcg4NLtjl/nIouHrWZywb/ymrozuT8hUIIYZnycu/nn6/o4nBNqImkGTPJWrMGrxYtiPr4I3xLtcTPG9klqk4F+0tJzr+K+vXrx5IlS8jJySE7O5tvv/2Wa6+9liVLlpCbm0tmZibff/+9o4sphEsrL/d+/vmyLg4eRUb6bfqeI9ddT85ffxH+5JPEfvN1mYFfuHDN31G6du3KmDFj6NGjB2Cen6dbt27ceuutdOrUifDwcJmTR4gamjSs1UVpHbg4J98w0If4UheAzskHeOTvb4nOSqHe8OFETH4Kj8hIu5fblShnHW4fFxenL+0k3bt3L23atHFQiZyP/D5EbVZRh+75nL/PuTM8sOt7BsZvJ6FeKEVjH2fomJEOLrljKaW2aK3jKttOav5CCKdUUU7+hg4RBP70DQFL3sdQZOS7ztfQZsJYbugZa+dSui4J/kIIl5K7fTunp00nYu9e/Pr2JfK5Z+nUpImji+VyXC74a61RSjm6GA7nrOk6IWzFePYsKXPnkf7ll7hHRBA1fz7+w4ZKPKgmlwr+3t7epKWlERISUqf/4Fpr0tLS8Pb2dnRRhLA5bTKR8e0Skl95haJz5wi+5x5CH3kEQz0ZQl0TLhX8o6OjOXXqFLLQi/lCGB0d7ehiCGFTefv3kzhtOrlbt+LTtSuRzz+Pd6uWji5WreBSwd/Dw4PYWOnQEQIsm97AVRVlZZP6+uuc+fBDDAEBNJg5k/qjRtaJefbtxaWCvxDCzJK5b1yR1prMFStIemkWxuRkAm++mbCJE3APCnJ00WodCf5CuKDK5r5xRQXHj5P44gyy163Dq00bol9bgE/nzo4uVq0lwV8IF1TZ3DeuxJSfT9o775L29tsoDw8inn6aoDtuR8naFzYlv10hXNCl0xuUft6VZK1dR+KLL1J44gQB11xD+OSn8AgPd3Sx6gQJ/kK4oMrmvqkOe3YgFyYmkjRrNpkrVuAZE0Pj99/Dr3dvmxxLlE2CvxAuyNrz0VenA7k6FwtdWMiZjz4mdeFCdFERYY+NI/i++3Dz9KxWuUX1SfAXwkVZcz76qnYgV+dikbN1K4kvTCP/wAHqDRhAxHPP4in3qjiMBH8hRJU7kKtysTCePUvynFfI+OYb3Bs0IPr1hdQbMqRO36XvDCT4CyGq3IFsycVCm0ykf/UVKa/OpSg7m5D77yP04YdxK17uVDiW3C4nhGDSsFb4eBgueq6iDuTKVtrK27OHY7ffTuLU5/Fq0YKm335D+BNPSOB3IhL8hRCM7BLFrNEdiAr0QQFRgT7MGt2h3Px9eReLp/pFkzjzJY7edDOFp+JpMHsWjT/8AK8WLexwFqIqJO0jhACq1oF82Wij+t5MDzhN48kvcTY1lcDbbiV8/HgM9evbssiiBiT4CyGq5fzFIv/IURJfnE7OHxs5FNKYef0fJdurFZOOZDGyi3MG/9o8KZ6lJPgLIarFlJdH6ltvcebd9zB6ePJelxv5rnFPTMoNnHiiudo6KV5VSfAXQlRZ5urVJL04g8L4eAJGXM/9Pr3Zm+9x0TbOOtFcbZwUrzok+AshSlSWDilMSCDxpZfI+vU3PJs1o/Hixfj17MG+ycvK3J8zTjRXmybFqwkJ/kIIoOJ0yA3twjjzwQek/PcN0JqwiRMJGXM3qnhaBleaaM6VympLMtRTCAGUnw75bvEPHBk9muRXXsWvVy+aLfuB0AcfKAn8UPX7BBzJlcpqS1ap+Sul3geuA5K11u3LeF0BC4BrgBxgjNZ6qzWObUsyIqD2kb9p+S5Ne9TPz+T+XT9w5ckt6Kgoot94A//Bg8p8r7UnmrMlVyqrLSmtdc13olR/IAv4oJzgfw3wKObg3xNYoLXuWdE+4+Li9ObNm2tctuq6tAkM5tpBRTe+COcmf9OK9Zm9kvj0XNy0iauPbeTuPT/hbSzg5w5XMnHxLNx86lZaxFUppbZoreMq284qaR+t9e/AmQo2uQHzhUFrrTcCgUqpBtY4tq1UNCJA2MaSbfH0mb2S2MnL6DN7JUu2xVt1/7X9b3ri3Al2pOyo9vsnDWtFu6wE5q5ZyNgd33C4fkMmDp1E8ymT6kTgP50hHb62EAWcLPX4VPFzp0tvpJR6EHgQoHHjxnYqWtlkRIB1VZZuscfYa0v/pq6YGtp0ehMTV08kxCeEb0d8i8HNUPmbii3ZFs/rS7cxdOO3/OfoH2R4+/Ny3D842K4Xk4a3dqpzt+bfRmtNek4hQX6eZOYVMuiV1fwxeQhBfnVjbQGnGu2jtX4beBvMaR9HlkVGBFiPJYHdHmOvLfmbllfWzcfPsGpfilNeED7b9xmz/5xNTEAMC4csrFrg33qKn+YvZtqO76ifn8XSpn34quM1PHdrD6c5v/OsUUEoMmm2HD/Lit2JLN+VSMNAb778V2/8vT2Yf2sXDIa6M820vUb7xAONSj2OLn7OacmIAOuxJN1ij5aWJX/T8sr68cYTxKfnorkQdKydlqqqQlMhMzbOYOammfSN6stH13xEI/9Glb+xWP7hw+iJDzP+z49I9g3isYHjeavjSNLwdMpUWE3Tdu+tO0rPl37llrf+4MM/jtMq0p9b4i78voa3jyTA26OCPdQu9qr5LwXGKqU+w9zhm6G1Pl3JexxKRgRYjyWB3R4tLUv+puWV9dJmaFVbJdZOJaXnpfP4msf5M/FP7ml/D491ecziGr8pJ4fUN/9H2qJFROHOa51uZHlMT7S6UBd0xvRmVSoIOQVGfj+QwordSUy9ri1Bfp74eRro2TSE4e0iGdgqDP86FOjLYq2hnp8CA4FQpdQp4HnAA0Br/T/gR8wjfQ5hHup5jzWOa2vWXCavLrMksNtiQfKyVPY3La+sZbE0QFq7P+Nw+mEeXfkoidmJvNT3Ja5vdr3F781cuZKkGTMpTEig/ujRjPfswf68y8OAM6Y3K/scZeUb+Xl3Iit2J7LmQAp5hSYCfT24o2djuvsFc1uPxtzWo5K+RGMBuEvO32Ja69sreV0Dj1jjWML1WBLYnaWlVVZZFZfX/MHyAFleumL859uZs2J/lc7z91O/8+TvT+Jt8Ob9Ye/TObyzRe8rOHWKpBkzyVq9Gq8WLWjy8Uf4duvGv8sZ/uqM6c2y/jbe7m6M6R0DQGpmPhO/2EFkgDe3xjViWLtIesQG426oILudnwXHN8CR1eZ/9cLgru9seh7Owqk6fEXtZGlgd4aWVlllHdQ6jK+3xFc7QFbUQrC0FaC1ZvHuxczdMpfWwa15bfBrRPpFVnpsU0EBZ95/n9Q3/wcGA+FPPknwP+9EeXiUe77Omt48X6ZZP+4lKTMfT4Mb+UYT206eBSAm1I9l4/rSJjIAN7dyOm6LCiF+64Vgf+pPMBnB4AWNr4AWQ+1zMk7AKjd52YKjb/IStVtVc/A1ydmfv3mqIlGBPqyfPLjM1wqKCpj2xzSWHl7KVU2uYkafGfh6VL4cYvbGjSROm07B0aP4Dx1KxNNT8Iis/ILhzB7+eAs/7kwEoG2DAIa1i+TqDpG0jPAv+w1aQ+oBc6A/vAqOrYOCTEBBg47QdBA0HWgO/B7Ol+qqDktv8pKav6hzqpODr0mrpKx0xaXKax2k5qYyftV4dqTs4OFOD/OvTv/CPFtK+QqTk0l++T+cW7YMj0aNaPT2W9Tr379aZXcUrTU7TmWwfFciGw6n8tW/euPp7kb/FmF0aRTEsHaRNA4p5wJ4LgGOrIGja8xBP7N4bElQLHS40RzwY/uDb7DdzscZSfAXdY6953MvnVoprwVQVv/BvjP7eHTlo6TnpfPKgFcYFjOswuNoo5Gzn35GyoIF6Px8Qh95hJAH7sfN27vmJ2EnR1KyWLzhGD/vSeJ0Rh7ubopezUJIy85n05EzLFx5iIT0XBZtOHah9ZV3Do6vN9fsj6yG1OKhn74hEDvAXLNvOgCCYhx3Yk5Igr+ocxxx9/b5lkN58wtd2n/wy/FfeGbdMwR4BrD46sW0DWlb4f5zd+zg9LRp5O/Zi1/v3kROfQ7PmBhbnIpV5RUWsf5QKo2CfWkZ4c/ZnAI+++skA1qG8cTQVlzZJoL6vh4X/d48MNIwYxvx335E2tojhJz9G3QRuPtAk17Q5R/m2n1Ee3Cr2q1Mrnh3d3VJ8Bd1TmVDBm0ZACrrYNVa87+//8cb29+gY1hHFgxaQKhPaLn7K0pPJ3nuPNK//BL3sDCi5s/Df9iwSlNDjpSVb2TVvmRW7E5k1b5ksguKuK9vLM9d15YujYLYNvUqfD1LhSat+fqnn7nDtJk+Hrvo6bYXP5VPkVbsO9uCkL4TzDX7Rj3B3ava5apryztKh6+ocyqa3RNw2MyfucZcnlv/HCuOrWBEsxFM7TUVL0PZwUybTGQs+Y7kOXMoOneO4DvvJPTRRzHU87NpGaurwGjC090NrTW9Zq0k8VweofU8uaptJMPbR9KraQie7qVq6RmnLozIObIGspMBOGxqwHpTe9aZ2rPR1JZM/Dg6+1qrlLG8jvmKOuOdkXT4CqfkDM3qimrffWavdMj6ronZiTy26jH2pu1lYreJjGk3ptzae97+AyROn07uli34dOlC5PNT8W7d2mZlq67TGbn8vDuJ5bsSSTyXx8rHB6CU4ulr2xAZ4E23JkEYzg/JzD0LB9fBkdVk7vkV/+xjAJyhPtnRffkgP5YfslpxmpCLjhFlxZvR6tpkjhL8hd04U7O6vNE7jggAO1J2MH7VeHIKc1g4eCEDGg0ALr9QPjmgMVes/ZYzixdj8PenwcwZ1B81ClXFvLat/bY3iddWHmLHyXQAWoTX49oODSgoMuHlbmBEp4ZgzIfjay/U7hO2gTZhNPiwzdiKNcY7WWdqz37dCJ/j7tzYLYr0LfFgw5vR6tpkjhL8hd1UdZSNI1oJ9g4A3x/+nhc2vECYbxjvXPUOzYOaA5dcKLUmZvcmgj9/hjO56fzeojf/bTEM/4NBTNpx2qH5aK01e09nsnx3IqO6RBEb6keB0YTWmknDWjGsXSTNw+uByQRJu4qD/So4/gcYc0EZIDoO+k+C2AEM+TSL49nGi46RW1jEqn0pzBrdwaafB3tNMeIsJPgLu6lKrdpRrQR7BQCTNvHa1td4b9d7xEXEMXfgXIK8g0peP3+hbJCdyr93fEv35P0cCWjAS93vZG9wDADnHNRyMpk0206eZfmuRJbvTuTkmVzcFDQO9iUmxJfh7SO5ukMDOHscjnwNa1abx9znpJl3ENYGut1tHpHTpDd4B5Ts+0TGsjKPmZCea/M7wF3pbmdrkOAvbO58Db68oQVl1artPRb/vOoGgLJaKeXt5/PNB5m9ZSpG7114ZPfmuvCnLwr8AClp57jj4CpuPbASo5uBt9qPYGnTPpgumbnTHr8TMHfYpmTlExXoQ1aBkdve3ohC0ad5CGMHNefKNhGEuGXDnu8upHLOHjW/2b+BedqE82PuA8pfxM/RqRdnmGLEXiT4C5sqa2RNaeXVqh3Z+VbVAFBWK2XSVztAQ6FJlzw35ZudrDm6j19SZ4NXMvmJI8g824tnT+/FoNxLjpm1bj1vr55LRGYKa6I6807760nzqV/u8W31O8ktKOL3gyms2JXIr3uTaBHhz9f/7k2AtweL7+1Bh3BP/JO3wJH34OPVcHoHoMHTH2L6whX/Ngf70JZg4dDTupZ6cSQJ/sKmyqrBnxdVQa3a0TXAqijrHAuLLm/nFHgc5JezH4G7JvfEPRTltAAu1N6vbehO0uzZZP60nKDIKKZ1+Rcbg5uXvL+ms4tWxWu/HeTN1YfJLSyivo8HV7WN5Oq2YSWTovU+shpObISifHDzgEY9YNDT5tp9VDcwVC+01LXUiyNJ8Bc2VV6tVEGFY6ddqQZoSc3bI/BPvCKXYCoIIffk3ejCCzduuZmK6LF5BUe+mIQuKiLssXEE33cft+1O4aQVZxctT2pWPr/sSWLF7kRevbkTIfW8aBzsy41dGzIqpoDOhdsxHP0Qfvgd8swjeIhoD93vh2aDoHEv8KpXozKUVpPUizMMJXYVEvyFTVW3Bu9KNcCKF4ApwitiGZ7BGzBmtaQg4Q500YW5dtqkHWPsjq9peu40Pv37Efnss3g2Ni84UlYQjGsSbJXfydnsAr7ZFs+KXYn8dfwMWkOTEF8SE04Skr+dkSdWM/LYGthxwvyGgGhocx3EDjTfTVsvvMrHtDVnGkrsCuQOX2FTFd1NW1u+kGWdo4dBgcrBvcEnuNc7SEFaX9zOXs+N3Rrx9ZZ4PLLOce/uZQw78SepPoHk/3s8Qx64xabTMhxKzsRo0rSODODkmRz6/WcVnSI8GBOdQH/DboKTNqCSdpk39q4PMf2KJ0UbCCHNLc7bO0ptuUO3puQO3zrIGZu8rlSDr66yznHMAF8+PvYaafmJ5CXcSLjqz6TRrbihUwOGHN5IvaVv4V2Qy/L2V9L08ce4oVfzSo5SdVprFv52iLfXHiEr3zx2vnPDeiwZ7Uejw6vY3+w3vE5vht2FYPA0z40z+FloOhgadgYL1wR2FuWl3+LTc4mdvKxWfvZqQmr+tURdqGG7ivXx65m0ZhIeBg/mDZxH14iuAOTt20fiC9PI3b4dn7huRE6dinfLllY9tta6pPUwbN7v7E86RzOVQB+3XfR120Uvtz34q+IgGdnRnLOPHWDO23ua58d3xkqEJSxZNKcufCek5l/HOGpcfG1hjYCntebjvR8zZ/Mcmgc257XBrxFVL4qirCxSFy7kzIcfYahfnwazZlF/5A1WS/HkG4vYcDiNFbsS+eNIGj8/0AqvE+t4OONjunv9TUN1BoATpjC+L7qCfT5dmT7+EfALuWxfrpw3t2TRHPlOXCDBv5aoa5NSWZM1Al5hUSEzN83k64NfM7jRYGb1m4WPuw/nfvqJpFmzMaakEHjrLYSPH48hMNAq5d6TcI7/rTnMxn3HaVe4i8EeuxnntRev+eabq/rremwwtWWhqQPrTe05oSMAUJkwvYzAD65dibg0/VZeTkO+E2YS/GsJVxoX72xqGvDO5J1hwqoJbE3eygMdHmBsl7EYj5/g5PQXyd6wAa+2bYhe+Bo+nTrVqJzpOQX8ujeZ1uHetNeHCNz2E/fsX848dRCDZxHa3RsV3QuamqdOuH5xKqcy8i/bT30f8+LtZbV2XL0SUXqEVHlpIHt+J5w5hSbBv5awdFy8M38YHaUmAW//mf08tuoxUnNTebnfywxvOJi0ha+T9s67KC8vIp59lqDbb0MZqtd5mnQuj593nWbnjk34x6+nl9pJC4/9YMqhAYoGDTqjmj0GTQegGvW8aBHyJ4bHM+nLHSV3GZ+XXWDk2SU7L7pn4HxrJ9DXg7M5hZeVwxUrEY6+V8TZU2gS/GsJS0bVOPuH0VGq22paeWIlk9dOpp5HPRYNX0TMnjMceWgEhSdPEnD99UQ8OQn3sLAql+dcXiEBBSmYDq9iy9LPGGb6m3+qdHCH/IBYPFvebg72Mf0qXIR8ZJcopn2/+7JgXlik+XTTSYouGeyRW1iEl7sbPh4Gl7i5rjKOHmnm7Ck0Cf61SGV3Rjr7h9FRqlpD1Frz3q73eG3ra7QLacfcts9QNO1tTv7yC55Nm9J40SL8ruhp8fG11hw4Hs/+TT+hD6+mU+E2AnQ8bsCVXsHkNxoAba+E2AF4BTWp0rmll1GLBy4L/Odl5BYy79bOtaZ16MiJ2pw9hSbBv5awJJ3j7B/GqrJWCqsqNcQ8Yx7Pb3ieH4/+yDWNhjHhUAvSp95iqijIAAAgAElEQVQNJhNhEycSMuZulKdn5Qc1FsCpPznwx/cUHlxF66IDtFKaHO3FVtWGjDa302nASDzD2+FZg8VaymvVGJQq8wLQMNCnTs1saUvO3g8nwb8WsDSd4+wfxqqwdgrLkoCXnJPM+FXj2Zm6k2e8R9Fj3nbOHlxGvSFDiJgyBc/oCt5vMlGYuItTm3+i8OBKmufuwM2YS3Pc+NvUlP+abmB9UQe26eYU4IHPLgOzWgYzMrJmq3SV16q5sVuUTeYJEhc4us+hMhL8awFL0zlV/TA6c+ewvVNYu1N3M27lOEjPYNHObvj+8iWmhg2JfuMN/AcPKvtN6ScwHlxF6t8r8EvYgH/RWWKBw7ohJ5qNIjG0F+P+8CO50Puyt1rrXCpq1VhrniBRNkf3OVRGgn8tYGk6pyofRmfvHLZnCmv50eU8t/YZrtvlxS2rgLythDz4IKH//hduPqVaTbln4ehaCg6uhCOr8cw4ijvgpgNZo9pzrkE/GnQZzhWdO+DjaeAfs1eSXFh+ea11LuW1aixp7ThzBcAVOHMKTYK/k6jJl6yidE5Z+7Vkkitn7xy2RwrLpE38d/t/+WXFW8xe6U3DE2fx7dmTyKnP4dWsGRTmFa9atQbjoVW4Je7ADRMF2pu9Xh3pPuxBaDaI5IIGDGsQgIfh4hROZcHd0ek4Z68AWEIuXuWT4O8EavolKy+dM6h1WLX36+ydw7bOp+YU5vDCz5OI/GQVs7ZqPIK9ifjPVAK6NUId/R7Wr4YTf4AxjyJlYHtRM9aZRnLAtxvR7fsztGM0xJiHYbYv5xgVTQXtDLlhZ68AVKY2XLxsySrBXyk1HFgAGIB3tdazL3l9DDAHiC9+6nWt9bvWOHZtUNMvWXnpnJrs19k7h22ZT43PjOe9V+9m1Pfx1M9RBF/ZlbA4MOx9DLaa58k5bmhCWMd/4tv6Sn7NbsauVBPD20fSJDGTV34+wDsbTlZapvLmogny9eD569tddo+GvWuwzl4BqIyrX7xsrcbBXyllAP4LXAWcAv5SSi3VWu+5ZNPPtdZja3q82sgaX7KycosTPt9e7f06+0gFsE0+dfvaTzk6cyY3HyvCGK6I7ZOMT/APZMaHs1534aeCVmwwtaNho1hmde9A24YBDAOGYQ7QT3+7y+KapqUXMEfVYJ29AlAZV7942Zo1av49gENa6yMASqnPgBuAS4O/KIetvmQ12a+jRyrYraZbkAMnN2La+ysbvlhG/c0FNPbQuPUqoFGfrvh0GE9SWC96v3uCnrEhDGsXyeR2ETSof/nvsDo1zdIXsPPnPOHz7Reds7VrsJb+bl2hAlARV7942Zo1gn8UcLLU41NAWbc33qiU6g8cACZorU+WsU2dZKsvmat+eW1a0zUVQcJ2OLLK3Fl7chMZJxSHtgcTkqnY0t6bta1v5FdjR67yiuJ/PbsRAWx9rnXJhGjlqUlNs6JztmYNtiq/W0dXAGrKVT//9mKvDt/vgU+11vlKqYeAxcBlQ06UUg8CDwI0Ll7HtC6w1ZesJvt1ZGeZVWu6WsOZI+Zgf3gVBYd/x7PwHAD7s5tg3NcKt4NppIbC/wZ25IDPXQxs3YD57SIZ1PrCOrWVBX6oWU2zonO2Zg22qr/bilJrzj6SxtUvXrZmjeAfDzQq9TiaCx27AGit00o9fBf4T1k70lq/DbwN5pW8rFA2l2Gr8cDV3a8jO8tqXNPNSoGja4pr92sgw9zIzPFpwLL8bqwvaEPggQxu2L8W1Bm+GmTA66ZxPNLgBvq2CMXbo3ozcA5qHcbHG09cNI+8pTXN8kb9xKfncucVja12N661WhGuMpLGmcfZO5o1gv9fQAulVCzmoH8bcEfpDZRSDbTWp4sfjgD2WuG4woYc2VlW5ZpuQTYc31A85n41lCxCHgix/aHvBPZ4d2HU54m0On2QR3Z8Q+OsZDa2dOf/BvmRrx9g05AHa1TmJdvi+XpL/EWBXwE3drMs+JQ31w7A11viubFbFKv2pdS4BmutVoSMpHF9NQ7+WmujUmossALzUM/3tda7lVLTgc1a66XAOKXUCMAInAHG1PS4wrYc2VlWaa62yAgJWy8E+5N/gqkQDF7Q+AoY8jyJIVfwQ0oY3WJD6dI4iMIdhxm38VUGn9pKYj1fXrrZnS3REeSeugsKy17VqirKCoYaWLUvxaL3lxf4wRxUV+1LsejmvMpYKw8uI2lcn1Vy/lrrH4EfL3luaqmfpwBTrHEsYR+O7Cy7LFdb35tpfTy5suAH+HQ1HFsH+ecABQ06Qq9H0E0HctCrHcv3Z7B8ayJ7Tp8FzjJhcBFNfv8R7/nz6ZeTy+fdGrF0UAK5eW3IPX4bmLyJssIFrabBMKqCG76qsp/KWCsPLiNpXJ/c4SvK5OjOspHN3Bh5dXLJ9AmsTDC/EBQL7UdD00GYmvQj0ehLw0AfTCbNHS/9Slp2AV0bB/H0Na25SqXBvGdI2rMHz57deblXDn/57Cc/dSAFKcMAN6td0GoaDCtbfNyaQdUaeXAZSeP6JPiLctm1syw/E46tvzAEM2Wf+XmfYGg6AJoOhNgBGOs34c9jZ1ixK5EV323HTcH6yYMxuClev6MrTUP9CNH5JM+bR/rnX+AeGorhxScZqz4nKTeZG6Oe4JfTjUjAuhe0mgbD82V4Yelu0nMvXoDFGYOqoysHouaUriDX6EhxcXF68+bNji6GsBVjAcRvvpC3P7UZdBG4e0OT3uZg33QgRHSA4sVMvth8klk/7uVsTiHeHm4MaBnGsHaRjOjUEHeDG1prMpZ8R/KcORSlpxP8zzs5MLork7Y8j7fBmwWDF9AprGaLqFfEWkMfnX0IpXBuSqktWuu4SreT4F/7OTKYXDh2Dn0DUniyZSId8reaa/mF2aDcoGEXiB0AzQZBdA/w8CYr38iqfcks353IuMEtaBXpz+8HUvhm6ymGt4+kf8swfD0vNFzzDx4kcdp0cjZvxqdzZyKmPsfnpj+Zu2UurYNb89rg14j0i7TLOQvhSJYGf0n71HKOHI+9fMMWNi3/kkl6B328dhNWkAG7INMvBv/Ot5sDfmw/8AkCILegiO93JLBiVyJrD6VSYDQRWs+Tk52jaBXpT/+WYfRvefGC6KbsbFLeeIMziz/A4OdH5IvT8R15PdM3vcjSw0sZ2mQoM/rOwMddOiKFKE2Cfy1n1/HYuenmkTjFqZzhaQcZ7gYpOoANpvasM7VnfVF7EvJCidrpw6SGrehZ4E1qWgYdoutj0prnluwitJ4Xd/ZswrB2EcTFBGNwU5cdSmtN5q+/kvTSLIynT1P/phsJf/xx0r2LuO+X+9mRsoOHOz3Mvzr9C6Uuf78QdZ0EfxtyhtytTcdjG/PNY+zP5+0TtoI2gYcfxPRhRmIP1po6sF83wnzL0wXx6blM+GI7WkPH6PosHdsXPy93fpkwgEbBPhUG7IKTJ0mcMYPsNb/j1bIlUa++gm/XruxN28u4ZeNIz0vn1QGvMjRmaM3PUYhaSoK/jTjL7e9WHY9tMpnvnj0f7I9vAGMuKANEdYN+T5jz9lFx4O7JT7NXVjh2XWsI8HZn3q2dS55rHOJb/uELCkh7913S3nobZTAQPvkpgu+8E+Xuzi/Hf+GZdc8Q4BnAB1d/QJuQNmXuwxkuyEI4Awn+NuIst7/XZL4ZAM4evxDsj66BnOJpmkJbQde7zCNyYvqCd8BFbzOZNDd2jeKN1YcxmsofVJCZZ6RZWL1Ki5G1fj1J01+k4Phx/IcPJ2LKZDwiItBa8+aON3lj+xt0DOvIgkELCPUJLXMfznJBFsIZSPC3EWe4/b1a883knIFja+Fw8Xj7s0fNz9eLhOZXQtNB5nH3AQ3LfPupszm8sfowP+9OIjUrH4ObwsvdjXyjqcztK2uBFCYlk/zybM79+BMeTRrT6J13qNevLwC5xlyeXfcsPx//mRHNRjC111S8DF7l7stZLshCOAMJ/jbiDLe/WzTfTGEenNx4Idif3mHeytPfXKPv+S9zKie0JZSRh88tKOL3gykE+njQs2kIbkrx3bZ4BrYKZ2i7CAa1DifA2+OyWjdU3ALRRiNnP/mElAWvoQsLCX10LCH334+blzm4J2YnMm7lOPad2cfEbhMZ025MpR27znBBFsJZSPC3EWe4/b2soOaGiZCM3bBuhznYn9gIxjxwczePsR84xZzKieoKhrLnsM/ILWTVvmRW7E5k9f4UcguLuK5jA3o2DaFhoA/bpg7F093tovdU5Y7QnG3bSJw2nfx9+/Dr14/I557Fs9T6DjtSdvDYysfIK8pj4eCFDGg0wKLfhzNckIVwFhL8i5XuCAz09UBrc5CrbqegM9z+bg52OTRWyfRz20kft130dttNoMqGX4HwdhB3nzmN06QPeJWfe8/ON+LnZf64/PO9Tfx9KoOIAC9u6hbN1e0j6REbXLLtpYH/vMqmizCePUvK3Lmkf/kV7pGRRC1YgP/Qqy6q0X9/+Hte2PAC4b7hvDv0XZoHNbf49+EMF2QhnIXc4cvlHYGX8vEwMGt0B9fJC2enwtE1HPvrR9yPrSFamdM8CTqYjboDjbtfS9zAkeAfUeFuTp7JYcXuRH7encSe0+fY/OyVeHsYWHcwFR9PA10aBeJWxhj8qtImExnffkvynFcoyswk+O67CXvkYdz8/EpdlLMJbvQbBfV+o3tkd+YOmEugd2CVjyWjfURtJ3f4VkFZufHSnL5TsCAHTpRazCTRPIIlxqs+CQ2782pKLMuyW5MfEMuk4a2Jq+Q8/jicxoxle9idYF7usHWkP/f1jSXfaMLbw0DfFmWPpqmOvP37SZw2ndytW/Hp1o3I56eyPNuPOQs3EZ+eiwK0Wx4+0Z9TUG8vpowruLbDs9UK/CArOwlxngR/LOvwc6pOwSIjnN5+YZnCk5ugqAAMntCoJwx+FpoOhgadaGhw53Hg8XJ2pbVmZ3wGy3clMrBVOD1ig/H3dsfT3Y0pV7dmWLtIYkL9rH8KWdmkvv46Zz78EENAAA1eeon6o0by3faEi1thHmn4Rn+Am1cKeYk3UHi2F3N/PsyNXZtYvUxC1CUS/Cm/I/DSbayhWmkHrSHtUKnx9mshP8P8WmQH6PmQeQhm417gWf5NUueZTJo/j51h+a5Eft6dSEJGHgY3RZCvJ91jgmjXMIC7e8UwZ8V+Zv+0z6rpEa01mcuXkzRrNsaUFAJvuYXwCeMxBJpr8qVbYQbfI3hHfYRSmtwT91KUY87vO9WFWAgXJcGfyhfSsFanYJVuMspMgqO/X5jf/ly8+fnAxtDuhpL57fGzLAWTbyziRFoOLSL8AXj0022cyy2kf8swJg5txZVtwgn09ax6Oasg/+hRkl6cQfaGDXi1bUP0wtfw6XTxFMvnA7tH4J94RS7BVBBCzsm70YUXzlNG5whRcxL8uXxkjjVG+5SlwpuM2taH4+sv1O6T95g38AkyL0Ie+7h5vH1QbJnj7cuSnW9k9f4Ulu9OZNW+ZOp5ubNh8mDc3BT/N6Y7saF+JSN4LC5nNX4Pprw80t5+m7R33kV5eRHxzDME3XE7ymC4bNsGgZ6keX2FZ/AGjFktyY2/HUwXgr2MzhHVIR39l5PgX8weHYGl0xXuGOmojtDXbRd9c3bBy4fAZLywCPmVL5hr9g06gdvlQbIyi9Yf5aWf9lFgNBHs58l1HRswrF1kyd2+7aPqW1ROS54v7dIv2QsRGcR+8iaFJ08ScN11RDz1JO5hYWW+NyM/g/AWH5GZtZWCtL7kJ18DuJk7fTGvcytfWlFVMq1H2ST424vW9A5IpUX2Zvq47eIKt734q1xMWrHfrRn0Gmuu2TfqCR5VS2skncvj592JLN+dyNPXtKFdw/q0jPTnHz0bM6xdJN3LmRa5PNW9Gar0lyw0J50xfy4mOmEnWQ0b0WLR/+F3xRXlvvdoxlEeXfko8TnxjIqeyG+nm1h9qUVRN8m0HmWT4G9L5xLMo3GKUzkfFySCBxwzRfBdUW/WmTqw3dCBySN70aaKH8LMvEI+2XSC5bsT2XYiHYCmYX6czTav/9q7WSi9m1VvSGZ1b4aas2I/BfkF3Hh4Lf/Y9zNu2sTiNsP5I244ayoI/BviN/DEmifwMHjw3tD36BrRlelDqlV0IS5ji2k9akMaSYK/NeVlFC9Cvtr8L3W/+XnfkJIO2p/z2jBtbVbJh2ayhR8arTX7EjM5l1tIz6YhuLu5Mf/XgzQPr8cTQ1syvH0kzcP9rXIapftA4tNzMShVUlMq/fqlAg/v4ZkdXxN7LpFNEW14s+NIkvxCUJnGcs/p470fM2fzHJoHNmfh4IU0rFf2hHFCVJe1p/WoLWkkCf41YSyAU39dCPbxW4oXIfcxL0Le5c7iRcjblyxCPhQY2sey3ZtMmm0n01mxO5HluxI5cSaHDlH1+f7Rvvh4GtgweTBBfp42ObXzH2JLPuTGM2dInvMKr6z9liSfQKb1HMPGBu1LXi/rS1ZQVMDMTTP55uA3DG40mFn9ZuHrUfkwVSGqytrTetSWNFKtDf42aZZpDUm7Sy1msh4Kc8yLkEd1g74Tihch7w7u5U8tXJEiky7Jz0/8YjtLtifgYVD0ahbKvwY046q2F6ZksFXgP6+yD7k2mUj/4kuS583DlJ1N2ohbGW/oQrq+8LEq60t2Ju8ME1ZNYGvyVh7o8ABju4zFTZU9H5AzqA1N/LrM2vNs1ZbZYWtl8Ldqsyz95IWFTI6sgexk8/MhLaDzP8yTosX0A5/qTTcAkFdYxNqDqSzflcjKfUn89Fh/Iut7c3NcIwa1Di+ZFtneKvqQ5+7eTeK06eT9/Te+PXoQOfU52jRvzguVBMoDZw8wbuU4UnNTubnxZD7/uQHzvvjJaYNqbWni13XWHM1XW2aHrZXB36rNst+mwc4vwS/cHOjPL2ZSP7rG5Tyams2cFftYvT+FnIIiArzdaRXpzw2vryM5M78kIDoi8EPZH3LfwlwePvIrx25eiyEoiIYvzyZgxIiSmTcr+pKtPLGSyWsn4+/hzz2xc3h9eT65heb9O2tQrS1NfGE9tWV22FoZ/K3aLOv3OPSdCOFtSm6uMqcBVla5CZmalc+ve5JoGOhD/5Zh+Hoa2HzsLKO6RDGsXSQpmfk8u2SX09QyL/qQa83AU9t4cNf3BBZkEXT77YSNfwxDQECl+9Fa896u91iwdQHtQ9qzYPACRr22yyWCam1p4tcV9kjROcN07dZQK4O/VZtl4RcvBF7VNEB8ei4rdpnH4G8+dgaThpu6RdO/ZRgRAd5senpISa25z+yVThUQzx/zw8/XcNP6T+mceojcpi2JfXkmPh3aV/JuszxjHs9veJ4fj/7INbHXMK33NLzdvV0mqNaWJn5dYM8UXW2YHdZ5e9lqYNKwVvh4XHxXrLWaZRWlAc5LPpdX8vNDH25m+g97yMgpZOyg5iwb15c5N3Useb30QiXOFhBNubn0Xv0lM5bNpmt+EhFTn6PL999YHPiTc5K5Z/k9/Hj0Rx7r+hiz+83G290bKD94OltQteVnSViXJd9NcUGtrPnboll2vjlZ3uyf8em5vLJiP8t3J3LqbA5bn7sKX093po1oT7CfJ7EWTIvsTLXMzFWrSJoxk8L4eOrfMILwSZNwD7X8prHdqbsZt3IcmYWZzB80nyGNL75ry1XyprWliV8XOFvlydlZJfgrpYYDCwAD8K7WevYlr3sBHwDdgDTgVq31MWscuzzWbJZVttLXeW+uOUzP2GDu6tWE8wukdWsSZPFxnCEgFsbHk/jSLLJ++w3P5s1o/MFi/Hr0qNI+fjr6E8+tf44Q7xA+vPpDWgVfXn5XCqq1oYlfFzhT5ckV1Dj4K6UMwH+Bq4BTwF9KqaVa6z2lNrsPOKu1bq6Uug14Gbi1pse2l8pW+nJTcGtcIyYNb01wDcbeOzIg6oIC0hYvJvWNNwEIe3wiIXffjfK0/HxM2sTr217nnZ3v0DW8K/MGzSPYO7jc7V0pqMpYf+fnDJUnV2KNmn8P4JDW+giAUuoz4AagdPC/AXih+OevgNeVUko76wLCpRSZdIULvVh7pklHBMTsTX+SOH06BYcP43/VlURMmYJHw6pNs5BTmMOUtVNYeXIlo1uM5tmez+JhcMwQVWuTsf6uwZVak87AGsE/CjhZ6vEpoGd522itjUqpDCAESLXC8W1KYa7Zm8q4TEUF+rB+8mC7l8lajKmpJP3nP5xb+j0eUVFEv/kG/oMGVWtfUzdMZfWp1TzZ/UnubHPnRR3Zrk7G+rsOV2pNOppTdfgqpR4EHgRo3Lixg0tj5uammHJ1G179eT95RlPJ867cnNRFRZz9/HNS5s3HlJdHyL8eIvShh3DzqX5u9NEujzKq+Sj6RFk4cZELkY5EURtZI/jHA41KPY4ufq6sbU4ppdyB+pg7fi+itX4beBsgLi7OaVJCD/RvSpi/V61oTubu3EXiCy+Qt3s3vr2uIPK5qXg1ja3xfpsENKFJQO1cVF06EkVtZI3g/xfQQikViznI3wbccck2S4G7gT+Am4CVrpDvL83Vm5NFGRmkLFjA2U8/wxAaQsNXXiHg2muqlZ6pa52f0pEoaqMaB//iHP5YYAXmoZ7va613K6WmA5u11kuB94APlVKHgDOYLxDCDrTWnFu6lKT/zKHo7FmC7ryTsHGPYvCv3tz/S7bFM+nLHRQWd4LEp+cy6csdQO3t/JSORFEbKWetgMfFxenNmzc7uhguLf/QIRKnTSfnr7/w7tSRBs8/j3fbtjXaZ+dpP5OeW3jZ84E+Hmx/fmiN9i2EqDml1BatdVxl2zlVh68oW1XTLKacHFLffJO0/1uEm58fkdOmEXjzTSi3ms/mUVbgr+h5IYRzkuDv5KoyxlxrTdZvv5H40ksYE05Tf/Rowp94HPfg8m+0Eq6hrvWzCNuT4F/MWb9c5Y0xf/yLHUz4fHtJWa8J0yTNmEnW6tV4tWxJ1Cev4Nu1q9XLE+Trwdmcy2v5Qb6144YuZyQ3mQlbqJWzelbV+S9XfHoumgtfriXbLh2xan/ljSUv0hoNJKdlsnnmXA5efS05f/5J+FNPEfvN1zYJ/ADPX98OD8PFI4Q8DIrnr29nk+MJma1S2IbU/HHuOzjLG2MO0Dn5AI/8/S3RWSlsbtKF2xbPxSMy0qblkZEv9ic3mQlbkOCPc3+5yhpjHpR3jgd3LmVg/HYS/EJ4ptcDbItoxT9tHPjPc/V7HlyN3GQmbEGCP8795Spd0048k8X1x/7gn3t+wt1UxEetruKLloMpNHgQ5QRlBeftO3FlcpOZsAUJ/jj/l2tklyiGqRROT3ud/L172RbRmoUdRnK6nnlxFWcpq3RM2oak2oQtSPDHub9cxrNnSZk7j/Qvv8Q9IoKo+fPZH94Ot58PoJysrM7cd+LqJNUmrE2CfzFn+3Jpk4mMb5eQ/MorFJ07R/CYMYSOHYuhnh8jgZFdox1dxMs4c9+JEOJiEvydUN7+/SROm07u1q34dO1K5PPP492qpaOLVSln7jsRQlxMxvk7kaKsbJJmv8zR0TdScPQoDWbOoMlHH7pE4Adz34mPh+Gi55ylP0IIcTGp+TsBrTWZK34madYsjElJBN58M2ETJ+AeZPni787AmftOhBAXk+BfBbYYxlhw/DiJL84ge906vFq3JnrBfHw6d7ZSie3P2fpOhBBlk+BvIWsPYzTl55P29jukvfMOysODiKefJuiO21Hu8icRQtie5PwtZM35VbLWruPI9SNI/e9/8b9yCE1//JHgu/4pgV8IYTcSbSxkjWGMhYmJJM1+mczly/GMiaHx++/h17u3tYpoM2Wlu0By+0K4Mgn+FqrJMEZdWMiZjz4mdeFCdFERYY+NI/i++3Dz9LRFUa2qrHTXpC93gILCogtLOcqdvEK4Fkn7WKi6wxhztm7l6I03kfzyy/h0j6PpD98T+u9/u0Tgh7LTXYUmXRL4z7PnFMNLtsXTZ/ZKYicvo8/slU4x9bYQrkZq/haq6jBG45kzJL/yKhnffIN7gwZEv76QekOGoJQqc/uqsOfkaVVJa9njTl6ZP0gI65DgXwWWDGPUJhPpX31F8qtzMWVnE3L/fYQ+/DBuvr5WKYO9g19F6wmUta2tyfxBQliHpH2sKG/PHo7dfjuJU5/Hu0ULmi75lvAnnrBa4Af7r+pUVrrLw01dtpqXve7klfmDhLAOqflbQVFWFikLXuPsxx9jCAyk4cuzCRgxwiopnkvZO/iVl+4q6zl71Lxl/iAhrEOCfw1orTn3448kz34ZY2oqgbfeQviECRjq17fZMR0R/MpLdzkizeLsay8I4Sok+FdT/tGjJL34Itkb/sC7XTui3/gvPh062Py4dT34yfxBQliHBP8qMuXlkfrWW5x59z2UtzcRzz1L0G23oQyGyt9sBRL8ZP4gIaxBgn8VZK5eTdKMmRSeOkXA9dcT8eQk3MPC7F4OCX5CiJqS4G+BwoQEkmbNIvOXX/Fs1ozGixfj17OHo4slhBDVJsG/Emc++JDkefNAa8ImTiRkzN0oF7k7VwghyiPBvxKm/Dz8evcm8ukpeERJqkUIUTtI8K9EyH33odzkXjghRO1So6imlApWSv2ilDpY/H+Z6w4qpYqUUtuL/y2tyTHtTQK/EKI2qmlkmwz8prVuAfxW/LgsuVrrzsX/RtTwmEIIIWqopsH/BmBx8c+LgZE13J8QQgg7qGnwj9Bany7+ORGIKGc7b6XUZqXURqWUXCCEEMLBKu3wVUr9CkSW8dIzpR9orbVSSpexHUATrXW8UqopsFIptVNrfbiMYz0IPAjQuHHjSgsvhBCieioN/lrrK8t7TSmVpJRqoLU+rZRqACSXs4/44v+PKKVWA12Ay4K/1vpt4G2AuLi48i4kQgghaqimaZ+lwN3FP98NfHfpBkqpIKWUV/HPoUAfYE8NjyuEEKIGahr8ZwNXKaUOAlcWP0YpFaeUeoRV8/oAAAUJSURBVLd4mzbAZqXUDmAVMFtrLcFfCCEcqEY3eWmt04AhZTy/Gbi/+OcNgO3nOhZCCGExuYNJCCHqIAn+QghRB8ncPqUs2RZfpxdJEULUHRL8iy3ZFn/R8ojx6blM+WYn4Ji1aoUQwpYk7VNszor9F62LC5BbWMScFfsdVCIhhLAdCf7FEtJzq/S8EEK4Mgn+xRoG+lTpeSGEcGUS/ItNGtYKHw/DRc/5eBiYNKyVg0okhBC2Ix2+xc536spoHyFEXSDBv5SRXaIk2Ash6gRJ+wghRB0kwV8IIeogCf5CCFEHSfAXQog6SIK/EELUQRL8hRCiDlJaO+dSuUqpFOC4o8tRLBRIdXQhrEDOw/nUlnOR83AeTbTWYZVt5LTB35kopTZrreMcXY6akvNwPrXlXOQ8XI+kfYQQog6S4C+EEHWQBH/LvO3oAliJnIfzqS3nIufhYiTnL4QQdZDU/IUQog6S4F+KUmq4Umq/UuqQUmpyGa97KaU+L359k1Iqxv6lrJwF5zFRKbVHKfW3Uuo3pVQTR5SzMpWdR6ntblRKaaWUU47SsOQ8lFK3FP9NdiulPrF3GS1lwWersVJqlVJqW/Hn6xpHlLMiSqn3lVLJSqld5byulFKvFZ/j30qprvYuo11oreWfOfVlAA4DTQFPYAfQ9pJtHgb+V/zzbcDnji53Nc9jEOBb/PO/XfU8irfzB34HNgJxji53Nf8eLYBtQFDx43BHl7sG5/I28O/in9sCxxxd7jLOoz/QFdhVzuvXAD8BCrgC2OToMtvin9T8L+gBHNJaH9FaFwCfATdcss0NwOLin78ChiillB3LaIlKz0NrvUprnVP8cCMQbecyWsKSvwfAi8DLQJ49C1cFlpzHA8B/tdZnAbTWyXYuo6UsORcNBBT/XB9IsGP5LKK1/h04U8EmNwAfaLONQKBSqoF9Smc/EvwviAJOlnp8qvi5MrfRWhuBDCDELqWznCXnUdp9mGs5zqbS8yhujjfSWi+zZ8GqyJK/R0ugpVJqvVJqo1JquN1KVzWWnMsLwJ1KqVPAj8Cj9imaVVX1O+SSZCWvOkwpdScQBwxwdFmqSinlBswFxji4KNbgjjn1MxBzK+x3pVQHrXW6Q0tVPbcDi7TWryqlegEfKqXaa61Nji6YuJjU/C+IBxqVehxd/FyZ2yil3DE3a9PsUjrLWXIeKKWuBJ4BRmit8+1Utqqo7Dz8gfbAaqXUMcy52aVO2Olryd/jFLBUa12otT4KHMB8MXA2lpzLfcAXAFrrPwBvzPPluBKLvkOuToL/BX8BLZRSsUopT8wduksv2WYpcHfxzzcBK3VxD5ETqfQ8lFJdgLcwB35nzS9XeB5a6wytdajWOkZrHYO572KE1nqzY4pbLks+V0sw1/pRSoViTgMdsWchLWTJuZwAhgAopdpgDv4pdi1lzS0F7ioe9XMFkKG1Pu3oQlmbpH2Kaa2NSqmxwArMoxre11rvVkpNBzZrrZcC72Fuxh7C3GF0m+NKXDYLz2MOUA/4sri/+oTWeoTDCl0GC8/D6Vl4HiuAoUqpPUARMElr7WwtSkvP5XHgHaXUBMydv2OcrYKklPoU88U2tLhv4nnAA0Br/T/MfRXXAIeAHOAex5TUtuQOXyGEqIMk7SOEEHWQBH8hhKiDJPgLIUQdJMFfCCHqIAn+QghRB0nwF+L/26kDAQAAAABB/taDXBDBkPwBhuQPMBRDqbkzKGFPBgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x107326320>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# However, as soon as more noise is on one variable or the other, the outcome becomes \n",
"# skewed in favor of the linear regression predicting the noisiest variable\n",
"n = 100\n",
"x = np.linspace(0, 1, n)\n",
"y = np.linspace(0, 1, n)\n",
"\n",
"x_o = x + np.random.normal(0, 0.1, n)\n",
"y_o = y + np.random.normal(0, 0.4, n)\n",
"\n",
"plt.plot(x, y, '--')\n",
"plt.scatter(x_o, y_o)\n",
"\n",
"lin_x = LinearRegression(fit_intercept=True)\n",
"lin_x.fit(x_o[:, None], y_o)\n",
"plt.plot(x, lin_x.predict(x[:, None]), label='y=f(x)')\n",
"\n",
"lin_y = LinearRegression(fit_intercept=True)\n",
"lin_y.fit(y_o[:, None], x_o)\n",
"plt.plot(lin_y.predict(y[:, None]), y, label='x=f(y)')\n",
"\n",
"w, b = odr_fit(x_o, y_o)\n",
"plt.plot(x, x*w + b, label='odr')\n",
"\n",
"plt.legend()"
]
},
{
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment