{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# REGRESIÓN LINEAL\n", "### Introducción\n", "Importamos los paquetes de Python que vamos a necesitar:" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [], "source": [ "import pandas as pd, numpy as np, matplotlib.pyplot as plt\n", "np.set_printoptions(suppress=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "y creamos el dataset de alturas y pesos:" ] }, { "cell_type": "code", "execution_count": 139, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
alturapeso
018087
116448
217280
316760
417685
518498
616155
719186
817866
917575
1016060
1117776
\n", "
" ], "text/plain": [ " altura peso\n", "0 180 87\n", "1 164 48\n", "2 172 80\n", "3 167 60\n", "4 176 85\n", "5 184 98\n", "6 161 55\n", "7 191 86\n", "8 178 66\n", "9 175 75\n", "10 160 60\n", "11 177 76" ] }, "execution_count": 139, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.DataFrame({'altura':[180, 164, 172, 167, 176, 184, 161, 191, 178, 175, 160, 177], \n", " 'peso':[87, 48, 80, 60, 85, 98, 55, 86, 66, 75, 60, 76]})\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Justificación del nombre \"regresión **lineal**\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vemos que los puntos (altura,2*altura+5) están alineados:" ] }, { "cell_type": "code", "execution_count": 140, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAUEUlEQVR4nO3df6jd933f8edLV5JtylIL+RpsS5GlxQrYJlGrU3EZ87IKr/YCi1lDOxUPApnRGtSSZrQDky7MK4EsTVa24a14aTazKXFU5rWuwaw2RAmGXGn3plJqGXm5kaNZdVffaPKGYJN8dd/743xVn0jn6p5r3XvPPV89H/BF3/P5/uDz4QOv+9HnfM73m6pCktQu64ZdAUnS8jPcJamFDHdJaiHDXZJayHCXpBYy3CWphRYN9yQ3Jzma5HiSE0meaMo/nOQ7Sf40yR8leV/PNY8nmUnyWpKHVrIBkqSrZbF17kkC/ERVnU+yAXgZ+DTwr4Ffr6pvJfkksL2q/kmSe4GvA3uAO4GXgJ1VdWklGyJJeteiI/fqOt983NBsBXwQ+HZT/iLw8Wb/EeCZqrpQVa8DM3SDXpK0StYPclKSMWAa+ADwZFUdSfIK8DHgD4FfALY2p98FTPZcfqYpW9Btt91Wd99999JqLkk3uOnp6R9V1Xi/YwOFezOlsivJrcB/SXI/8EngXyX5HPAccLE5Pf1ucWVBkv3AfoD3v//9TE1NDVIVSVIjyemFji1ptUxVvQ0cBh6uqpNV9XNVtZvuHPsPmtPO8O4oHmAL8Gafez1VVZ2q6oyP9/3DI0l6jwZZLTPejNhJcgvwIHAyye1N2TrgN4HfbS55DtiX5KYk24F7gKMrUXlJUn+DTMvcATzdzLuvAw5V1fNJPp3kQHPOs8C/B6iqE0kOAa8Cc8ABV8pI0upadCnkauh0OuWcuyQtTZLpqur0O+YvVCWphQx3SWohw12ShmT69Dme/OYM06fPLfu9B1rnLklaXtOnz/HoVya5ODfPxvXrOPjYBLu3bVq2+ztyl6QhmDx1lotz88wXvDM3z+Sps8t6f8NdkoZgYsdmNq5fx1hgw/p1TOzYvKz3d1pGkoZg97ZNHHxsgslTZ5nYsXlZp2TAcJekodm9bdOyh/plTstIUgsZ7pLUQoa7JLWQ4S5JLWS4S1ILGe6S1EKGuyS1kOEuSS1kuEtSCxnuktRChrsktZDhLkktZLhLUgsZ7pLUQoa7JLXQouGe5OYkR5McT3IiyRNN+a4kk0mOJZlKsqfnmseTzCR5LclDK9kASdLVBnlZxwVgb1WdT7IBeDnJC8A/A56oqheSfBT4IvA3k9wL7APuA+4EXkqys6ourVAbJElXWHTkXl3nm48bmq2a7X1N+U8Cbzb7jwDPVNWFqnodmAH2IElaNQO9Zi/JGDANfAB4sqqOJPk14L8m+RLdPxJ/rTn9LmCy5/IzTZkkaZUM9IVqVV2qql3AFmBPkvuBTwGfqaqtwGeA32tOT79bXFmQZH8zVz81Ozv73movSeprSatlqupt4DDwMPAJ4Nnm0O/z7tTLGWBrz2VbeHfKpvdeT1VVp6o64+PjS6y2JOlaBlktM57k1mb/FuBB4CTdwP5Ic9pe4PvN/nPAviQ3JdkO3AMcXe6KS5IWNsic+x3A0828+zrgUFU9n+Rt4F8mWQ/8P2A/QFWdSHIIeBWYAw64UkaSVleqrpoOX3WdTqempqaGXQ1JGilJpquq0++Yv1CVpBYy3CWphQx3SWohw12SWshwl6QWMtwlqYUMd0lqIcNdklrIcJekFjLcJamFDHdJaiHDXZJayHCXpBYy3CWphQx3SWohw13Sqpg+fY4nvznD9Olzw67KDWGQNzFJ0nWZPn2OR78yycW5eTauX8fBxybYvW3TsKvVao7cJa24yVNnuTg3z3zBO3PzTJ46O+wqtZ7hLmnFTezYzMb16xgLbFi/jokdm4ddpdZzWkbSitu9bRMHH5tg8tRZJnZsdkpmFRjuklbF7m2bDPVV5LSMJLWQ4S5JLbRouCe5OcnRJMeTnEjyRFP+jSTHmu2HSY71XPN4kpkkryV5aCUbIEm62iBz7heAvVV1PskG4OUkL1TV37t8QpIvA/+72b8X2AfcB9wJvJRkZ1VdWv7qS5L6WXTkXl3nm48bmq0uH08S4BeBrzdFjwDPVNWFqnodmAH2LGutJUnXNNCce5KxZtrlLeDFqjrSc/gB4C+q6vvN57uAN3qOn2nKrrzn/iRTSaZmZ2ffW+0lSX0NFO5VdamqdgFbgD1J7u85/Eu8O2oHSL9b9LnnU1XVqarO+Pj4UuosSVrEklbLVNXbwGHgYYAk64GfB77Rc9oZYGvP5y3Am9dVS0nSkgyyWmY8ya3N/i3Ag8DJ5vCDwMmqOtNzyXPAviQ3JdkO3AMcXd5qS5KuZZDVMncATycZo/vH4FBVPd8c28ePT8lQVSeSHAJeBeaAA66UkaTVlaqrpsNXXafTqampqWFXQ5JGSpLpqur0O+YvVCWphQx3SWohw12SWshwl6QWMtwlqYUMd0lqIcNdklrIcJekFjLcpRaZPn2OJ785w/Tpc8OuiobMF2RLLTF9+hyPfmWSi3PzbFy/joOPTfhC6huYI3epJSZPneXi3DzzBe/MzTN56uywq6QhMtyllpjYsZmN69cxFtiwfh0TOzYPu0oaIqdlpJbYvW0TBx+bYPLUWSZ2bHZK5gZnuEstsnvbJkNdgNMyktRKhrsktZDhLkktZLhLUgsZ7pLUQoa7JLWQ4S5JLWS4S0Pkg760UvwRkzQkPuhLK2nRkXuSm5McTXI8yYkkT/Qc+9UkrzXlX+wpfzzJTHPsoZWqvDTKfNCXVtIgI/cLwN6qOp9kA/BykheAW4BHgA9V1YUktwMkuRfYB9wH3Am8lGRnVV1amSZIo+nyg77emZv3QV9adouGe1UVcL75uKHZCvgU8IWqutCc91ZzziPAM03560lmgD3Ad5a57tJI80FfWkkDfaGaZCzJMeAt4MWqOgLsBB5IciTJt5L8THP6XcAbPZefacquvOf+JFNJpmZnZ6+vFdKI2r1tEwd+9gMGu5bdQOFeVZeqahewBdiT5H66o/5NwATwG8ChJAHS7xZ97vlUVXWqqjM+Pv6eGyBJutqSlkJW1dvAYeBhuiPyZ6vrKDAP3NaUb+25bAvw5rLUVhoRLnHUsC06555kHHinqt5OcgvwIPDP6c7D7wUOJ9kJbAR+BDwHfC3Jv6D7heo9wNEVqr+05rjEUWvBIKtl7gCeTjJGd6R/qKqeT7IR+GqSV4CLwCeaL19PJDkEvArMAQdcKaMbSb8ljoa7Vtsgq2W+B/xUn/KLwN9f4JrPA5+/7tpJI8gljloL/IWqtMxc4qi1wHCXVoDvMtWw+eAwaYlcCaNR4MhdWgJXwmhUOHKXlsCHfWlUGO7SElxeCTMWXAmjNc1pGWkJXAmjUWG4S0vkShiNAqdlpIarYNQmjtwlXAWj9nHkLuEqGLWP4S7hKhi1j9MyEq6CUfsY7lLDVTBqE6dlJKmFDHdJaiHDXZJayHCXpBYy3CWphQx3SWohw12SWshwl6QWMtwlqYUWDfckNyc5muR4khNJnmjK/2mSP0tyrNk+2nPN40lmkryW5KGVbIAk6WqDPH7gArC3qs4n2QC8nOSF5tjvVNWXek9Oci+wD7gPuBN4KcnOqrq0nBWXJC1s0ZF7dZ1vPm5otrrGJY8Az1TVhap6HZgB9lx3TSVJAxtozj3JWJJjwFvAi1V1pDn0K0m+l+SrSS4/ceku4I2ey880ZVfec3+SqSRTs7Oz19EESdKVBgr3qrpUVbuALcCeJPcD/xb4q8Au4M+BLzenp98t+tzzqarqVFVnfHz8PVVe7eTr7qTrt6RH/lbV20kOAw/3zrUn+XfA883HM8DWnsu2AG9eZz11g/B1d9LyGGS1zHiSW5v9W4AHgZNJ7ug57e8CrzT7zwH7ktyUZDtwD3B0eauttvJ1d9LyGGTkfgfwdJIxun8MDlXV80n+Y5JddKdcfgj8Q4CqOpHkEPAqMAcccKWMBnX5dXfvzM37ujvpOqTqWgtfVken06mpqalhV0NrxPTpc77uThpAkumq6vQ75mv2tOb4ujvp+vn4AUlqIcNdklrIcJekFjLcJamFDHdJaiHDXZJayHCXpBYy3CWphQx3SWohw12SWshwl6QWMtwlqYUMd0lqIcNdklrIcJekFjLcJamFDHdJaiHDXZJayHCXpBYy3CWphQx3SWohw12SWmjRcE9yc5KjSY4nOZHkiSuO/3qSSnJbT9njSWaSvJbkoZWouCRpYesHOOcCsLeqzifZALyc5IWqmkyyFfhbwP+4fHKSe4F9wH3AncBLSXZW1aUVqL8kqY9FR+7Vdb75uKHZqvn8O8A/7vkM8AjwTFVdqKrXgRlgz/JVWZK0mIHm3JOMJTkGvAW8WFVHknwM+LOqOn7F6XcBb/R8PtOUXXnP/UmmkkzNzs6+x+pLkvoZKNyr6lJV7QK2AHuSfAj4LPC5Pqen3y363POpqupUVWd8fHwpddY1TJ8+x5PfnGH69LlhV0XSEA0y5/6XqurtJIfpTr1sB44ngW7ofzfJHroj9a09l20B3lyW2uqapk+f49GvTHJxbp6N69dx8LEJdm/bNOxqSRqCQVbLjCe5tdm/BXgQ+JOqur2q7q6qu+kG+k9X1f8EngP2JbkpyXbgHuDoirVAf2ny1Fkuzs0zX/DO3DyTp84Ou0qShmSQkfsdwNNJxuj+MThUVc8vdHJVnUhyCHgVmAMOuFJmdUzs2MzG9et4Z26eDevXMbFj87CrJGlIUnXVdPiq63Q6NTU1NexqtML06XNMnjrLxI7NTslILZdkuqo6/Y4tac5da9/ubZsMdUk+fkCS2shwl6QWMtwlqYUMd0lqIcNdklrIcJekFjLcJamFDHdJaiHDXZJayHCXpBYy3CWphQx3SWohw12SWshwX0G+8k7SsPjI3xXiK+8kDZMj9xXiK+8kDZPhvkIuv/JuLPjKO0mrzmmZFbJ72yYOPjbhK+8kDYXhvoJ85Z2kYXFaRpJayHCXpBZaNNyT3JzkaJLjSU4keaIp/60k30tyLMkfJ7mz55rHk8wkeS3JQyvZAEnS1QYZuV8A9lbVh4FdwMNJJoDfrqoPVdUu4HngcwBJ7gX2AfcBDwP/JsnYitRektTXouFeXeebjxuararq//Sc9hNANfuPAM9U1YWqeh2YAfYsY50lSYsYaM49yViSY8BbwItVdaQp/3ySN4BHaUbuwF3AGz2Xn2nKJEmrZKBwr6pLzfTLFmBPkvub8s9W1VbgIPArzenpd4srC5LsTzKVZGp2dva91V6S1NeSVstU1dvAYbpz6b2+Bny82T8DbO05tgV4s8+9nqqqTlV1xsfHl1INSdIiBlktM57k1mb/FuBB4GSSe3pO+xhwstl/DtiX5KYk24F7gKPLW21J0rUM8gvVO4CnmxUv64BDVfV8kv+c5IPAPHAa+GWAqjqR5BDwKjAHHKiqSytTfUlSP6m6ajp81XU6nZqamhp2NSRppCSZrqpOv2P+QlWSWshwl6QWMtwlqYVuiHD3XaaSbjStf5677zKVdCNq/cjdd5lKuhG1Ptx9l6mkG1Hrp2V8l6mkG1Hrwx18l6mkG8/IT8u4EkaSrjbSI3dXwkhSfyM9cncljCT1N9Lh7koYSepvpKdlXAkjSf2NdLiDK2EkqZ+RnpaRJPVnuEtSCxnuktRChrsktZDhLkktZLhLUgulqoZdB5LMAqev4xa3AT9apuoMi21YG2zD2mAbBrOtqsb7HVgT4X69kkxVVWfY9bgetmFtsA1rg224fk7LSFILGe6S1EJtCfenhl2BZWAb1gbbsDbYhuvUijl3SdKPa8vIXZLUY82He5KvJnkryStXlP9qkteSnEjyxZ7yx5PMNMceWv0aX20pbUhyd5L/m+RYs/3ucGr94/q1Ick3eur5wyTHeo6NRD8s1IYR64ddSSabek4l2dNzbM31AyytHSPWFx9O8p0kf5rkj5K8r+fY6vZFVa3pDfgbwE8Dr/SU/SzwEnBT8/n25t97gePATcB24AfA2Ii14e7e89bK1q8NVxz/MvC5UeuHa7RhZPoB+GPgbzf7HwUOr+V+eA/tGKW++G/AR5r9TwK/Nay+WPMj96r6NvC/rij+FPCFqrrQnPNWU/4I8ExVXaiq14EZYA9DtsQ2rEkLtAGAJAF+Efh6UzRK/QD0bcOatEAbCrg8QvxJ4M1mf032Ayy5HWvSAm34IPDtZv9F4OPN/qr3xZoP9wXsBB5IciTJt5L8TFN+F/BGz3lnmrK1aKE2AGxP8idN+QPDquASPAD8RVV9v/k8Sv1w2ZVtgNHph18DfjvJG8CXgMeb8lHrh4XaAaPTF68AH2v2fwHY2uyvel+MarivBzYBE8BvAIeakVf6nLtWlwMt1IY/B95fVT8F/CPga73zdmvUL/HjI95R6ofLrmzDKPXDp4DPVNVW4DPA7zXlo9YPC7VjlPrik8CBJNPAXwEuNuWr3hejGu5ngGer6ygwT/c5Dmd49y8lwBbW7n/t+rah+W/bWYCqmqY7N7dziPW8piTrgZ8HvtFTPEr90LcNI9YPnwCebfZ/n3f/uz9S/cAC7Rilvqiqk1X1c1W1m+5g4QfNoVXvi1EN9z8A9gIk2QlspPuAnueAfUluSrIduAc4OrRaXlvfNiQZTzLWlO+g24ZTQ6vl4h4ETlbVmZ6yUeoH6NOGEeuHN4GPNPt7gctTS6PWD33bMUp9keT25t91wG8Cl1f2rH5fDPsb5wG+kf463f+WvUP3r98/oBuE/4nu/NZ3gb0953+W7l/L12i+eR/2tpQ20P0C5gTdb9a/C/ydYdd/oTY05f8B+OU+549EPyzUhlHqB+CvA9NNXY8Au9dyPyy1HSPWF58G/nuzfYHmh6LD6At/oSpJLTSq0zKSpGsw3CWphQx3SWohw12SWshwl6QWMtwlqYUMd0lqIcNdklro/wO9/ip1AOet1wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "df['2*altura + 5'] = 2*df.altura.values+5\n", "plt.plot(df.altura.values, df['2*altura + 5'].values, '.'); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cálculo del error RMSE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "El RMSE entre el peso y 2*altura+5 es:" ] }, { "cell_type": "code", "execution_count": 143, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
alturapeso2*altura + 5xi-yi(xi-yi)^2
018087365-27877284.0
116448333-28581225.0
217280349-26972361.0
316760339-27977841.0
417685357-27273984.0
518498373-27575625.0
616155327-27273984.0
719186387-30190601.0
817866361-29587025.0
917575355-28078400.0
1016060325-26570225.0
1117776359-28380089.0
media78220.3
RMSE279.7
\n", "
" ], "text/plain": [ " altura peso 2*altura + 5 xi-yi (xi-yi)^2\n", "0 180 87 365 -278 77284.0\n", "1 164 48 333 -285 81225.0\n", "2 172 80 349 -269 72361.0\n", "3 167 60 339 -279 77841.0\n", "4 176 85 357 -272 73984.0\n", "5 184 98 373 -275 75625.0\n", "6 161 55 327 -272 73984.0\n", "7 191 86 387 -301 90601.0\n", "8 178 66 361 -295 87025.0\n", "9 175 75 355 -280 78400.0\n", "10 160 60 325 -265 70225.0\n", "11 177 76 359 -283 80089.0\n", "media 78220.3\n", "RMSE 279.7" ] }, "execution_count": 143, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_rmse = df.copy()\n", "df_rmse['xi-yi'] = df_rmse.peso - (2*df_rmse.altura + 5)\n", "df_rmse['(xi-yi)^2'] = (df_rmse.peso - (5+2*df_rmse.altura))**2\n", "df_rmse.loc['media', ['altura', 'peso', '2*altura + 5', 'xi-yi', '(xi-yi)^2']] = ['', '', '', '', np.mean(df_rmse['(xi-yi)^2'])]\n", "df_rmse.loc['RMSE', ['altura', 'peso', '2*altura + 5', 'xi-yi', '(xi-yi)^2']] = ['', '', '', '', np.sqrt(df_rmse.loc['media', '(xi-yi)^2'])]\n", "df_rmse['(xi-yi)^2'] = np.round(df_rmse['(xi-yi)^2'], 1)\n", "df_rmse" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternativamente, se descarga el error MSE de Scikit-Learn y se hace su raíz cuadrada:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "279.6789826449841" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.metrics import mean_squared_error as mse\n", "rmse = lambda x, y: np.sqrt(mse(x, y))\n", "rmse(df['peso'], df['2*altura + 5'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Y vemos que coincide el número (en la tabla se ha redondeado a un decimal). \n", "### Fijar un modelo lineal\n", "Calculamos ahora la óptima recta de regresión para altura vs. peso:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import LinearRegression as lm\n", "X, y = df.altura.values.reshape(-1, 1), df.peso.values\n", "fitted_mod = lm().fit(X, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para calcular los coeficientes que se multiplican con las variables (en este caso solamente A) usamos el atributo *coef_*, y para el coeficiente libre, *intercept_*:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([1.34224736]), -160.21547955772695)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fitted_mod.coef_, fitted_mod.intercept_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Así pues, la recta óptima es $peso = 1.34224736 \\cdot altura - 160.21547955772695$.\n", "### Un dataset para el cual un modelo lineal es inapropiado:\n", "Tomemos ahora el dato data1.csv, y representemos el dato contenido en él:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZzElEQVR4nO3dfYxd9X3n8fdnxrgpadk4NiAHP8XLgwSoTfEs8hZtCwlNaGvFq01TGZEV20KtVt40BFUBbyVQIyGhbMsSackfXsKSqo6padwFsWoDoSAUKYMz45AWkxosw8AEb2zMZHezaW2P57t/3Dtwfbl37sN5PvfzkqyZOfdhvnc893N+8z2/8zuKCMzMrF7Gii7AzMzS53A3M6shh7uZWQ053M3MasjhbmZWQ8uKLgBg1apVsWHDhqLLMDOrlOnp6bci4vxOt5Ui3Dds2MDU1FTRZZiZVYqkmW63uS1jZlZDDnczsxpyuJuZ1ZDD3cyshhzuZmY15HA3M6shh7uZWY6mZ+Z44JnDTM/MZfp9SjHP3cxsFEzPzHHTg5Ocml9g+bIxdt+6mU3rV2TyvTxyNzPLyeSRE5yaX2Ah4PT8ApNHTmT2vRzuZmY52bxxJcuXjTEuOGfZGJs3rszse7ktY2aWk03rV7D71s1MHjnB5o0rM2vJgMPdzCxXm9avyDTUF7ktY2ZWQz3DXdJDko5JerFl20ckTUp6QdKUpKtbbtsp6bCkQ5I+kVXhZmbWXT8j94eBG9q2fQn4k4j4CHBX82skXQ5sA65oPuYrksZTq9bMzPrSM9wj4jng7fbNwHnNz/8F8Gbz863AIxFxMiJeBQ4DV2NmZrka9oDqbcA3Jf0pjR3ELze3XwRMttxvtrntPSRtB7YDrFu3bsgyzMysk2EPqP4B8PmIWAt8Hvhqc7s63Dc6PUFE7IqIiYiYOP/8jleJMjOzIQ0b7jcD+5qfP8q7rZdZYG3L/dbwbsvGzMxyMmy4vwn8avPzjwKvND9/HNgm6WckfRi4BNifrEQzMxtUz567pD3AtcAqSbPA3cDvAV+WtAz4Z5q984g4KGkv8BIwD+yIiDMZ1W5mZl30DPeIuLHLTZu63P8e4J4kRZmZWTI+Q9XMrIYc7mZmNeRwNzOrIYe7mVkNOdzNzGrI4W5mVkMOdzOzGnK4m5nVkMPdzKyGHO5mZjXkcDczqyGHu5lZTqZn5njgmcNMz8xl/r2GvRKTmZkNYHpmjpsenOTU/ALLl42x+9bNbFq/IrPv53A3M8vB5JETnJpfYCHg9PwCk0dOvLN988aVqQe9w93MLAebN65k+bIxTs8vcM6yMVacuzzTkbzD3cwsB5vWr2D3rZvfGal3Gsk73M3MKmjT+hVnBXjrSH7zxpWpfq+es2UkPSTpmKQX27Z/VtIhSQclfall+05Jh5u3fSLVas3MSmiYWTCLI/nbP35ZJgdX+xm5Pwz8V+DPFzdIug7YCvxCRJyUdEFz++XANuAK4EPAtyRd6uuomlldJZkF0z6ST1PPkXtEPAe83bb5D4B7I+Jk8z7Hmtu3Ao9ExMmIeBU4DFydYr1mZqXSbRZMJ1WY534p8G8k3QP8M/BHEfFd4CJgsuV+s81t7yFpO7AdYN26dUOWYWaWn+mZufdMXWyfBdOtd16Vee7LgBXAZuBfAXslbQTU4b7R6QkiYhewC2BiYqLjfczMyqJbOLfPgukW2FnPjmk3bLjPAvsiIoD9khaAVc3ta1vutwZ4M1mJZmbFWyqc++md9zvCT8uw4f4/gI8Cz0q6FFgOvAU8Dnxd0n00DqheAuxPo1AzsyIlDed+R/hp6RnukvYA1wKrJM0CdwMPAQ81p0eeAm5ujuIPStoLvATMAzs8U8bM6iCNcM5ydkw7NTK5WBMTEzE1NVV0GWZmlSJpOiImOt3mJX/NzGrI4W5mNqA856sPy2vLmJkNIO/56sPyyN3MbACDnJFaJIe7mdkAFqdEjoueUyKLbN+4LWNmNoB+p0QW3b5xuJuZDaif+ep5LzfQzm0ZM7MMDNK+yYJH7mZmCXVaLTLv5QbaOdzNzBJYqree53ID7dyWMTNLoKxTIx3uZmYJFN1b78ZtGTOzBIrurXfjcDczS6jI3no3bsuYmdWQw93MrIYc7mZmNdQz3CU9JOlY85J67bf9kaSQtKpl205JhyUdkvSJtAs2M7Pe+hm5Pwzc0L5R0lrg14DXW7ZdDmwDrmg+5iuSxlOp1MysTRUumlGUnrNlIuI5SRs63PRfgC8Aj7Vs2wo8EhEngVclHQauBr6TvFQzs3cVvepi2Q3Vc5f0SeCHEfH9tpsuAt5o+Xq2ua3Tc2yXNCVp6vjx48OUYWYjrKxnhpbFwOEu6Vzgj4G7Ot3cYVt0ep6I2BURExExcf755w9ahpmNuLKeGVoWw5zE9C+BDwPflwSwBjgg6WoaI/W1LfddA7yZtEgzs3ZpnhnaaVXHqhs43CPiH4ALFr+W9BowERFvSXoc+Lqk+4APAZcA+1Oq1czsLIOcGdotwOvau+8Z7pL2ANcCqyTNAndHxFc73TciDkraC7wEzAM7IuJMivWamQ1sqQAv+opJWelntsyNPW7f0Pb1PcA9ycoyM0vPUgG+2Ls/Pb9Qq969Fw4zs9pbKsDLuqpjUoroOJklVxMTEzE1NVV0GWZWY3U8aCppOiImOt3mkbuZjYQyLsubJS8cZmZWQw53M7MacribmdWQw93MrIYc7mZmNeRwNzOrIYe7mVkNOdzNzGrI4W5mVkMOdzOzGnK4m5nVkMPdzKyGHO5mZjXkcDczq6Ge4S7pIUnHJL3Ysu0/S/pHSX8v6a8lfaDltp2SDks6JOkTWRVuZmbd9TNyfxi4oW3bU8CVEfELwMvATgBJlwPbgCuaj/mKpPHUqjUzs770DPeIeA54u23bkxEx3/xyEljT/Hwr8EhEnIyIV4HDwNUp1mtmZn1Io+f+u8DfND+/CHij5bbZ5rb3kLRd0pSkqePHj6dQhplVzfTMHA88c5jpmbmiS6mdRJfZk/THwDywe3FTh7t1vEhrROwCdkHjGqpJ6jCz6pmemeOmByc5Nb/A8mVj7L5180hdBi9rQ4/cJd0MbAFuinevsj0LrG252xrgzeHLM7M6mp6Z4/5vvczJ0wssBJyeX2DyyImiy6qVoUbukm4A7gB+NSJ+2nLT48DXJd0HfAi4BNifuEozq43WEXvQGGGes2yMzRtXDvw8k0dOsHnjSo/4O+gZ7pL2ANcCqyTNAnfTmB3zM8BTkgAmI+L3I+KgpL3ASzTaNTsi4kxWxZtZ9UweOcGp+caIfUxwzcWruO36SwcKaLd0eusZ7hFxY4fNX13i/vcA9yQpyszqa/PGlSxfNsbp+QXOWTY2cLDD2TuIxZaOw/1siQ6ompkNatP6Fey+dXOilkr7DmLQls4o0LvHQoszMTERU1NTRZdhZhXinjtImo6IiU63eeRuZpW0af2KkQ31fnjhMDOzGnK4m5nVkMPdzAoxzNIDnR7jJQw6c8/dzHI3zDz1To8BPN+9C4/czSx3neapD/OYYZ5nVHjkbma56zRPvdfUxm5z2z3fvTPPczezQrSGOfTXXum0Axjl+e6e525mpdM6T/2BZw73tZxAp7ntnu/emXvuZla4xZbLuIZbIdLeyyN3M+tLlu2PNNabsbM53M2spzyW2HV7JV1uy5hZT55yWD0Od7MaSvusTffEq8dtGbMamZ6ZY9+BWR6deoP5hUitheKeePX0c5m9h2hcCPtYRFzZ3PZB4C+BDcBrwG9HxFzztp3ALcAZ4A8j4puZVG5mZ1nsi5883bg2KaR7lSL3xKuln7bMw8ANbdvuBJ6OiEuAp5tfI+lyYBtwRfMxX5E0nlq1ZtbVYl98MdhFNVsoXggsHf1cQ/U5SRvaNm+lcdFsgK8BzwJ3NLc/EhEngVclHQauBr6TTrlm1k3r6fnj42P81qY1fOqqNbmNttOYKukLX6dn2J77hRFxFCAijkq6oLn9ImCy5X6zzW1mlrFh++JlCmVf+Do9aR9QVYdtHRevkbQd2A6wbt26lMswG02D9sXLFsq+8HV6hg33H0la3Ry1rwaONbfPAmtb7rcGeLPTE0TELmAXNBYOG7IOM0ugbKHsWTnpGTbcHwduBu5tfnysZfvXJd0HfAi4BNiftEgzy0YZQ9mzctLRz1TIPTQOnq6SNAvcTSPU90q6BXgd+DRARByUtBd4CZgHdkTEmYxqN7OEHMr15fXczew9RnmN9Crxeu5m1jdPR6wHry1jZmfxImH14HA3s7N4kbB6cFvGrABl7Wkv1nXXliuY++mp0tVn/XO4m+WsrD3tstZlw3Fbxixnefa0Oy3C1W1hLvfa68Ujd7Oc5XWKfaeRONB1dN5PXWVtJ9l7OdzNcpbXKfbdRuLdlhvoVZfbNtXicDcrwDBncw46au42El9qdL5UXV6xsVoc7mYV0DpqXtbnWu3dRuLD/tXgFRurxeFulkBePejWUfOp+QX2PP86+w7M9myNLN622JJZHJkPU6tXbKwWh7vZkNLoQfe7c1gcNS9eHzXorzXSrcZhd0peHKw6HO5mQ0ragx5k57A4at53YJZHp97gzEL01RrpdlDVB0brz+FuNqSkPehBdw6Lo+Z/d9WavkfdnWr0gdHR4HA3G1LSHvSwO4dBWiPdavSB0frzeu5mBSrqpCCfjFQPXs/drKT6GYVnEcQ+MFp/DnezEvNZoTasRAuHSfq8pIOSXpS0R9L7JH1Q0lOSXml+9G+i1V63xbiS8mJeNqyhw13SRcAfAhMRcSUwDmwD7gSejohLgKebX5uVTlqBvDi6/rMnD3HTg5OpBnzaF87Iaidk5ZO0LbMM+FlJp4FzgTeBncC1zdu/BjwL3JHw+5ilKs12R5ZTC9M8K9QtntEy9Mg9In4I/CnwOnAU+N8R8SRwYUQcbd7nKHBBp8dL2i5pStLU8ePHhy3DRlTSEWia7Y6sL0u3af0Kdlx3ceIgdotntAw9cm/20rcCHwZ+DDwq6TP9Pj4idgG7oDEVctg6bPSkMQJNcxGsvNdcGXb2jBf+Gi1J2jLXA69GxHEASfuAXwZ+JGl1RByVtBo4lkKdZu9Iow2SdiDnNbUwyY7NC3+NliTh/jqwWdK5wD8BHwOmgP8H3Azc2/z4WNIizVqlNQKt4lzvpDu2Kr5mG87Q4R4Rz0v6K+AAMA98j0ab5eeAvZJuobED+HQahZotqsMI1K0Vy5qXH6gQnzLeWdV+LkmPGVTt9Vp2vPxADXga29kWA27Fucv54hMHK/VzcWvF8uBwrwgv0/qu1h3dmMSZhej74hVl4NaK5cHhXhEOhHe17uggGB8TEe9evKLsbYs6HDOw8nO4V0RWgVD2IOykfUd315YrmPvpqXd2eFVoX7m1YllzuFdI2oFQ1T7+Uju6B5457PaVGQ73kVblPn63HZ3bV2YNDveMVKHdUccgrGM/uwq/S1Y+nueegSq1Oxwc5Val3yXL31Lz3BNdrMM6q9Lqe2mtOJimUV5zvP21V+l3ycrFbZkM1LHdkZe8Rqpl/Iul02v375INy+GegTr2fbPSHrJpHeRdKrzL2uro9Np3XHexf5dsKA73jHgec29ZjVR7hXdZZgm174C6vXb/LtkwHO7WUR5ti6xGqr3Cuwytjm47II/SLS0Od3uPvNoWWY1Ue4V3GUK02w7Io3RLi8Odch5cK9KgbYthf35ZhWw/z1t0iJbhrwert5EP97IeXCvSIMGT9OeXVcgWHd69lOGvB6u3kQ/3vEapVTJI8JTl4GQVlX0HZNWWKNwlfQB4ELgSCOB3gUPAXwIbgNeA346I0p6NkucotYy67az6DZ4V5y5nTALC7YU+jcIAwYqXdOT+ZeBvI+K3JC0HzgX+E/B0RNwr6U7gTuCOhN8nM4uj1G8cmEU97lu3UWoal3v74hMHObPQWFP9ri1X9PX41qsoLS7VW+Wf4yDqOECwcho63CWdB/wK8B8AIuIUcErSVuDa5t2+BjxLicN90b4Ds5yaX+AbB2a7vuEGPQhW9hFa0p3V4uMDiAjmfnqq52MWw+3k6cbjxkStQy6rk7TMekkyct8IHAf+u6RfBKaBzwEXRsRRgIg4KumCTg+WtB3YDrBu3boEZSTX7xtukF50FUZoSWdsDPP41h0CUOuQ83ICVqQk4b4MuAr4bEQ8L+nLNFowfYmIXcAuaKwKmaCOxAZ5w/Xbi+53h9HrNPksR/5JZ2wM8/jFn/Wp0wss0Bi5Vznklvo/8nICVqQk4T4LzEbE882v/4pGuP9I0urmqH01cCxpkVnrFVLDhGzrDmN8fIwf/vifmJ6ZO+vxS43u8xr5J52xMejjW3/WVe+59/o/8nICVqShwz0i/pekNyRdFhGHgI8BLzX/3Qzc2/z4WCqVZqzbG27YkF0MsX0HZnl06g0e2f86+9r6+UuN7uvcm61LuPX6P/JcditS0tkynwV2N2fKHAF+h8Ya8Xsl3QK8Dnw64fcoVJKQXTyANr8QHR+/VDvIvdny6+f/qC47MqueROEeES8Ana4C8rEkz1smvd7AvVo2Sz1+qZGdR33ll0U7zywtvsxem05vyG5v0n5bNn6Tj54qzJay6lvqMnsjv/xAq25vyG5/Wg8yhdJv7OLluZOt8zETqwaHe4tB35Dui1dH3iNp/25Y0RzuLQZ9Q7ovXh15j6T9u2FFc7i3GPQM1MX77bju4hyrtGEUMZJ2O86KVPlwb++jJu2r9vOGrPLBslE9uNt+8tTkkRPvbDero0qHe3vI3rXlCr74xMHMQ7eqB8uy2CmVcWex1DLGQGV3zGaDqHS4t4fs37x4NJfQrerBsrR3SmX8C6ZXTVXdMZsNaqzoApJYDNnx5uJTv37lapaNCQHjY8osdBf/xL/945eVItD61f7zSvrz6RSURetVU9o/A7OyqvTIvf0AKADNqwI1Pmb7vYsO9UFbImnP4CjjXzC9avIsFhsVtTpD9YFnDvNnTx5iIWBccPvHL8t1Jkuey/eWpSVSpZ67Wd2MzBmqeYwkh1mKIIsgLkvvuAx/wbQrY01meatVuGf9J/dSIZ338r2bN65k2Zg4fSYyPb5gZtVUq3CHbEdtS4V0Icv35nR8wcyqp3bhvlTbJOmIvkzL904eOcH8mca1SM+cya4t4/61WTXVKty7tU3S6nn3Cuml/mpI+y+KvI4vlOGgrZkNrlbh3q1tkmbPuywH6/KY0tf6czs1v8D933qZ266/tBSv38yWlvgkJknjkr4n6Ynm1x+U9JSkV5ofc0uCbieo1PXElU3rV7Djuov7CtvpmTkeeOYw0zNzfT//4s9tTLAQ8O1X3uKmBycHeg4zK0biee6Sbqdxqb3zImKLpC8Bb0fEvZLuBFZExB1LPUeaV2LKsudeVUnaK9Mzc9z/rZf59itvERRz/oCZdZbZPHdJa4DfBO4Bbm9u3gpc2/z8a8CzwJLhnqZubZM82ill3YEkvcj3bddfyndfe7tUZ6Ka2dKS9tzvB74A/HzLtgsj4ihARByVdEGnB0raDmwHWLduXcIyijc9M8eN/23ynQDc83vlOfiY9OCrT9k3q56hw13SFuBYRExLunbQx0fELmAXNNoyw9ZRFvsOzHJqfgFoHHzcd2C2NCGYRjiX5UCymfUnycj9GuCTkn4DeB9wnqS/AH4kaXVz1L4aOJZGoWXXvncq297K4Ww2WoaeLRMROyNiTURsALYBfxcRnwEeB25u3u1m4LHEVVbAp65aw/LxxnLDy8fFp65aU3RJZjbCspjnfi+wV9ItwOvApzP4HqWzaf0K9mz/1wO1Psp6ANbMqi+VcI+IZ2nMiiEiTgAfS+N5y6ZXGA/S+vDZn2aWpVqdoZqltMO4LEv2mlk9Vfoye3lK+5JydT1r1szKwSP3PqW9UJfnjptZlmp1mb2s+QComZXJyFxmL2ueK25mVeGeu5lZDTnczcxqyOFuZlZDDnczsxpyuJuZ1ZDD3cyshkoxz13ScWBmiIeuAt5KuZwq8OseLX7do2WQ170+Is7vdEMpwn1Ykqa6TeCvM7/u0eLXPVrSet1uy5iZ1ZDD3cyshqoe7ruKLqAgft2jxa97tKTyuivdczczs86qPnI3M7MOHO5mZjVU2XCXdIOkQ5IOS7qz6HryIOkhScckvVh0LXmRtFbSM5J+IOmgpM8VXVMeJL1P0n5J32++7j8puqY8SRqX9D1JTxRdS14kvSbpHyS9ICnxBS4q2XOXNA68DPwaMAt8F7gxIl4qtLCMSfoV4CfAn0fElUXXkwdJq4HVEXFA0s8D08C/HYH/awHvj4ifSDoH+DbwuYiYLLi0XEi6HZgAzouILUXXkwdJrwETEZHKiVtVHblfDRyOiCMRcQp4BNhacE2Zi4jngLeLriNPEXE0Ig40P/+/wA+Ai4qtKnvR8JPml+c0/1VvJDYESWuA3wQeLLqWKqtquF8EvNHy9Swj8IYfdZI2AL8EPF9sJflotiZeAI4BT0XESLxu4H7gC8BC0YXkLIAnJU1L2p70yaoa7uqwbSRGNaNK0s8B3wBui4j/U3Q9eYiIMxHxEWANcLWk2rfiJG0BjkXEdNG1FOCaiLgK+HVgR7MNO7SqhvsssLbl6zXAmwXVYhlr9py/AeyOiH1F15O3iPgx8CxwQ8Gl5OEa4JPN/vMjwEcl/UWxJeUjIt5sfjwG/DWN9vPQqhru3wUukfRhScuBbcDjBddkGWgeWPwq8IOIuK/oevIi6XxJH2h+/rPA9cA/FltV9iJiZ0SsiYgNNN7XfxcRnym4rMxJen9zwgCS3g98HEg0K66S4R4R88B/BL5J4wDb3og4WGxV2ZO0B/gOcJmkWUm3FF1TDq4B/j2NEdwLzX+/UXRROVgNPCPp72kMZp6KiJGZFjiCLgS+Len7wH7gf0bE3yZ5wkpOhTQzs6VVcuRuZmZLc7ibmdWQw93MrIYc7mZmNeRwNzOrIYe7mVkNOdzNzGro/wOUzQnZ3afo+QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "df = pd.read_csv('data1.csv')\n", "plt.plot(df.col1.values, df.col2.values, '.'); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fijamos mediante mínimos cuadrados un modelo lineal tomando como variable explicativa a col1, y de salida a col2:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "X, y = df.col1.values.reshape(-1, 1), df.col2.values\n", "fitted_mod = lm().fit(X, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si pintamos la variable de salida real (col2) vs. la predicción del modelo:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAewklEQVR4nO3df5wcdZ3n8ddnZgiCCxJIYCOTZIj82sCpy4xczh8rKiC6SPYBygZFAWFzejkV2Tsg7h4/3MVjxUO9E9dHFgO4QEKQLHByKrIGOE8HnI5hJSgYQyYMyZoQBlZRk8zM5/7obuh0qvpHdVV1dfX7+XjwyExVd/Wnw8y7v/nU91tl7o6IiORLT7sLEBGR+CncRURySOEuIpJDCncRkRxSuIuI5FBfuwsAmDFjhg8MDLS7DBGRjlIoFJ5z95lB+zIR7gMDA4yMjLS7DBGRjmJmo2H71JYREckhhbuISA4p3EVEckjhLiKSQwp3EZEcUriLiOSQwl1EJEWF0XFuWLOBwuh4oq+TiXnuIiLdoDA6zoduHGbXxBTT+nq47aIFDM6dnshraeQuIpKS4Y072DUxxZTD7okphjfuSOy1FO4iIilZMO8QpvX10GuwT18PC+YdkthrqS0jIpKSwbnTue2iBQxv3MGCeYck1pIBhbuISKoG505PNNTL1JYREcmhuuFuZsvNbJuZPV6x7Y1mNmxm68xsxMxOrNi31Mw2mNmTZvbupAoXEZFwjYzcbwZOq9r2eeBqd38jcEXpe8xsPrAIOK70nK+aWW9s1YqISEPqhru7Pww8X70ZOLD09WuALaWvFwIr3X2nuz8NbABOREREUhX1hOrFwHfN7AsUPyDeXNp+ODBc8bix0ra9mNliYDHAnDlzIpYhIiJBop5Q/TjwaXefDXwa+HppuwU81oMO4O7L3H3I3Ydmzgy8S5SIiEQUNdzPA1aXvr6TV1ovY8Dsisf180rLRkREUhI13LcAby99/U7gF6Wv7wUWmdm+ZnYEcBTwaGsliohIs+r23M1sBXASMMPMxoArgb8AvmxmfcDvKfXO3X29ma0CngAmgCXuPplQ7SIiEqJuuLv7OSG7BkMefw1wTStFiYhIa7RCVUQkhxTuIiI5pHAXEckhhbuISA4p3EVEckjhLiKSQwp3EZEcUriLiOSQwl1EJIcU7iIiOaRwFxFJSWF0nBvWbKAwOp74a0W9WYeIiDShMDrOh24cZtfEFNP6erjtogUMzp2e2Osp3EVEUjC8cQe7JqaYctg9McXwxh0vb18w75DYg17hLiKSggXzDmFaXw+7J6bYp6+H6ftPS3Qkr3AXEUnB4Nzp3HbRgpdH6kEjeYW7iEgHGpw7fY8ArxzJL5h3SKyvVXe2jJktN7NtZvZ41fZPmNmTZrbezD5fsX2pmW0o7Xt3rNWKiGRQlFkw5ZH8Jacek8jJ1UZG7jcDXwG+Ud5gZu8AFgKvd/edZnZoaft8YBFwHPBa4AEzO1q32hORvGplFkz1SD5OdUfu7v4w8HzV5o8D17r7ztJjtpW2LwRWuvtOd38a2ACcGGO9IiKZEjYLpt2iLmI6GnibmT1iZg+Z2ZtK2w8Hnql43Fhp217MbLGZjZjZyPbt2yOWISKSnqD2S3kWTK9Rt3feCYuY+oDpwALgTcAqM5sHWMBjPegA7r4MWAYwNDQU+BgRkawIa79Uz4IJa7N0yiKmMWC1uzvwqJlNATNK22dXPK4f2NJaiSIi7Vdr6mIjvfOkpz5Wi9qWuRt4J4CZHQ1MA54D7gUWmdm+ZnYEcBTwaByFioi0UzPtlySe36y6I3czWwGcBMwwszHgSmA5sLw0PXIXcF5pFL/ezFYBTwATwBLNlBGRPGi0/ZLU85tlxUxur6GhIR8ZGWl3GSIiHcXMCu4+FLRPl/wVEckhhbuISJPSnNIYla4tIyLShLSnNEalkbuISBOyuiK1msJdRKQJWV2RWk1tGRGRJmR1RWo1hbuISJOyuCK1mtoyIiIJSHtFajWN3EVEWlQYHd+rTZP2itRqCncRkRbU6q0neTOOetSWERFpQVanRircRURa0O7eehi1ZUREWtDu3noYhbuISIva2VsPo7aMiEgOKdxFRHJI4S4ikkN1w93MlpvZttIt9ar3/RczczObUbFtqZltMLMnzezdcRcsIiL1NTJyvxk4rXqjmc0GTgE2V2ybDywCjis956tm1htLpSIiVTrhphntUne2jLs/bGYDAbu+CFwK3FOxbSGw0t13Ak+b2QbgROBHrZcqIvKKdl91Mesi9dzN7AzgWXd/rGrX4cAzFd+PlbYFHWOxmY2Y2cj27dujlCEiXSyrK0OzoulwN7P9gb8CrgjaHbDNg47j7svcfcjdh2bOnNlsGSLS5bK6MjQroixieh1wBPCYmQH0A2vN7ESKI/XZFY/tB7a0WqSISLVmV4YGXbmxkX2dqulwd/efAoeWvzezTcCQuz9nZvcCt5vZ9cBrgaOAR2OqVURkD42uDK3Vn89r776RqZArKJ4QPcbMxszswrDHuvt6YBXwBPAdYIm7T8ZVrIhIFLX683nt3TcyW+acOvsHqr6/BrimtbJEROJT7s/vnpjaqz9fa1/iXtgMr54J++wX+6HNPfB8Z6qGhoZ8ZGSk3WWISI5lpufuDg9fB2tKY+C3XAynXB3pUGZWcPehoH26KqSIdIVa/flUruq4+RFYfure2990USIvp3AXEUnS8tNgc8A6zkufhv0PTuxlFe4iIkm46jXB2//bc9C7T+Ivr3AXEYnLcxvgK4PB+656MdVSFO4iIq36H8fCr7fuvX3Bf4LT/nv69aBwFxGJLqz18pdPwQGHpVtLFYW7iEgzXnoOrntd8L6UWy+1KNxFRBrxtbfBv/5L8L4MhXqZwl1EpJaw1ssHV8HR2b3ZnMJdRKTa1CR8NmQO+pUvgAVd3TxbFO4iImUPXQdr/jZ4XwZbL7Uo3EVEwlov8/8Mzr4l3VpionAXke4VFupLx2DfA9KtJWYKdxHpLqM/hJveE7yvw1ovtSjcRaQ7hI3SIVehXlY33M1sOXA6sM3djy9tuw54H7AL+CVwgbu/UNq3FLgQmAQ+6e7fTah2EZH6wkL9Y/8P/vD4dGtJUSMj95uBrwDfqNj2PWCpu0+Y2d8BS4HLzGw+sAg4juI9VB8ws6N1qz0RSdVvtsMXjgzel8NRepBGbrP3sJkNVG27v+LbYeD9pa8XAivdfSfwtJltAE6keA9WEZE9xH4HpBtPgbFHg/d1SaiXxdFz/yhwR+nrwymGfdlYadtezGwxsBhgzpw5MZQhIp2kMDrOh24cZtfEFNP6erjtogXRAz6s9fK+L8Pg+ZFr7GQthbuZ/RUwAdxW3hTwsMCbtLr7MmAZFO+h2kodItJ5hjfuYNfEFFMOuyemGN64o7lwz8Eq0iRFDnczO4/iidZ3+St32R4DZlc8rB/YEr08Ecmjwug4W174HX09xuSUs09fDwvmHdLYk+/9BKz9RvC+Lmu91BIp3M3sNOAy4O3u/tuKXfcCt5vZ9RRPqB4FhDTARKQbVbZj+np7+PMTZ3PWCf31R+0hrZd1fiSTH30g+Rtcd5hGpkKuAE4CZpjZGHAlxdkx+wLfs+I/fYbd/WPuvt7MVgFPUGzXLNFMGRGpVNmOmZyc4vCD9qsdzCGhPrjza+zwA+k1uKTZlk4XaGS2zDkBm79e4/HXANe0UpSI5NeCeYcwra+H3RNT4e2YJ+6BVR8JPsBVL1IYHeelG4fprXWMLmevtMvbZ2hoyEdGRtpdhoikJHQKZBOrSGOfRtmBzKzg7kNB+3T5ARFJ3eDc6Y2F+oe+CUed0tgxZA8KdxFpjxc2w5f+XfA+zXppmcJdRNLVZRfwaheFu4ikIyzUjz8L3r888mHVew+mcBeR5ExNwWeDA7dw/kaGnx4vhnKDh6sO8lgvYZAzCncRid+tZ8GGB4L3laYyNhvKQc9p+RIGOaZwF5H4NNhPDwrl8vaw9krQcxqaM9+lFO4i0rqwUP/0enhN/16bq0N5+v7T6o7kg4J8cO70l0fw6rnvSeEuItGsWwF3fyx4X51ZL9Wh3Eh7JSzINd89mMJdRJoT01TG6lBupL2iIG+cwl1EGhMW6md8BU74cEuHVnslfgp3EQk3Pgpffn3wvpgXHGlUHi+Fu4jsrUbr5djJlcUTnimWI81TuIvkUORVmyGh/ut9D+MN//ZFphx6TfPJO4HCXSRHCqPjrF47xp0jzzAx5Y0tEKqxipQrnoeeXp4aHWfajcOaT95BGrkT03KK90rd5u7Hl7YdDNwBDACbgLPdfby0bylwITAJfNLdv5tI5SKyh/IKzp27p16+K33NVZvXHQUvbQs+WFU/XSc8O08jI/ebga8AlXekvRz4Z3e/1swuL31/mZnNBxYBx1G8h+oDZna0brUnkrzyXPFysBsEj7IjTmVM64SnLgQWj0Zus/ewmQ1UbV5I8b6qALcAD1K8YfZCYKW77wSeNrMNwInAj+IpV0TCVK7g7O3t4f2D/XveeDos1D/+IzhsfsuvH0co60Jg8Ynacz/M3bcCuPtWMzu0tP1wYLjicWOlbSKSsMDWyYN/Bzd9LvgJpVF6YXSc4TUbMhHKuhBYfOI+oWoB2wJv0mpmi4HFAHPmzIm5DJHu9HLrpMHWS9ZCWRcCi0/UcP+Vmc0qjdpnAeWzMmPA7IrH9QNbgg7g7suAZVC8QXbEOkSkUlioD30UTv/iXpuzFso6cRufqOF+L3AecG3pz3sqtt9uZtdTPKF6FPBoq0WKSA3bn4Ib3hS8r84q0iyGslaqxqORqZArKJ48nWFmY8CVFEN9lZldCGwGPgDg7uvNbBXwBDABLNFMGZGExHABr7BQjnJyVKGcLebe/o7I0NCQj4yMtLsMkc6Q8A2mNWOlc5hZwd2HgvZphapIJ3CHqw8K3vfX26Bv39heSjNW8kHhLpJlCY/Sg2jGSj4o3EWyqA2hDq/02q84/TjGf7tLM1Y6mMJdpA1CT1iGhfr598HAWxOvSb32/FC4i6SsOkS/P/hDXrvufwY/OMFRejX12vNF4S6SsnKIbtz3g8UN6wIeFFOoB/0LIexfDeq154umQoqkLaz1Mu8d8JG7Y3uZoDYLULP1Um9+u67YmC2aCinSbjt+Cf/rhOB9CbVegtosQM3WS62FSOrJdxaFu0iSYpz10uyoOazNErX1op58Z1G4i7Sg6VkvEGmkXjlq7gu6VnuA8qUF7lo79vLlWlu5Box68p1F4S4S0V5tigv/PYM3HxH84M9shWn7R36tylHzrokpVjyymdVrxxpqjaxeO8auiSnuqnp8uU3TzLVjdMXGzqFwF4lor1kvNwc8qM4ovdFWS3nUXL4/qtNYaySs7x61d66Lg3UOhbtIREseGmRJ2CVdGmi9NHOCsjxqXr12jDtHnmFyyhtqjQS1UtQ77w4Kd5FmhfXTF62AY9/b8GGaDdnyqPnME/obbo2EtVLUO88/hbtIIx7+Anz/b4L3RZzKGPUEZbOtkerHq3feHbSISaSWFK6dXi9ktXBIwmgRk0izwkJ95h/BkuHYXqbeKFwLhySqlsLdzD4NXETx5P1PgQuA/YE7gAFgE3C2u4+3VKVIGl4cgy8eF7wvplkvzdLJT4kqcrib2eHAJ4H57v670r1TFwHzgX9292vN7HLgcuCyWKoViVE5kJc8NBj+oJhnvTQr7oVDavF0j1bbMn3Afma2m+KIfQuwlOINtQFuAR5E4S4ZUxgdZ/CmAUJjvYl+epKj6zhPfqrF010ih7u7P2tmXwA2A78D7nf3+83sMHffWnrMVjM7NOj5ZrYYWAwwZ86cqGVIl2ppBHrVa4JD/fJn4FUHNl1L0svy41o4pBZPd2mlLTMdWAgcAbwA3Glm5zb6fHdfBiyD4myZqHVI94k0Ar1mFuz+beCuYydXFo8RIdgh/amFUT/YdG2Y7tJKW+Zk4Gl33w5gZquBNwO/MrNZpVH7LGBbDHWKvKypEWiNqYyFCzYxvHEHt8UQyGkty2+ltaL57d2llXDfDCwws/0ptmXeBYwALwHnAdeW/ryn1SJFKjU0Ag0L9bP/EeafAcAgjV80Kytaba3o2jDdo5We+yNm9k1gLTAB/IRim+UPgFVmdiHFD4APxFGoSFnoCHTd7XD3x4OflOK9SBuh1ookTStUO4imsYVIeBVp3FqdtaKfAynTCtUc0DS2PZWnMgaa82b46LdTracZaq1IGnraXYA0Juy63F3npeeKUxkDgv11O2/nhrcXMh3s8EprpddQa0USo5F7h+j6XmuN1suRu1bg/sr1zbPettCsFUmDeu4dJInQynoQ1gr11+28nX36erji9OMY/+2ulz/w1L6SbqGee07E3WvNdB8/LNQvG4X9DqIwOs4lAR9KN6zZoFWYIijcu1rmlqP/wzvh2ULwvqpZL2EfdF3fvhIpUbgnJPPtDjIUhDFOZcxjP7sTfpYke9RzT0Cm2x1V2hocYaH+57fCH70v3VoyqpN+liR96rmnLHPtjhpSnzP98/8DK88J3lcapRdGxxles6ErR6rVH7ad9LMk2aJwT0Bm2h1Z0mDrJa2RahZbHUHvXT9LEpXCPQF57PtGFhbqA2+D87+11yg9rpFqrfDOSqujkVH6knccqZ8liUThnpCuXiL++xfh2pAbsNQZpccxUq0X3llodTTz3rv6Z0kiU7hLoEhtiyZnvSQ1Uq0X3llodWiULklTuJPN/ms7Nd22iDiVMamRar3wzkLbTKN0SVrXh3tW+q9Z0nDbIizUL30a9j+47uskFbKNHLfdIZqFDxjJt64P9yz0X7Om5sh31UfgieCba71u5+1ccuoxLGkg2MuSCtl2h3cjOqFG6VwthbuZHQTcCBwPOPBR4EngDmAA2ASc7e7jLVWZoGb7r3lr4QS9n8BRZY3Wy7GTKzVVTyRjWlqhama3AP/X3W80s2nA/sBngOfd/VozuxyY7u6X1TpOu1eoNhrYeWvhNPR+6tyLtDA6zl1rxzDgzBP6dUehBnTr+5b4JbJC1cwOBP4EOB/A3XcBu8xsIXBS6WG3AA8CNcO93cq/YOUbYIT9wuWthRP6fp66H24PufVtjamMZ57Q39Drlj8QvlkYY2IyHx+UjcrbAEGyq5W2zDxgO3CTmb0BKACfAg5z960A7r7VzA4NerKZLQYWA8yZEzInOiWN/sLlrYVT/X6WPDQID4U8uMGpjPXeZ/nveufuKcr/ZszDB2UYXU5A2qWVcO8DTgA+4e6PmNmXgcsbfbK7LwOWQbEt00IdLWv0F66ZGQ6dMEIrv5/Qe5HOegP8x4dDnx9lvnj577r8P9zo7FvNNbsSNgtz7KU7tBLuY8CYuz9S+v6bFMP9V2Y2qzRqnwVsa7XIpNX7hav+BW4kpBv9wKgXDomN/He9BJ97LYNB+xq8zG6U6XyVf9e9PcYHhmY33avPiigrYbVQSdISOdzd/V/N7BkzO8bdnwTeBTxR+u884NrSn8Hz5jKkVkhFHYHvEWK9PTz7wu8ojI43fOzERv4xXjsdmp/Ol6f53VFXwmoKpKSh1XnunwBuK82U2QhcAPQAq8zsQmAzEHJmLlvCfuGi9kjLIbZ67Rh3jjzDykc3s3rt2B4hXevYsfdmYw71VuQl3DphJax0r5bC3d3XAUHTcN7VynGzpNmWTaXyCbSJKQ8M6VrHjq03Gxbq//WX8OoZ0Y4pQP3wzvoJdck33YmpAWG/pI20TsqPKYd09WMS6bn/74uhcFPwvpRH6d2qE06oS+fTnZiaELZiM2rLpt7orlaLoun2RYZaL1mU5khaUx6l3RTuFZodbTXaOkm8xxwW6u+/CY4/M7nX7SBpj6Q15VHaTeFeodnRVltPmG0ehuXvDt6nUfpe0h5J62SqtJvCvUKU0VbqMz/UeomkHSPpvMwKks7U8SdUq/uorfZVm7mIWKqjsrBQ738TXPRAw4fp5hkc3fzeJZ9qnVDt6HCv7qNecfpxfPZb6xPvq6bWv939e7jmsOB9EUbpSdSdxcCsV1MWaxaJIrezZar7qN9+fGsqfdXE+7cJtV7irjuL0/3q1ZTFmkWS0NPuAlpR7qP2WvHiU+85ftYe3yfVV61+3dhe56rXhAf7VS+23FOPu+6gD4t2q1dTFmsWSUJHj9zDZiR8+/GtvOf4WYmNyGKfCREW6H/5JBzwh6FPa7a9EHfdWZzuV6+mLNYskoSO7rlXa/c/uZtabfrA1fCD64MP1MAIvd3vtbKOrPWv1XOXbpHbnnu1NOYyR7kUQeW+jft+MPzgTbRdsrICMovT/erVlMWaReKWq3BP+p/ctQK83hUef967CHoDDnrmjfD65i+cuWDeIfT1GLsnnd4eU3tBRPaQq3BPelVgrQAP/GDZ9jP46gKWBB0sjgVHZoCX/hQReUWuwr2WOPqstf5lUPnBsuShQQi5KGPhgk2xfOgMb9zBxGTxdnWTk8m1ZdSfFulMuQr3sLZJXCcf617h8aaB4NvWHf0e+ODK4mMivK8gacz6yMpJWxFpXq7CPaxtEufJx71Oxk1NwmcPDn7wlS8k1jJJ48JUlX9vuyam+NIDT3HxyUcr4EU6QMvhbma9wAjwrLufbmYHA3cAA8Am4Gx3H2/1dRoRNppNZJT7nc/A8A3B+1K6gFczsz6itFfKf2/lgP/BL57jx5ue1whepAO0PM/dzC6heKu9A0vh/nngeXe/1swuB6a7+2W1jhHnnZhqTVWMZZQbtuBo+gB86rHox01QK+2Vwug4X3rgKX7wi+dwoNfgklOPYck7jky2aBGpK7F57mbWD/wpcA1wSWnzQuCk0te3AA8CNcM9TmGj2ZbnNoeF+mWjsN9BQHZPPrbSlhqcO52LTz6aH296Xqs6RTpIq22ZLwGXAgdUbDvM3bcCuPtWMzs06IlmthhYDDBnzpwWy0jIlnWw7O3B+6paL4XRcc75h1fulbriL7LTumi1LaUbT4h0nsjhbmanA9vcvWBmJzX7fHdfBiyDYlsmah2J+NvDYOL3e2+f+1a44L7Ap6xeO8auiSmgePJx9dqxzIRgHOGsVZ0inaWVkftbgDPM7L3Aq4ADzexW4FdmNqs0ap8FbIuj0FSEtV6W/BhmHl3zqdWfTtn6tFI4i3SbyJf8dfel7t7v7gPAIuD77n4ucC9wXulh5wH3tFxlknb9NvxSu+XL7NYJdoCzTuhnWq9hwLRe46wT+us+pzA6zg1rNlAYTWUykYh0kSTmuV8LrDKzC4HNQPMXTknDT26FewIvDBBpKuPg3OmsWPwfGm59aIGQiCQplnB39wcpzorB3XcA74rjuIkIa72cfx8MvLXmU+vNhmmm9ZGVqzqKSD7laoVqKHe4+qDgfQ2uIo17pK2bRohIkvId7r9aD3//5uB9TbZe4h5pa3qhiCQpn+F+9xJYd+ve28+9C448OdIhkxhpawaLiCQlX+Eetujor7dB374tHVojbRHpJJ0f7pMT8OR9MPw12PzDPffFfAEvjbRFpFN0drg/W4BV58OLm+GgOXDqNfDH5758rRcRkW7V2eF+8DyYcSSc9jk45r3QE3STUhGR7tPZ4b7fdPjwP7W7ChGRzIl8+QEREckuhbuISA4p3EVEckjhLiKSQwp3EZEcUriLiOSQwl1EJIcU7iIiOWTu7b/bp5ltB0YjPHUG8FzM5XQCve/uovfdXZp533PdfWbQjkyEe1RmNuLuQ+2uI216391F77u7xPW+1ZYREckhhbuISA51ergva3cBbaL33V30vrtLLO+7o3vuIiISrNNH7iIiEkDhLiKSQx0b7mZ2mpk9aWYbzOzydteTBjNbbmbbzOzxdteSFjObbWZrzOxnZrbezD7V7prSYGavMrNHzeyx0vu+ut01pcnMes3sJ2b2rXbXkhYz22RmPzWzdWY20vLxOrHnbma9wFPAKcAY8GPgHHd/oq2FJczM/gT4DfANdz++3fWkwcxmAbPcfa2ZHQAUgD/rgv/XBrza3X9jZvsAPwA+5e7DbS4tFWZ2CTAEHOjup7e7njSY2SZgyN1jWbjVqSP3E4EN7r7R3XcBK4GFba4pce7+MPB8u+tIk7tvdfe1pa9/DfwMOLy9VSXPi35T+naf0n+dNxKLwMz6gT8Fbmx3LZ2sU8P9cOCZiu/H6IJf+G5nZgPAHwOPtLeSdJRaE+uAbcD33L0r3jfwJeBSYKrdhaTMgfvNrGBmi1s9WKeGuwVs64pRTbcysz8A7gIudvd/a3c9aXD3SXd/I9APnGhmuW/FmdnpwDZ3L7S7ljZ4i7ufALwHWFJqw0bWqeE+Bsyu+L4f2NKmWiRhpZ7zXcBt7r663fWkzd1fAB4ETmtzKWl4C3BGqf+8Eninmd3a3pLS4e5bSn9uA/6JYvs5sk4N9x8DR5nZEWY2DVgE3NvmmiQBpROLXwd+5u7Xt7uetJjZTDM7qPT1fsDJwM/bW1Xy3H2pu/e7+wDF3+vvu/u5bS4rcWb26tKEAczs1cCpQEuz4joy3N19AvjPwHcpnmBb5e7r21tV8sxsBfAj4BgzGzOzC9tdUwreAnyY4ghuXem/97a7qBTMAtaY2b9QHMx8z927ZlpgFzoM+IGZPQY8Ctzn7t9p5YAdORVSRERq68iRu4iI1KZwFxHJIYW7iEgOKdxFRHJI4S4ikkMKdxGRHFK4i4jk0P8HKNU6sW7AJtMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(df.col1.values, df.col2.values, '.')\n", "plt.plot(df.col1.values, fitted_mod.predict(X)); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "mientras que el RMSE es:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "16.554928774453874" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rmse(df.col2.values, fitted_mod.predict(X))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sin embargo, si fijamos otro tipo de regresión (en este ejemplo se fija una regresión polinómica), el error RMSE es menor que el RMSE anterior, que era el menor error posible con regresión lineal:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZxU5ZX/8c/pDUUUERCBZhFFIxBFGhHXoGI0xojGXYw4YphEE4Mm45KMZhsTk4yO+UUyE6LEDSWgGI2KOy4xtobGBUGFFmxoQNkaUQR6qfP741ZD0VRvVbfW/r5fr351961bVaeKrnMP5z73eczdERGR/FKQ6QBERCR8Su4iInlIyV1EJA8puYuI5CEldxGRPFSU6QAAevTo4QMHDsx0GCIiOaWiomKdu/eMd1tWJPeBAwcyb968TIchIpJTzKyqudvUlhERyUNK7iIieUjJXUQkDym5i4jkISV3EZE8pOQuIpKHlNxFRNKooqqGKXMrqaiqSenzZMU4dxGRjqCiqobxd5ZTWx+hpKiA6ZePpmxAt5Q8lyp3EZE0KV+6ntr6CBGHuvoI5UvXp+y5lNxFRNJk9KDulBQVUGhQXFTA6EHdU/ZcasuIiKRJ2YBuTL98NOVL1zN6UPeUtWRAyV1EJK3KBnRLaVJvpLaMiEgeajW5m9k0M1tjZu/GbBtuZuVm9paZzTOzUTG33WBmlWb2gZmdkqrARUSkeW2p3O8GTm2y7bfAz919OHBT9HfMbAhwATA0ep8/mllhaNGKiEibtJrc3f1lYEPTzcBe0Z+7AquiP48DZrj7NndfBlQCoxARkbRK9ITqZOBpM/tvggPE0dHtfYHymP2qo9t2YWaTgEkA/fv3TzAMERGJJ9ETqt8Frnb3fsDVwF3R7RZnX4/3AO4+1d1HuvvInj3jrhIlIiIJSjS5TwBmR3+exY7WSzXQL2a/Una0bEREJE0STe6rgK9Efz4RWBL9+THgAjPrZGb7A4OBN5ILUURE2qvVnruZPQiMAXqYWTXwU+DbwO/NrAjYSrR37u4LzWwmsAioB65094YUxS4iIs1oNbm7+4XN3FTWzP43AzcnE5SIiCRHV6iKiOQhJXcRkTyk5C4ikoeU3EVE8pCSu4hIHlJyFxHJQ0ruIiJ5SMldRCQPKbmLiOQhJXcRkTyk5C4ikoeU3EVE0qSiqoYpcyupqKpJ+XMluhKTiIi0Q0VVDePvLKe2PkJJUQHTLx9N2YBuKXs+JXcRkTQoX7qe2voIEYe6+gjlS9dTuG0j/6yu5cgD9g090astIyKSBqMHdaekqIBCg+KiAnqU1FM8/Sz6zL2a8XeWh96qUeUuIpIGZQO6Mf3y0ZQvXc/ogXvT7fHLGMBH3NpwNnUeVPJhVu9K7iIiaVI2oFuQwJ/+Cax/iZt9Ai/5CIqLChg9qHuoz9WWZfamAacDa9x9WMz27wPfI1hO7wl3vza6/QZgItAAXOXuT4casYhIlqmoqgkq8kHdW6++502D1+6AUZM4deiP2but92untlTudwN3APc2bjCzE4BxwKHuvs3M9o1uHwJcAAwF+gDPmdlBWkdVRPJVu0bBVD4HT/wIBn8VTvk1ZYVFKRsx0+oJVXd/GdjQZPN3gVvcfVt0nzXR7eOAGe6+zd2XAZXAqBDjFRHJKvFGwcT18QIaZlzCus4H8OaRt0FharviiY6WOQg4zsxeN7OXzOyI6Pa+wIqY/aqj23ZhZpPMbJ6ZzVu7dm2CYYiIpE+8i5CajoKJ2zv/dCW1957DmrpOfGPDVVx4z7spv5Ap0UNHEdANGA0cAcw0s0GAxdnX4z2Au08FpgKMHDky7j4iItmiufbLTqNg4vXOt26CB86DbZ8xse5GVvs+FNaHPzqmqUSTezUw290deMPMIkCP6PZ+MfuVAquSC1FEJPPitV8ak/P2UTBNNdTBrAmw9n0+GvsXls4pptAjKRkd01Siyf1vwInAi2Z2EFACrAMeAx4ws9sITqgOBt4II1ARkUxqbL/U1bcxObvD41fDhy/AGXdw0IgzmN63HaNqktSWoZAPAmOAHmZWDfwUmAZMM7N3gVpgQrSKX2hmM4FFBEMkr9RIGRHJB622X5p68RZ48z44/loY8a3tj5HqpN7IgpycWSNHjvR58+ZlOgwRkXBU3AN/vwqGXwzj7gCLdzoyeWZW4e4j492muWVERMK0+OmgHXPgWPjG7SlL7K1RchcRaadm52VfWQGzLoX9hsG590BhcUbiA80tIyLSLs1ekbr+Q5h+HuzRAy6aBZ26ZDROVe4iIu0Q94rUzetg+jngEbh4NuzZK9NhKrmLiLRH0ytSj+63G0w/FzatggtnQI/B2/dN57J6TaktIyLSDrFDIo8asCeHv/odWP02nH8/9D9y+37pXlavKSV3EZF2KhvQjbJ+XeGRSdsvUuJLp+20T0tXtKaD2jIiIu3lDs/8BBbMgpN+uv0ipVhtmlAshVS5i4i016u3Q/kf4cjvwrFXx12so91XtIZMyV1EpD3m3wfP/QyGnQOn/IqK5Rub7a2nc7qBptSWERFpq/f+HkwrcMCJcOb/QkFB2xfrSDMldxGRtqh8Hh66DPqOhPPug6ISIPO99eaoLSMi0prl5TBjPPQ4GMbvfPVppnvrzVFyFxFpyeq3g4uUuvaFbz0Cu++9yy6Z7K03R20ZEZHmrF0M950Fu3WFSx6FLj0zHVGbKbmLiMRTUwX3jgMrDBJ719JMR9QuasuIiDT12cdBYq/bDJc+Cd0PyHRE7dZq5W5m08xsTXRJvaa3/cjM3Mx6xGy7wcwqzewDMzsl7IBFRFLqiw1BK+bzNTD+4WBu9hzUlrbM3cCpTTeaWT/gZGB5zLYhwAXA0Oh9/mhmhaFEKiLSROizLm77DO4/O5ib/cIHod8R4TxuBrTalnH3l81sYJyb/ge4Fng0Zts4YIa7bwOWmVklMAp4LflQRUR2CH3Wxbot8MAFO2Z4HPSV8ILNgIROqJrZGcBKd3+7yU19gRUxv1dHt8V7jElmNs/M5q1duzaRMESkAwv1ytD6Wpg5AapehbP+tMsMj7mo3cndzDoDPwFuindznG0e73Hcfaq7j3T3kT175s7wIhHJDqFdGdpQD4/8Oyx5Gk6/DQ49N9xAMySR0TIHAPsDb1uwqncpMN/MRhFU6v1i9i0FViUbpIhIU6FcGRppgL99FxbO5tVBP2C3nmdRFn6oGdHu5O7uC4B9G383s4+Ake6+zsweAx4ws9uAPsBg4I2QYhUR2Ul7rgzdZVreSAQe+z4smMn/RM7nD+8dScmS8rSvmJQqrSZ3M3sQGAP0MLNq4Kfufle8fd19oZnNBBYB9cCV7t4QYrwiIu22y8nXiaMoe+fn8NZ03hjw7/xh8VcytmJSqrRltMyFrdw+sMnvNwM3JxeWiEh4dj752kDJM9fBqllw3A8pPOB7lCx9nbr6SFbN6pgsXaEqInmv8eRrXX0DNxXfz5dXzYGjr4ITb6TMLCtndUyWkruI5L2yAd2YPvFICp67icOr5wTL4538CwgGhWTlrI7JUnIXkfznTtmS30P1fXDE5XDqr7cn9nyl5C4i+c0dnvlPeO0OKPs3+Nrv8j6xg5K7iOQzd3jqenj9/2DUJPjabztEYgcldxHJV5EIzPkP+NedMPoKOOVXHSaxg5K7iOSjSASeuBoq7g5GxcScPO0olNxFJL9EGuDvV8Gb98Ox18BJN3W4xA5K7iKSTyIN8OiV8PaD8JXrYMwNHTKxg5K7iOSLhnr423dgwSwY82MYc12mI8ooJXcRyX0NdTD727DwkaANc9wPMx1Rxim5i0hua6iDhy6D9x6Dk38Jx1yV6YiygpK7iOSuuq0w61JYPAdO+TUcdUWmI8oaSu4ikptqN8OMi2Dpi/D1W4NpBWQ7JXcRyT1bP4Xp50H1G3Dm/8LwizIdUdZRcheR3LJ5Hdx3Fqx5D86ZBkPPynREWUnJXURyx6cr4b4zYeNyuPBBGHxypiPKWgWt7WBm08xsjZm9G7Ptd2b2vpm9Y2aPmNneMbfdYGaVZvaBmZ2SqsBFpINZVwnTToFNq+Hi2UrsrWg1uQN3A6c22fYsMMzdDwUWAzcAmNkQ4AJgaPQ+fzSzwtCiFZGOadWbQWKv2wKX/h0GHpPpiLJeq8nd3V8GNjTZ9oy710d/LQdKoz+PA2a4+zZ3XwZUAqNCjFdEOpplL8Pd34DiznDZ09Dn8ExHlBPaUrm35jJgTvTnvsCKmNuqo9t2YWaTzGyemc1bu3ZtCGGISK6pqKphytxKKqpq4u+w6DG4/2zoWgoTn4YeB6Y3wByW1AlVM/sJUA9Mb9wUZzePd193nwpMBRg5cmTcfUQkf1VU1TD+znJq6yOUFBUw/fLRO69jWnEPPD4Z+pbBRTOh8z6ZCzYHJVy5m9kE4HRgvLs3JudqoF/MbqXAqsTDE5F8VFFVw+3PLWZbXYSIQ119hPKl64Mb3eGV24Jpew84ES55VIk9AQlV7mZ2KnAd8BV3/yLmpseAB8zsNqAPMBh4I+koRSRvxFbsTlBhFhcVMHpQ92CRjWf+E8qnwLBzgguUikqafZzypesZPaj7zhW/AG1I7mb2IDAG6GFm1cBPCUbHdAKetWCu5HJ3/467LzSzmcAignbNle7ekKrgRST3lC9dT219ULEXGBxzYA8mjz2Isr6dg5kd330IjvxOMFdMQfzmQqstHWk9ubv7hXE239XC/jcDNycTlIjkr9GDulNSVEBdfYTiooIgsfcqgOnnwrKXYOzP4JjJLS6yEXuAaGzpKLnvTFeoikhalQ3oxvTLR+9oqXTbAn85F9a+D2f+HwyPV0/urOkBYvSg7mmIPLfYjnOhmTNy5EifN29epsMQkXRb814w1HHrJjj/3uAEahup5w5mVuHuI+PdpspdRDJj2SswYzwU7w7/9iT0PrRddy8b0K3DJvW2COMiJhGR9nn7r8HMjnv1hsufbXdil9apcheR9HGHF38NL/0GBh4H598Hu6v6TgVV7iKSHnVbgrVOX/oNDB/P/OPvZEr5+uanHogj3nQFrU5h0EGpcheR1PvsE5hxIaycD2N/TkXpJYy/6/V2jVOPN7Yd0Hj3ZqhyF5HU+ngB/PnEYGTM+ffBsZMpX7Zhl3HqrYk3tj3eNgmocheR1PlgDjw0EXbrCpc9Bb0PA+KPU29taGNzY9s13j0+jXMXkfC5w2t3wDM3Bgn9whnByJgYsckc2tZeiXcA6Mjj3TXOXUTSp74WnvwhzL8XhowLrjot6bzLbrHj1KfMrWzTdALxxrZrvHt8Su4iEp4vNsDMS+CjV+C4H8EJP2l28q9Ymk4gfEruItImrbY/1lXCA+fBpyvgrD/BYRe0+bF3mW9GlXjSlNxFpFWtTrG75Dl4+DIoKIIJf4f+o9v9HGqvhEtDIUWkVc0OOXSHV26F6edA137w7RcSSuwSPlXuInko7BEkcXvi2z6HR6+ARY/CsLPhjD9AyR4hRC9hUHIXySMVVTXMnl/NrHkrqI94aFdt7tIT77IB7ro4mIP95F/C0d9vcXENSb+2LLM3jWAh7DXuPiy6bR/gr8BA4CPgPHevid52AzARaACucvenUxK5iOyksS++rS5YmxTCXaVoe098yXPw58vACuDih9s1B7ukT1t67ncDpzbZdj3wvLsPBp6P/o6ZDQEuAIZG7/NHMysMLVoRaVZjX7wxsRuEO6zQHV65Laa/PjcliV0TgYWjLWuovmxmA5tsHkewaDbAPcCLwHXR7TPcfRuwzMwqgVHAa+GEKyLNie2LFxYWcE5ZKWePKA1nBMrWT+FvV8D7j8PQb8K4O3bpr4fR59fC1+FJtOfey91XA7j7ajPbN7q9L1Aes191dJuIpFiiY8VbTcofL4C/fisYv37Kr2D0Fbv018NKylr4Ojxhn1CNd0Yl7uQ1ZjYJmATQv3//kMMQ6ZjaO1a81aT85v3wxA+DBTUufaLZYY5hJWVdqRqeRJP7J2bWO1q19wbWRLdXA/1i9isFVsV7AHefCkyFYOKwBOMQkSQ0m5TrtsCT/wFv3gf7Hw9nT4MuPZt9nLCSsq5UDU+iyf0xYAJwS/T7ozHbHzCz24A+wGDgjWSDFJHUiJuU138IsyYE7ZjjfhidH6blcRFhJmVdqRqOVqf8NbMHCU6e9gA+AX4K/A2YCfQHlgPnuvuG6P4/AS4D6oHJ7j6ntSA05a9I5uzUc/9sLjx2VZDMz/oTHNx0oJxkk5am/NV87iIStGGeugEq/gJ9R7Lg6Nt5ec3uao1kOc3nLiLNW7sYZl0KaxbC0VdRMfj7jJ9WoeGIOU4Th4l0VO7BghpTvwKffwzjH4Kv/pLyjzZpXdI8oMpdpCPashEenwwLH4GBx8E3p8JefQANR8wXSu4iGZDRdT+Xvw4PXw6bVsJJN8Exk7ePhmmM66bTh1LzRa167jlMyV0kzTJ2iX1DPbz8O3j5t8HcMJc9Bf1GZT4uSQn13EXSrNmFL1KgcRKudxe8CdNOgZduYf2gM5k69F4qIoMzFpeknip3kTRLV087qMRf48zIC+xfdC/1nXZj+Zg7OO35HtQuWknJy6t3qs7bEldG20nSLkruImmWrkvs335/MXfwO8YWz+efkaEsHvE7Nnsvaus/iDsHTGtxqW2TW5TcRTIgkUvs21U1L3qMb83/AZGCz/iv+ouZbqdx/yFDAFqszluKSzM25hYld5EcEFs1F7U0V/uWGnjyWlgwk+Lew1k4+nd027AP98ccEBL9X4OGSOYWJXeRJKSrBx1bNdfWR3jw9eXMnl+9c2vk/Sfg8avhi/Uw5sdw3DUMLSxmazRG2FGZJxKrZmzMLUruIgkKowfd1oNDY9XcuD6qE9Ma6dEAc66Fdx+GXl+G8bOg92EtxpjoQUkzNuYOJXeRBCXbg27PwaGxap49v5pZ81bQEHGKi4yv8U+Y8jPYuglO+E84djIUFrcYI6ATox2AkrtIgpLtQbf34NBYNX9zRCnvvP8BZ668jW4vPQ19RsC4KdBrSJti1InRjkHJXSRByfagEzo4uFNWM4eyN28Ipukd+3M46ntQGP+j3FyMOjGa/zSfu0gGtav3vWEpPH4NLJ0L/Y4MqvUeg1u+TxjPK1lL87mLZKm2nKCsWLaG2pdv58jld1JQWAJfvxXKLoOCxGcP0YnR/KfkLpLFPnj9afZ88hoOsmqe9lH0PucPHHrIlzIdluSApCYOM7OrzWyhmb1rZg+a2W5mto+ZPWtmS6LfVR5I3mucoKuiqiacB/x8LTzyXQ6ecx57sJXLa3/IFbWTeeVj1WPSNgn/pZhZX+AqYIi7bzGzmcAFwBDgeXe/xcyuB64HrgslWpEQhdV3DnXOlUhDsI7p87+A2i9YfegVnP7mKDZ5SSgnP9Vr7ziSLQOKgN3NrA7oDKwCbgDGRG+/B3gRJXfJMmEm5NCGFlZXwJM/glXzYf/j4bRb6d3zIO4sy8KDkGS9hJO7u680s/8GlgNbgGfc/Rkz6+Xuq6P7rDazfePd38wmAZMA+vfvn2gY0kElW4GGOdY76TlXPl8Dz/0c3rofuuwHZ98Fw84GMyC8k58a396xJNOW6QaMA/YHNgKzzOzitt7f3acCUyEYCploHNLxhFGBhjkJVsLj3etr4V9/hhdvCcasH/MDOP4/oNOeLd4t0QObJv7qWJJpy4wFlrn7WgAzmw0cDXxiZr2jVXtvYE0IcYpsF0YFGvYkWO2qrt1hybPw9I9h/RI4cCycekubxqwnc2DTxF8dSzLJfTkw2sw6E7RlTgLmAZuBCcAt0e+PJhukSKywKtCMjPX+eAE8c2NwIVL3A+GiWXDQV9t892QPbBrf3nEk03N/3cweAuYD9cCbBG2WLsBMM5tIcAA4N4xARRrlZAW6aRW8cDO8NR1235sVo27i8U6nMarTfpS142HUWpG20vQDOUTD2OLL6vdl22fw6v+Df/4BvAGO/A5vDZzIBfe9l/A5g6x+vZJWmn4gD2gY284aE1y3ziX84vGF2fe+1NfC/Hvgpd/A5rUw7Bw46UboNpBX51aqtSIpp+SeIzSMbYfYA12BGQ0R33nxiky+L5EGWDAL5v4KNlbBgGPhwr9C6Y7mi1orkg5K7jlCCWGH2AMdOIUFhrtvf18y0rZwh/cfhxf+C9a+D/sdCuMfhgNP2j5evVFOnjOQnKPkniNSlRBysX/b9EB30+lDqfmidvsBL63tK3f48PngZOmq+dB9MJx7DxxyRouzNqq1Iqmm5J5Dwk4IudrHb+lANyXJfnabucOSZ+CVW2HF69C1P5xxBxx2YbMLZ4ikk/4KO7Bc7uM3d6BLefuqoS5YiPrV38OaRdC1XzC/+uGXQFFJuM8lkgQl9xTJhXZHPvbxU9bPrt0M8++F16bApytg3yFw1lQY9s2dFqROhVz4W5Lso3HuKZBL7Q4ljlZsXg9v/AnemApbaqD/0XDs1TD45F1OlKZCLv0tSfppnHua5VK7IxtP7GXFAaemCl67A+bfB/Vb4OCvw7GTod+olD5t09eeS39Lkl2U3FMgH9sd6ZKuSrXZA8jHC4J++ruzwQrg0PPhmKug58GhxxAvpqavXX9Lkigl9xTQOOa2S1Wl2lL1v0sSnXgkZb4QXr0dKp+Dki4w+rsw+gro2jesl9qqeK/9yhMO1N+SJETJPUWysd2RbVJVqbZW/TcmUTzCCZE3KH34l/DZQtijJ5x4IxwxEXZP/b9d0wNQc69df0uSCCV3iSsdfe9UVaqtVf/H9C1iU/EcLrJnGGCfsJUB8PXbYPhFULx7mC+xWc0dgFSlS1iU3GUX6ep7p6pSjfu47sEVpPOmMXzBwwwv2MKqvQ5j6YgbGXT8RVBQGNbLapPmDkCq0iUsSu5kyeiMLNLevnei71+qKtXYxz2mbxHDP3kInrwHPlkAxZ3hsPNh5ET69D40lOdLhE6USqp1+OSuccS7ak/iSfb9S0ml6k6ZfUDZxnvg1UeCoYy9DwtaL18+B3brGu7zJUAtGEm1Dp/c01Wl5pL2JJ6sGoe9cQW8PQPefhA2fAgle8JhF0DZBOhzeGZiaoFaMJJKSSV3M9sbuBMYBjhwGfAB8FdgIPARcJ671yQVZQqls0rNRs0drNqaeLp1LqHADPDMtBe2boL3/g7vzIBlrwAezKF+3DUw5Ezo1CW98bRBRygQJPOSrdx/Dzzl7ueYWQnQGfgx8Ly732Jm1wPXA9cl+Twp01ilPjy/mtYuJs+qKjUEyR6sKqpq+MXjC2mIBHOq33T60DbdP3YVpcapetv1PtZuhiXPwsJHYPFTUL8Vuu0PY24I+undBrb9sdIsHwsEyU4JJ3cz2ws4HrgUwN1rgVozGweMie52D/AiWZzcG82eX01tfYSH51c3+4Fr70mwbK/Qkj1YNd7fAXen5ovaVu/TmNy21QX3KzDaluS2bgqm2F30N1jyXNBH79wDRlwCXz4PSkemZa6X9tJ0ApIpyVTug4C1wF/M7DCgAvgB0MvdVwO4+2oz2zfenc1sEjAJoH///kmEkby2fuDa04vOhQot2REbidw/9oAAtPyeb6mBD+bAokfhwxegoRa67AeHXwxDzggm8criudM1nYBkUjKfjCJgBPB9d3/dzH5P0IJpE3efCkyFYFbIJOJIWns+cG3tRbf1gNHaZfKprPyTHbGRyP0b3+vauggRgsp9p/d80ypY/HTQR1/2EkTqYa9SOOLbQUIvHdXiCkfp1tK/kaYTkExKJrlXA9Xu/nr094cIkvsnZtY7WrX3BtYkG2SqtZakEkmysQeMwsICVm7cQkVVzU73b6m6T1fln+yIjfbeP/a97ta5hI2bt3BilxV8qfIPMOeZYOIuCHroR10JQ8ZBnxFZ23Jp6d9I0wlIJiWc3N39YzNbYWYHu/sHwEnAoujXBOCW6PdHQ4k0xZr7wCWaZBuT2Oz51cyat4IZbyxndpN+fkvVfd72Zt0p27OGsq7/DFotH86FrRvBCqHfkTD2ZzD4FNj3kKxM6LFa+zfSWHbJpGQblt8HpkdHyiwF/g0oAGaa2URgOXBuks/RvMaFRlKYBJJJso0n0OojHvf+LbWD8qY3G4kEy9Etfw2q/hl8/2x1cFuX/eDg0+DAk4KvNEzWFaa2/BupSpdMSSq5u/tbQLxVQE5K5nHbbMNSuHMs9B8dVH39j4I+w6GoU2hP0doHuLWWTUv3b6myy9mqr74WVr8NVa8GiXx5eVCZA+zZBwYcAwOOCr73/FLWV+ctSUU7TyQsub3M3rpK+Mf/BElkw4fBtsJOwdWIfUcEvdq+I2CfQW1OIvE+kM19SNvassnrD/nWT2Hl/CCJV70K1fOCYYoA3Q+EAUcHo1oGHAV7D8jpZN4euTBaSnJf/i6z1+NAOHNK8PPna2DF60GSqf4XzPsL1P8xuG23rkHC7z0ceg2DXkOhx+BdFjZu7gPZ3H+t2zOEMi8+2J99DKvfgY/fjn5/B2o+Cm6zguC9LZsQTehHQZe4o2AzJp0H2bw9ZyI5I7eTe6wu+8Ih3wi+ABrqYe17QVW5an7w/bUpEKkLbi8oDpZO6zU0WMl+30NYVNmJ+vp6Il7Qpg9k3vTFY7nD5rWwvhLWfhB8rVkUfG1eu2O/bvsHk3Ed/q2gFVZ6RFZMyNWcdFfSefm3ITklt9sy7VVfC+uXwCeL4JN34ZOFQdLatHL7Ltu8mOW+L9XWiyFDD6NX/y/BPvsHLYW9+0HJHjs9ZE62XBrqg9e8cXnwVfNR0NZaH/2q/WzHvsWdg974vofAfl+G/Q6F/YZldSKPZ8rcSm595gMiDoUG13z1YK484cCUPmdO/m1ITsnftkx7FZUElXqvoew0iOeLDbBuMaxbQs3SBRR8vJgjG1bTeclMWLR558fYfR/Yqw+fFvdkZf1e9Nm3H1f2HQCf9YBlPaFzd+i8TzDyI8QTuy1yh9rP4Yv10a8NMT/HfG1aHVwk9PnH4JEd97cC2Ls/7HMA9BsVfO9+QPA/m71Ks+qioURlopLOm3ac5KScr9ybVkehVkvusHkd1CyDmir4dAV8uoKNnyxn5fJKerCRHnxKoTXzHhbtDp32jH51geI9oHg3KNotSPxFuwcHHCuA2GnLYk86RkJDhBcAAAdPSURBVOqhoS649L6hFuq2QO0XULc5+L5tU5C4G5qZ18UKdxxw9tyPdQU9WFbblZ79DmLgoIODpN61NH0HogxKesIykSyTt5V70z7qTacP5RePLwyvr2oGXXoGX/1Gbd88fW4lt1YG/8Uvtgg3jNmXy4bvGfSkv1gfzImyZUMwkmTbZ8GkV7Wbg6+tm6B+bTCTYf224HtsFU3MgcI9OOlbWBJ8LygODg7FewQJu2s/6LQX7NE9msCbfu0Dnbpur7x3er+WFTD9oOGUdU8uuWVj66GlaYwBjWKRDiGnk3vTEQlz3l2dlhEKO00tUFTEYV86CHplf4IIewRHNg73ay0mjWKRjiKnm6mNSbYwOvnU14b1pqjAMKCwwFLWV228eOWarx6cFQmtrZq+X8m+P/ESZaa1FlPY74FItsrpyr3pFYJAtF/tKb9YJhtOlrW3JRL2Va/ZONyvtZhy9spfkXbK+ROqsTIx3C1WOqfvzZaWSC713EXyTd6eUG0qHZVkIlMRpCIRZ0vvOBv+B9NUNsYkkm55ldxT/V/ulpJ0uqfvHT2oO0UFRl2Dp/T8gojkprxK7pDaqq2lJJ2R6XvTdH5BRHJP3iX3ltomyVb02TR9b/nS9dQ3BGuRNjSkri2j/rVIbsqr5N5c2ySsnndrSbql/zWE/T+KdJ1fyIaTtiLSfnmV3Jtrm4TZ886Wk3XpGNIX+77V1ke4/bnFTB57UFa8fhFpWdIXMZlZoZm9aWaPR3/fx8yeNbMl0e9pywTNXaCSrxeulA3oxpUnHNimZFtRVcOUuZVUVNW0+fEb37cCg4jDP5asY/yd5e16DBHJjKTHuZvZNQRL7e3l7qeb2W+BDe5+i5ldD3Rz9+taeowwp/xNZc89VyXTXqmoquH25xbzjyXrcDJz/YCIxJeyce5mVgp8HbgZuCa6eRwwJvrzPcCLQIvJPUzNtU3S0U7J1gNIsot8Tx57EP/6aENWXYkqIi1Ltud+O3AtsGfMtl7uvhrA3VebWdy11sxsEjAJoH///kmGkXkVVTVc+Ofy7QnwwW9nz8nHZE++6pJ9kdyTcHI3s9OBNe5eYWZj2nt/d58KTIWgLZNoHNli9vxqauuDqXtr6yPMnl+dNUkwjOScLSeSRaRtkqncjwHOMLPTgN2AvczsfuATM+sdrdp7A2vCCDTbNT06ZdvRSslZpGNJeLSMu9/g7qXuPhC4AHjB3S8GHgMmRHebADyadJQ54OwRpZQUBtMNlxQaZ48ozXRIItKBpWKc+y3ATDObCCxnp8VK81fZgG48OOmodrU+svUErIjkvlCSu7u/SDAqBndfD5wUxuNmm9aScXtaH7r6U0RSKa+uUE2lsJNxtkzZKyL5KaeX2UunsJeUy9erZkUkO6hyb6OwJ+rS2HERSaW8WmYv1XQCVESySYdZZi/VNFZcRHKFeu4iInlIyV1EJA8puYuI5CEldxGRPKTkLiKSh5TcRUTyUFaMczeztUBVAnftAawLOZxcoNfdseh1dyzted0D3L1nvBuyIrknyszmNTeAP5/pdXcset0dS1ivW20ZEZE8pOQuIpKHcj25T810ABmi192x6HV3LKG87pzuuYuISHy5XrmLiEgcSu4iInkoZ5O7mZ1qZh+YWaWZXZ/peNLBzKaZ2RozezfTsaSLmfUzs7lm9p6ZLTSzH2Q6pnQws93M7A0zezv6un+e6ZjSycwKzexNM3s807Gki5l9ZGYLzOwtM0t6gYuc7LmbWSGwGDgZqAb+BVzo7osyGliKmdnxwOfAve4+LNPxpIOZ9QZ6u/t8M9sTqADO7AD/1gbs4e6fm1kx8A/gB+5enuHQ0sLMrgFGAnu5++mZjicdzOwjYKS7h3LhVq5W7qOASndf6u61wAxgXIZjSjl3fxnYkOk40sndV7v7/OjPnwHvAX0zG1XqeeDz6K/F0a/cq8QSYGalwNeBOzMdSy7L1eTeF1gR83s1HeAD39GZ2UDgcOD1zEaSHtHWxFvAGuBZd+8Qrxu4HbgWiGQ6kDRz4BkzqzCzSck+WK4md4uzrUNUNR2VmXUBHgYmu/umTMeTDu7e4O7DgVJglJnlfSvOzE4H1rh7RaZjyYBj3H0E8DXgymgbNmG5mtyrgX4xv5cCqzIUi6RYtOf8MDDd3WdnOp50c/eNwIvAqRkOJR2OAc6I9p9nACea2f2ZDSk93H1V9Psa4BGC9nPCcjW5/wsYbGb7m1kJcAHwWIZjkhSInli8C3jP3W/LdDzpYmY9zWzv6M+7A2OB9zMbVeq5+w3uXuruAwk+1y+4+8UZDivlzGyP6IABzGwP4KtAUqPicjK5u3s98D3gaYITbDPdfWFmo0o9M3sQeA042MyqzWxipmNKg2OAbxFUcG9Fv07LdFBp0BuYa2bvEBQzz7p7hxkW2AH1Av5hZm8DbwBPuPtTyTxgTg6FFBGRluVk5S4iIi1TchcRyUNK7iIieUjJXUQkDym5i4jkISV3EZE8pOQuIpKH/j8U9OEreTBNMAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "RMSE: 10.255197731739395\n" ] } ], "source": [ "from sklearn.preprocessing import PolynomialFeatures\n", "plt.plot(alturas, pesos, '.')\n", "alturas_sorted = np.array(sorted(alturas))\n", "pol_feat = PolynomialFeatures(2).fit(X)\n", "Xpol = pol_feat.transform(X)\n", "Xpol_predict = pol_feat.transform(alturas_sorted.reshape(-1, 1))\n", "plt.plot(alturas_sorted, lm().fit(Xpol, y).predict(Xpol_predict)); plt.show()\n", "print('RMSE:', rmse(y, lm().fit(Xpol, y).predict(Xpol)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Un último ejemplo - Repaso, multicolinealidad y VIF\n", "Abordemos un último dataset, que en este caso relaciona el precio de unas viviendas en función de sus características principales:" ] }, { "cell_type": "code", "execution_count": 126, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
MetrosHabitacionesBalconBañosZonaPrecio
0803No1Periferia220000
17621Centro210000
25011Centro165000
39732Centro260000
412542Periferia380000
5834No2Centro270000
61003No1Periferia220000
\n", "
" ], "text/plain": [ " Metros Habitaciones Balcon Baños Zona Precio\n", "0 80 3 No 1 Periferia 220000\n", "1 76 2 Sí 1 Centro 210000\n", "2 50 1 Sí 1 Centro 165000\n", "3 97 3 Sí 2 Centro 260000\n", "4 125 4 Sí 2 Periferia 380000\n", "5 83 4 No 2 Centro 270000\n", "6 100 3 No 1 Periferia 220000" ] }, "execution_count": 126, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv('casas.csv')\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "En primer lugar tenemos que preprocesar las columnas Balcón y Zona:" ] }, { "cell_type": "code", "execution_count": 127, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
MetrosHabitacionesBalconBañosZonaPrecio
0803011220000
1762110210000
2501110165000
3973120260000
41254121380000
5834020270000
61003011220000
\n", "
" ], "text/plain": [ " Metros Habitaciones Balcon Baños Zona Precio\n", "0 80 3 0 1 1 220000\n", "1 76 2 1 1 0 210000\n", "2 50 1 1 1 0 165000\n", "3 97 3 1 2 0 260000\n", "4 125 4 1 2 1 380000\n", "5 83 4 0 2 0 270000\n", "6 100 3 0 1 1 220000" ] }, "execution_count": 127, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.preprocessing import LabelEncoder as le\n", "df.Balcon = le().fit_transform(df.Balcon)\n", "df.Zona = le().fit_transform(df.Zona)\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Acto seguido fijamos un modelo lineal y calculamos su score $R^2$ para medir su calidad:" ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.9814930926346503" ] }, "execution_count": 128, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X, y = df.drop(columns='Precio').values, df.Precio.values\n", "fitted_mod = lm().fit(X, y)\n", "fitted_mod.score(X, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Los coeficientes de cada variable explicativa son" ] }, { "cell_type": "code", "execution_count": 129, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ -896.13970588, 75193.01470588, 98419.11764706, 7412.68382353,\n", " 49218.75 ])" ] }, "execution_count": 129, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fitted_mod.coef_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "y como las columnas son" ] }, { "cell_type": "code", "execution_count": 130, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Index(['Metros', 'Habitaciones', 'Balcon', 'Baños', 'Zona'], dtype='object')" ] }, "execution_count": 130, "metadata": {}, "output_type": "execute_result" } ], "source": [ "col_names = df.drop(columns='Precio').columns\n", "col_names" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "y el coeficiente libre o intercepto es" ] }, { "cell_type": "code", "execution_count": 131, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "21888.786764705786" ] }, "execution_count": 131, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fitted_mod.intercept_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "concluimos que la fórmula lineal para los precios de las casas según este modelo es:" ] }, { "cell_type": "code", "execution_count": 132, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'Precio = 21888.79 + -896.14·Metros + 75193.01·Habitaciones + 98419.12·Balcon + 7412.68·Baños + 49218.75·Zona'" ] }, "execution_count": 132, "metadata": {}, "output_type": "execute_result" } ], "source": [ "intercept = np.round(fitted_mod.intercept_, 2)\n", "formula = f'Precio = {intercept}'\n", "for i in range(len(col_names)):\n", " col = col_names[i]\n", " coef = np.round(fitted_mod.coef_[i], 2)\n", " formula += f' + {coef}·{col}'\n", "formula" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nos preguntamos ahora si hay multicolinealidad entre las columnas de este dataset, causando que el modelo sea poco robusto. Calculamos los VIF para cada variable explicativa:" ] }, { "cell_type": "code", "execution_count": 133, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
columnaVIF
0Metros156.572935
1Habitaciones231.443140
2Balcon16.875664
3Baños96.266766
4Zona7.242168
\n", "
" ], "text/plain": [ " columna VIF\n", "0 Metros 156.572935\n", "1 Habitaciones 231.443140\n", "2 Balcon 16.875664\n", "3 Baños 96.266766\n", "4 Zona 7.242168" ] }, "execution_count": 133, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from statsmodels.stats.outliers_influence import variance_inflation_factor as VIF\n", "vif_data = pd.DataFrame()\n", "vif_data['columna'] = col_names\n", "vif_data['VIF'] = [VIF(df.drop(columns='Precio').values, i) for i in range(len(col_names))]\n", "vif_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vemos que todos los VIFs superan el umbral de 5 comentado en las explicaciones previas. Si removemos la columna con mayor VIF y repetimos el cálculo:" ] }, { "cell_type": "code", "execution_count": 134, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
columnaVIF
0Metros60.849074
1Balcon3.360900
2Baños32.900099
3Zona7.201042
\n", "
" ], "text/plain": [ " columna VIF\n", "0 Metros 60.849074\n", "1 Balcon 3.360900\n", "2 Baños 32.900099\n", "3 Zona 7.201042" ] }, "execution_count": 134, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = df.drop(columns = 'Habitaciones')\n", "vif_data = pd.DataFrame()\n", "vif_data['columna'] = df.drop(columns = 'Precio').columns\n", "vif_data['VIF'] = [VIF(df.drop(columns='Precio').values, i) for i in range(len(vif_data['columna']))]\n", "vif_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Aunque en general los VIFs han bajado, todavía hay algunos superiores a 5, así que removemos otra vez la columna con mayor VIF:" ] }, { "cell_type": "code", "execution_count": 135, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
columnaVIF
0Balcon2.461538
1Baños3.384615
2Zona1.615385
\n", "
" ], "text/plain": [ " columna VIF\n", "0 Balcon 2.461538\n", "1 Baños 3.384615\n", "2 Zona 1.615385" ] }, "execution_count": 135, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = df.drop(columns = 'Metros')\n", "vif_data = pd.DataFrame()\n", "vif_data['columna'] = df.drop(columns = 'Precio').columns\n", "vif_data['VIF'] = [VIF(df.drop(columns='Precio').values, i) for i in range(len(vif_data['columna']))]\n", "vif_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Y ahora todos los VIFs son inferiores a 5, y por tanto hemos evitado la multicolinealidad en el dataset. Si fijamos un nuevo modelo lineal con este dataset reducido:" ] }, { "cell_type": "code", "execution_count": 136, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Score R^2: 0.9151431858859627\n" ] } ], "source": [ "X, y = df.drop(columns='Precio').values, df.Precio.values\n", "fitted_mod = lm().fit(X, y)\n", "print('Score R^2:', fitted_mod.score(X, y))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Aunque el modelo ha disminuido un poco su calidad, ahora no existe multicolinealidad por lo que el modelo es más robusto. Calculemos cómo queda la fórmula final:" ] }, { "cell_type": "code", "execution_count": 137, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'Precio = 41428.57 + 31785.71·Balcon + 107321.43·Baños + 78214.29·Zona'" ] }, "execution_count": 137, "metadata": {}, "output_type": "execute_result" } ], "source": [ "col_names = df.drop(columns='Precio').columns\n", "intercept = np.round(fitted_mod.intercept_, 2)\n", "formula = f'Precio = {intercept}'\n", "for i in range(len(col_names)):\n", " col = col_names[i]\n", " coef = np.round(fitted_mod.coef_[i], 2)\n", " formula += f' + {coef}·{col}'\n", "formula" ] } ], "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }