Aula7-Regressao-Linear-Simples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regressão linear simples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Francisco Aparecido Rodrigues, francisco@icmc.usp.br.<br>\n",
    "Universidade de São Paulo, São Carlos, Brasil.<br>\n",
    "https://sites.icmc.usp.br/francisco <br>\n",
    "Copyright: Creative Commons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos supor que temos o conjunto de dados abaixo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAELCAYAAAA7h+qnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPBUlEQVR4nO3db4xldX3H8fdnWf8N1mjD9I/g7mBjrZa0wd40Ko21YhsaKfqgtZqRAG2dNGkVjYn/9gEx6T5oNFYfGOMEEK1XTIu0EmKt1j+xTZR2FqyCq5aAO66ijDVV6yQi4dsH926ZGRiZgdnzO3vP+5WQM/c3Z3c/OWHvZ865Z883VYUkadj2tQ4gSWrPMpAkWQaSJMtAkoRlIEkC9rcO8HCdccYZtbCw0DqGJJ1Sjhw58t2qmt+6fsqWwcLCAisrK61jSNIpJcmxB1v3MpEkyTKQJFkGkiQsA0kSloEkiY7LIMnVSe5OcuuGtbcm+UqSLyb5hyRP7DKTJJ0qxmNYWIB9+ybb8Xjvfu+uzwyuAS7YsvYJ4Jyq+jXga8CbOs4kSb03HsPSEhw7BlWT7dLS3hVCp2VQVZ8Fvrdl7eNVde/05eeBs7rMJEmngkOHYH1989r6+mR9L/TtM4M/Af5pu28mWUqykmRlbW2tw1iS1Nbq6u7Wd6s3ZZDkEHAvsO1JT1UtV9Woqkbz8w/419SSNLMOHNjd+m71ogySXAJcCCyWo9ck6QEOH4a5uc1rc3OT9b3QvAySXAC8AbioqtYfan9JGqLFRVhehoMHIZlsl5cn63shXf4gnuRa4PnAGcB3gCuY3D30GOC/p7t9vqr+/KF+r9FoVD6oTpJ2J8mRqhptXe/0qaVV9fIHWb6qywySpAdqfplIktSeZSBJsgwkSZaBJAnLQJKEZSBJwjKQJGEZSJKwDCRJWAaSJCwDSRKWgSQJy0CShGUgScIykCRhGUiSsAwkSVgGkiQsA0kSloEkCctAkoRlIEnCMpAkYRlIkrAMJElYBpIkLANJEh2XQZKrk9yd5NYNaz+b5BNJ/mu6fVKXmSTpoYzHsLAA+/ZNtuNx60R7r+szg2uAC7asvRH4ZFU9Dfjk9LUk9cJ4DEtLcOwYVE22S0uzVwidlkFVfRb43pblFwPvm379PuAlXWaSpJ/m0CFYX9+8tr4+WZ8lffjM4Oer6i6A6fbnttsxyVKSlSQra2trnQWUNFyrq7tbP1X1oQx2rKqWq2pUVaP5+fnWcSQNwIEDu1s/VfWhDL6T5BcBptu7G+eRpP93+DDMzW1em5ubrM+SPpTBDcAl068vAT7SMIskbbK4CMvLcPAgJJPt8vJkfZakqrr7w5JrgecDZwDfAa4A/hH4O+AAsAr8UVVt/ZD5AUajUa2srJy8sJI0g5IcqarR1vX9XYaoqpdv863zu8whSdqsD5eJJEmNWQaSJMtAkmQZSJKwDCRJWAaSJCwDSRKWgSQJy0CShGUgqceGMGGsLzp9HIUk7dSJCWMnBsucmDAGs/eQuD7wzEBSLw1lwlhfWAaSemkoE8b6wjKQ1EtDmTDWF5aBpF4ayoSxvrAMJPXSUCaM9YV3E0nqrcVF3/y74pmBJMkykCRZBpIkLANJEpaBJAnLQJKEZSBJwjKQJGEZSJKwDCRJ9KgMkrw2yW1Jbk1ybZLHts4kda0vk736kkPd6UUZJDkTeDUwqqpzgNOAl7VNJXXrxGSvY8eg6v7JXl2/Efclh7rVizKY2g88Lsl+YA74VuM8Uqf6MtmrLznUrV6UQVV9E3gbsArcBXy/qj6+db8kS0lWkqysra11HVM6qfoy2asvOdStXpRBkicBLwbOBp4MnJ7kFVv3q6rlqhpV1Wh+fr7rmNJJ1ZfJXn3JoW71ogyAFwJ3VtVaVf0EuB54buNMUqf6MtmrLznUrb6UwSrw7CRzSQKcDxxtnEnqVF8me/Ulh7qVqmqdAYAkbwH+GLgXuAX4s6r68Xb7j0ajWllZ6SqeJM2EJEeqarR1vTdjL6vqCuCK1jkkaYj6cplIktSQZSBJsgwkSZaBJAnLQJKEZSBJwjKQJGEZSJKwDCRJWAYS4GQvqTePo5BaOTHZ68RAlxOTvcCHs2k4PDPQ4DnZS7IMJCd7SVgGkpO9JCwDycleEpaB5GQvCe8mkoDJG79v/hoyzwwkSTsvgyQXJrE8JGkG7ebN/SPAN5P8dZJnnKxAkqTu7aYMfglYBl4K3Jrkc0lemeQJJyeaJKkrOy6Dqvp6VV1RVWcDvwvcDvwNcFeSv03yOycrpCTp5HpYnwFU1aeq6mLgl4EjwCLwL0nuTPLaJN6lJEmnkIdVBkl+O8k1wFeBc4B3Ab8H/D3wFuD9exVQknTy7fgn+CQHgUum/y0AnwGWgOur6sfT3T6Z5HPAB/Y2piTpZNrN5Zw7gG8B1wBXV9Wd2+x3G/DvjzCXJKlDuymDPwA+VlX3/bSdquprgB8mS9IpZDd3E330oYrgkUjyxCTXJflKkqNJnnOy/ixJ0mZ9uuvnnUzOPP4wyaOBuYf6BZKkvdGLMpj+w7XnAZcCVNU9wD0tM0nSkPTlWUNPBdaA9ya5JcmVSU7fulOSpSQrSVbW1ta6TylJM6ovZbAfeBbw7qo6F/gR8MatO1XVclWNqmo0Pz/fdUZJmll9KYPjwPGqumn6+jom5SBJ6kAvyqCqvg18I8nTp0vnA19uGEmSBqUXHyBPvQoYT+8kugO4rHEeSRqM3pRBVX0BGLXOIUlD1IvLRJKktiwDSZJlIEmyDCRJWAaSJCwDSRKWgSQJy0CShGUgScIyUA+Mx7CwAPv2TbbjcetE0vD05nEUGqbxGJaWYH198vrYsclrgMXFdrmkofHMQE0dOnR/EZywvj5Zl9Qdy0BNra7ubl3SyWEZqKkDB3a3LunksAzU1OHDMDe3eW1ubrIuqTuWgZpaXITlZTh4EJLJdnnZD4+lrnk3kZpbXPTNX2rNMwNJkmUgSbIMJElYBpIkLANJEpaBJAnLQJKEZSBJwjKQJGEZSJLoWRkkOS3JLUlubJ1FkoakV2UAXA4cbR1CkoamN2WQ5CzgRcCVrbNI0tD0pgyAdwCvB+7bbockS0lWkqysra11l0ySZlwvyiDJhcDdVXXkp+1XVctVNaqq0fz8fEfpJGn29aIMgPOAi5J8HfgQ8IIkH2gbSZKGoxdlUFVvqqqzqmoBeBnwqap6ReNYkjQYvSgDSVJbvRt7WVWfAT7TOIYkDYpnBpIky0CSZBlIkrAMJElYBpIkLANJEpaBJAnLQJKEZSBJwjJoYjyGhQXYt2+yHY9bJ5I0dL17HMWsG49haQnW1yevjx2bvAZYXGyXS9KweWbQsUOH7i+CE9bXJ+uS1Ipl0LHV1d2tS1IXLIOOHTiwu3VJ6oJl0LHDh2FubvPa3NxkXZJasQw6trgIy8tw8CAkk+3ysh8eS2rLu4kaWFz0zV9Sv3hmIEmyDCRJloEkCctAkoRlIEnCMpAkYRlIkrAMJElYBpIkLANJEj0pgyRPSfLpJEeT3Jbk8taZhsCJa5JO6Muzie4FXldVNyf5GeBIkk9U1ZdbB5tVTlyTtFEvzgyq6q6qunn69Q+Bo8CZbVPNNieuSdqoF2WwUZIF4Fzgpgf53lKSlSQra2trXUebKU5ck7RRr8ogyeOBDwOvqaofbP1+VS1X1aiqRvPz890HnCFOXJO0UW/KIMmjmBTBuKqub51n1jlxTdJGvSiDJAGuAo5W1dtb5xkCJ65J2ihV1ToDSX4L+FfgS8B90+U3V9VHt/s1o9GoVlZWuognSTMjyZGqGm1d78WtpVX1b0Ba55CkoerFZSJJUluWgSTJMpAkWQaSJCwDSRKWgSQJy0CShGUgScIykCQxsDJwspckPbhePI6iC072kqTtDebMwMlekrS9wZSBk70kaXuDKQMne0nS9gZTBk72kqTtDaYMnOwlSdsbzN1EMHnj981fkh5oMGcGkqTtWQaSJMtAkmQZSJKwDCRJQKqqdYaHJckacOxh/vIzgO/uYZxTncfjfh6LzTwem83C8ThYVfNbF0/ZMngkkqxU1ah1jr7weNzPY7GZx2OzWT4eXiaSJFkGkqThlsFy6wA94/G4n8diM4/HZjN7PAb5mYEkabOhnhlIkjawDCRJwyuDJBck+WqS25O8sXWeVpI8JcmnkxxNcluSy1tn6oMkpyW5JcmNrbO0luSJSa5L8pXp/yfPaZ2plSSvnf49uTXJtUke2zrTXhtUGSQ5DXgX8PvAM4GXJ3lm21TN3Au8rqqeATwb+IsBH4uNLgeOtg7RE+8EPlZVvwL8OgM9LknOBF4NjKrqHOA04GVtU+29QZUB8JvA7VV1R1XdA3wIeHHjTE1U1V1VdfP06x8y+Yt+ZttUbSU5C3gRcGXrLK0leQLwPOAqgKq6p6r+p22qpvYDj0uyH5gDvtU4z54bWhmcCXxjw+vjDPwNECDJAnAucFPbJM29A3g9cF/rID3wVGANeO/0stmVSU5vHaqFqvom8DZgFbgL+H5Vfbxtqr03tDLIg6wN+t7aJI8HPgy8pqp+0DpPK0kuBO6uqiOts/TEfuBZwLur6lzgR8AgP2NL8iQmVxDOBp4MnJ7kFW1T7b2hlcFx4CkbXp/FDJ7u7VSSRzEpgnFVXd86T2PnARcl+TqTy4cvSPKBtpGaOg4cr6oTZ4vXMSmHIXohcGdVrVXVT4Drgec2zrTnhlYG/wE8LcnZSR7N5EOgGxpnaiJJmFwPPlpVb2+dp7WqelNVnVVVC0z+v/hUVc3cT387VVXfBr6R5OnTpfOBLzeM1NIq8Owkc9O/N+czgx+m728doEtVdW+SvwT+mckdAVdX1W2NY7VyHnAx8KUkX5iuvbmqPtowk/rlVcB4+oPTHcBljfM0UVU3JbkOuJnJXXi3MIOPpfBxFJKkwV0mkiQ9CMtAkmQZSJIsA0kSloEkCctAkoRlIEnCMpAkYRlIj9h0CMzxJO/fsn5Dkq8lmWuVTdopy0B6hKbP+f9T4OIkLwFIchmT2QiXVtV6y3zSTvg4CmmPJHkP8BLgAuDTwHuq6g1tU0k7YxlIe2Q6G+KLTJ55fzvwG1X147appJ3xMpG0R6rqf4EbgccAV1kEOpV4ZiDtkSQj4HPAl4CDwK9O5wJIvWcZSHsgyWOZPO/+DuClwH8yGRx0UdNg0g55mUjaG38F/ALwyundQ5cAL0pyadNU0g55ZiA9QknOAz4LXFxVH9yw/lbglcA5VXW8VT5pJywDSZKXiSRJloEkCctAkoRlIEnCMpAkYRlIkrAMJElYBpIk4P8AjNEyvizrObsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline \n",
    "import numpy as np\n",
    "\n",
    "# define os dados\n",
    "x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) \n",
    "y = np.array([1, 3, 2, 5, 7, 8, 8, 9, 10, 12]) \n",
    "plt.plot(x, y, 'bo')\n",
    "plt.xlabel(\"x\", fontsize = 15)\n",
    "plt.ylabel(\"y\", fontsize = 15)\n",
    "plt.show(True)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos que é uma relação entre as variáveis $X$ e $Y$. Queremos determinar um modelo que melhor descreva essa relação. Podemos supor um modelo de regressão linear simples, como:\n",
    "$$\n",
    "y_i \\approx \\beta_0 + \\beta_1x_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como há outros fatores, além de $x_i$ que afetam os valores de $y_i$, podemos escrever:\n",
    "$$\n",
    "y_i = \\beta_0 + \\beta_1 x_i + \\epsilon_i,\n",
    "$$\n",
    "onde $\\epsilon_i$ é uma variável aleatória que indica o erro na aproximação. O objetivo dos métodos de regressão é encontrar o melhor valor de $\\beta_0$ e $\\beta_1$ que minimizem o erro no ajuste. Ou seja, queremos encontrar a linha no plano $x-y$ que melhor se ajuste aos dados observados."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimando os coeficiente através do método dos momentos ou dos mínimos quadrados, obtemos:\n",
    "    $$\n",
    "    \\begin{cases}\n",
    "    \\hat{\\beta}_1=\\frac{\\sum_{i=1}^n (x_i- \\bar{x})(y_i - \\bar{y})}{\\sum_{i=1}^n (x_i - \\bar{x})^2} = \\frac{S{xy}}{S{xx}}\\\\\n",
    "    \\hat{\\beta}_0= \\bar{y}-\\hat{\\beta_1}\\bar{x}\n",
    "    \\end{cases}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos implementar uma função para realizar a estimação."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statistics import variance \n",
    "import math \n",
    "\n",
    "def estimate_coef(x, y): \n",
    "    # número de observações/pontos\n",
    "    n = np.size(x) \n",
    "  \n",
    "    # médias de x e y\n",
    "    m_x, m_y = np.mean(x), np.mean(y) \n",
    "  \n",
    "    # calculating cross-deviation and deviation about x \n",
    "    SS_xy = np.sum(y*x) - n*m_y*m_x \n",
    "    SS_xx = np.sum(x*x) - n*m_x*m_x \n",
    "  \n",
    "    # calcula os coeficientes de regressão\n",
    "    b_1 = SS_xy / SS_xx \n",
    "    b_0 = m_y - b_1*m_x \n",
    "  \n",
    "    return(b_0, b_1) \n",
    "\n",
    "# função para mostrar os dados e o ajuste linear\n",
    "def plot_regression_line(x, y, b): \n",
    "    # mostra os dados\n",
    "    plt.scatter(x, y, color = \"b\", marker = \"o\", s = 50) \n",
    "  \n",
    "    # prediz os valores\n",
    "    y_pred = b[0] + b[1]*x \n",
    "  \n",
    "    # mostra a reta de regressão\n",
    "    plt.plot(x, y_pred, color = \"r\") \n",
    "  \n",
    "    plt.xlabel('x', fontsize = 15) \n",
    "    plt.ylabel('y', fontsize = 15) \n",
    "    plt.show(True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assim, aplicando ao conjunto de dados:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated coefficients:\n",
      "b_0 = 1.2363636363636363  \n",
      "b_1 = 1.1696969696969697\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAELCAYAAAA7h+qnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAf+0lEQVR4nO3dd3SVVbrH8e8DWIFRRuN4xYKOjICoiBlFwK6jg6jYANu1MxaUsRcso14cB1FhRCkLFAsicxFEEZAiqCjFUASkCGIBRY0FCVFpee4fOxjMTTDAydlvzvl91nKFvCee86yz4Pyy9373fszdERGR7FYtdgEiIhKfwkBERBQGIiKiMBARERQGIiIC1IhdwJbaddddvV69erHLEBGpUqZPn/6Nu+eUvl5lw6BevXrk5eXFLkNEpEoxs0/Luq5pIhERURiIiIjCQEREUBiIiAgKAxERIc13E5nZU0Br4Gt3b1x87WHgNGAN8BFwqbuvSGddIiJJV1AAgwfDokVQvz60awe1a6fu+dM9MhgAnFLq2ligsbsfDHwI3JHmmkREEm3SJKhbF/7+d+jaNXytWzdcT5W0hoG7vwV8V+raGHdfV/ztFGDPdNYkIpJkBQXQqlX4WlgYrhUWllxftSo1r5O0NYPLgFGxixARSYrBg6GoqOzHiorC46mQmDAws87AOmDgJn6mg5nlmVlefn5++ooTEYlk0aKSEUFphYWweHFqXicRYWBmFxMWli/wTbRec/e+7p7r7rk5Of/vaA0RkYxTvz7UrFn2YzVrwv77p+Z1ooeBmZ0C3Aac7u4/xq5HRCRJ2rWDar98Uju78M0vj1WrFh5PhbSGgZkNAiYDB5jZMjO7HOgJ1AbGmtksM+udzppERJKsdm0YORIa1fyUEdXPYArN2GXHn365XqtWal4nrfsM3P28Mi73T2cNIiJVytq1tJzcnbn+D9ZVhzFH30/X9tvQ9vzUBQFU4SOsRUQy3pQp8Le/wezZ2Omns83jj3Pq3ntXyktFXzMQEZFSVqyAq6+G5s3hu+9g2DAYPhwqKQhAYSAikhzuMGgQNGgAffuGrcbz5kGbNpX+0pomEhFJgsWL4ZprYOxY+POfYdQoOPTQtL28RgYiIjGtXg0PPACNG8PUqdCzJ0yenNYgAI0MRETiefNNuOoqWLAA2raFxx6DPfaIUopGBiIi6fbNN3DppXDssWFkMGpUOGQoUhCAwkBEJH3c4emnwwLx88/DHXfA3LlwSumT/dNP00QiIukwb164XfStt6BlS+jdGw48MHZVv9DIQESkMv30E3TuDE2awJw50K9fWCtIUBCARgYiIpXn9dfD7aJLlsBFF0G3brDbbrGrKpNGBiIiqfbll3DeeWEtoEYNeOMNePbZxAYBKAxERFKnqAh69QoLxMOGwX33wezZcNxxsSv7TZomEhFJhfffD4fKTZ0KJ5wQQqF+/dhVVZhGBiIiW2PVKrj5ZjjsMPj443DL6NixVSoIQCMDEZEtN3w4XHcdLF0KHTrAQw9BnTqxq9oiGhmIiGyupUvDSaJt2sBOO8GkSdCnT5UNAlAYiIhU3Lp14fyghg1hzJgwEpgxA1q0iF3ZVtM0kYhIRUybFhaIZ82CVq3C6aL77hu7qpTRyEBEZFN++AE6doRmzeDrr2HIEBgxIqOCABQGIiJlcw8niTZoEG4T7dgR5s+Hs88Gs9jVpZymiURESluyBK69FkaPhqZN4dVXITc3dlWVSiMDEZEN1qyBBx8Mh8hNmgTdu4dNZBkeBKCRgYhIMGlSWCCeNw/OOgt69IA994xdVdpoZCAi2e3bb+GKK+Coo8Ju4ldfhZdeyqogAIWBiGQr93CSaIMGMGAA3HJLGBW0bh27sig0TSQi2WfhwtB1bMKEcMtonz5w8MFl/mhBQbipaNGicNxQu3ZQu3aa602DtIaBmT0FtAa+dvfGxdd+DwwG6gGfAG3d/ft01iUiWeLnn8MC8b/+BTvuGFpPXnklVCt7kmTSpLC/rKgICguhZk248UYYOTJ0rswk6Z4mGgCU7vx8OzDe3esD44u/FxFJrXHj4KCD4IEH4JxzYMGCsGBcThAUFIQgKCgIQQDh64brq1alsfY0SGsYuPtbwHelLp8BPFP852eANumsSUQy3FdfwYUXwkknhe/HjIGBA+EPf9jk/zZ4cBgRlKWoKDyeSZKwgPwHd18OUPy13L5wZtbBzPLMLC8/Pz9tBYpIFVRUBH37hgXi//wH7r47dB3bEAq/YdGikhFBaYWFsHhxCmtNgCSEQYW5e193z3X33JycnNjliEhSzZkTbhX929/gkENCF7L774cddqjwU9SvH9YIylKzJuy/f4pqTYgkhMFXZvZfAMVfv45cj4hUVYWFcNtt4QiJhQvh6afDHUMNG272U7VrV+5yAtWqhcczSRLC4BXg4uI/XwwMj1iLiFRVr70WjpHo2hX++7/DAvEll2zxoXK1a4e7hmrXLhkh1KxZcr1WrdSVngTpvrV0EHAssKuZLQPuBR4C/mNmlwOfAeemsyYRqeI+/xw6dQq7hhs2hDffhKOPTslTt2wJX3wRFosXLw5TQ+3aZV4QQJrDwN3PK+ehE9JZh4hkgPXr4Ykn4K67YO1a6NIlNKbfdtuUvkytWnD55Sl9ykTSDmQRqXqmTw+Lw9Onw1/+Ak8+CX/8Y+yqqrQkrBmIiFTMypVhSujww2HZMhg0KPQcUBBsNY0MRCT53GHoULj+eli+HK66KhwrsfPOsSvLGBoZiEiyffIJnHZaOEIiJwcmTw7TQgqClFIYiEgyrV0bbhM98MCwV6BbN8jLgyOOiF1ZRtI0kYgkz+TJYYF4zhw4/XR4/HHYe+/YVWU0jQxEJDm+/z6EQPPm4c/DhsHw4QqCNFAYiEh87vDCC+FQuX794IYbQtexNjrEOF00TSQicS1eHLqOjRsHf/5zuFX00EOB7OkylgQKAxGJY/XqsEDcpQtstx307BluGa1eHciuLmNJoDAQkfSbODF88C9cCG3bwmOPwR57/PLwxl3GNtjQW6BVq3BeUCaeDxST1gxEJH3y8+Hii+G442DNGhg1KswDbRQEkH1dxpJAYSAila+oCPr3DwvEL7wAd9wBc+fCKaVbogfZ1mUsCTRNJCKVa968MCX09tthsr9377CRbBM2dBkrKxAysctYEmhkICKV46efoHNnaNIkjAL69Qu9Bn4jCCD7uowlgcJARFLv9dehceNwmFz79qHr2OWXl/8JX0q2dRlLAk0TiUjqLF8eNowNHgx/+hOMHw/HH79FT5VNXcaSQGEgIltv/Xro0ycsDK9eDffdFxrTb7fdVj1ttnQZSwKFgYhsnVmzwnlC06bBCSdAr15hBViqFK0ZiMiWWbUKbroJcnNDz4Hnn4exYxUEVZRGBiKy+YYPh+uug6VLoUMHeOghqFMndlWyFTQyEJGKW7o0nCTapg3stFM4QKhPHwVBBlAYiMhvW7cunB/UsCGMGRNGAjNmQIsWsSuTFNE0kYhs2rRpYYF41qxwSlzPnrDvvrGrkhTTyEBEyvbDD9CxIzRrBl9/Df/7vzBihIIgQykMROTX3OE//wlTQr16hUCYPx/OOQfMYlcnlUTTRCIJEr2z15IlcO21MHo03+zdlBcueIUdD86lnYEajGU2c/fYNQBgZjcAVwAOzAEudfefy/v53Nxcz8vLS1d5IpWurM5e1aqlqbPXmjXQrRs88ADrrAadi/6H3tWvZeWPNdJbh1Q6M5vu7rmlrydimsjM6gLXA7nu3hioDrSPW5VI+mzc2WvDsc2FhSXXV62qxBefNCn0HO7cmbV/aUXjavPpuroTK3+skd46JKpEhEGxGsAOZlYD2BH4InI9ImkTpbPXt9/CFVfAUUeFT/lXX+WZ015iGXumtw5JhESEgbt/DnQDPgOWAz+4+5jSP2dmHcwsz8zy8vPz012mSKVJa2cvd3j22dB1bMAAuOWW0ICmdWt1GMtiiQgDM6sDnAHsC+wB1DSzC0v/nLv3dfdcd8/NyclJd5kilWZDZ6+ypLSz18KF4TC5iy8OTzpjBnTt+suLp60OSZxEhAFwIvCxu+e7+1pgKNA8ck0iaVPpnb1+/hnuvRcOPhhmzgytJ995J3yfzjoksZISBp8BzcxsRzMz4ARgfuSaRNKmUjt7jR8fPvTvvz/sFViwIOwoLuNTXx3Gslci9hm4+1QzGwLMANYBM4G+casSSa+Ud/b66qtwxPTAgeHJxoyBk05Kfx1SJSRmn8Hm0j4DkXIUFYXm87fdFlZ9b78d7rwTtt8+dmWSAOXtM0jEyEBEUmTOHLjqKnj3XTjmmLA20KBB7KqkCkjKmoGIbI3CwjASaNo03DE0YABMmKAgkArTyECkqnvttXCe0KefwmWXhVtFd9kldlVSxWhkIFJVff55uDuodWvYcUd4803o319BIFtEYSBS1axfD//+dzhi+rXXoEuX0Hjm6KNjVyZVmKaJRKqS6dPDHoHp0+Hkk+GJJ+CPf4xdlWQAjQxEqoKVK6FTJzj88DA99OKLMGqUgkBSRiMDkSRzh6FD4frrYflyuPrqMC20886xK5MMozAQIQEdxsryySeh5eRrr8Ehh4RQOOKIyEVJplIYSNYrq8PYjTdG7Oy1di107w7/+EfoOfzII2FkUEP/XKXy6G+XZLWNO4xtsOE8/1atwhk9aT2TZ/LksEA8Zw6ccUa4a2jvvdNYgGQrLSBLVovSYaws338fQqB58/DnYcPg5ZcVBJI2CgPJatE7e7nDCy+EYyP69w/zU/PnQ5s2lfzCIr+maSLJahs6e5UVCJXe2WvRIrjmGhg3Ltwy+vrr0KRJJb6gSPk0MpCsFqWz1+rVodHMQQfBtGlh49i77yoIJCqFgWS1tHf2mjgx3CZ6771hKmjBgjA6qF49xS8ksnk0TSRZLy2dvfLz4eab4dlnYd99w+7hU05J4QuIbB2FgQjhg//yyyvhiYuK4Omn4dZbw/2rd94JnTuHU0ZFEqTC00Rm1trMNK0kUlHz5sGxx8IVV8CBB4aTRbt0URBIIm3Oh/tw4HMz+5eZNaysgkSqvJ9+Cr/9N2kCH3wQbhmdOBEaNYpdmUi5NicM/gj0BdoCc81sspldaWa/q5zSRKqg11+Hxo3hwQfh/PPDAvFll5V/y5JIQlT4b6i7f+Lu97r7vsBJwGLgMWC5mT1nZsdVVpEiibd8ObRvHxaFt9km9B8eMABycmJXJlIhW/Triru/4e4XAX8CpgMXAOPM7GMzu8HMtDAt2WH9+rBPoEGDcHzE/ffD+++HtQKRKmSLPrTN7BjgUuBsYC3wBPAycDJwH/Bn4PwU1SiSTLNmhfOEpk2DE0+EJ58MW5pFqqAKh4GZ7QNcXPxfPWAi0AEY6u6ri39svJlNBp5PbZkiCbJqVdg01qNHaD4/cCCcd144blqkitqckcES4AtgAPCUu39czs99AEzbyrpEkmn4cLjuOli6NIwK/vlPqFMndlUiW21zwuA0YLS7l3Pgb+DuHwJaTJbMsnRpCIHhw8OZQoMHw5FHxq5KJGU2526ikb8VBFvDzHY2syFmtsDM5puZ/qVJfOvWwaOPQsOGMHYsdO0K06crCCTjJOmunx6Ekcc5ZrYtoG2aEte0aWEqaNYsOPVU6NkT6tWLXZVIpUjETpjijWtHA/0B3H2Nu6+IW5VkrR9+gGuvhWbNwgFzL70Er76qIJCMlogwAPYD8oGnzWymmfUzs5qlf8jMOphZnpnl5efnp79KyWzuYS2gQQPo3Ts0oZ8/H846S3cKScZLShjUAJoCvdz9UKAQuL30D7l7X3fPdffcHO3slFRasgT++tewi7hu3TBF1L17aGwgkgWSEgbLgGXuPrX4+yGEcBCpXGvWhHOEDjwwdBvr0QOmToXDDotdmUhaJWIB2d2/NLOlZnaAuy8ETgDmxa5LMtzbb8NVV4Wjps85J4wE6taNXZVIFIkIg2LXAQOL7yRaQjjuQiT1vv02NJt56inYZx8YMSLcLSSSxRITBu4+C8iNXYdkMPfQdvLmm2HFCrjtNrj77pLmxyJZLDFhIFKpFiyAq68OTWaaNw93Cx10UOyqRBIjKQvIIpXj55/hnnvgkEPC5rG+fcNagYJA5Fc0MpDMNW5cGA0sXgwXXgiPPAK77Ra7KpFE0shAMs9XX8EFF8BJJ4XNYmPHwnPPKQhENkFhIJmjqAj69Ak7iIcMCdNDs2eHxjMiskmaJpLMMGdOOFRu8uTQcrJ3bzjggNhViVQZGhlI1VZYGPYMHHooLFoEzzwDb7yhIBDZTBoZSNU1YgR07AiffgpXXAEPPRTaUIrIZlMYSHQFBeGw0EWLQj/5du1+43y4ZcugUycYOjScKfT229CyZdrqFclECgOJatIkaNUqrP0WFobNwDfeCCNHlvH5vn59aDBz112hA9mDD8JNN8G220apXSSTKAwkmoKCEAQFBSXXCgvD11at4IsvoFat4gfy8sIC8YwZcMop8MQTsN9+aa9ZJFNpAVmiGTw4jAjKUlQUHmflytBk5ogjQjoMHhyGDQoCkZRSGEg0ixaVjARKKyx0th8xJDSi79kz7CResADatlXXMZFKoDCQaOrXL/vA0H34hFHVW3PBy+eGXcNTpoRA2Gmn9BcpkiUUBhJNu3ZQbaO/gTVYyy10ZR6NaLn+TVb/81F47z04/PB4RYpkCYWBRFO7dpj+r10bjt/+XaZzGF25jfE1TmbekPlsd/sNUEP3OIikg8JAomrZ6Du+PedvjP+5BXv/bgVjOw7nuO+HcfjZe8UuTSSr6NcuicMdBg6EG29km+++g5tuYud//IOTfrmXVETSSWEg6ffhh3DNNTB+fLhldOzY0HxGRKLRNJGkz+rVcP/9cPDBYRPZk0/CO+8oCEQSQCMDSY8JE8JegYULoX17eOwx2H332FWJSDGNDKRy5efDxRfD8cfD2rUwejQMGqQgEEkYhYFUjqIi6N8/dB0bNAg6d4a5c+Hkk2NXJiJl0DSRpN4HH8BVV4UjSY86KnQda9QodlUisgkaGUjq/Pgj3HEHNGkC8+aFkcHEiQoCkSpAIwNJjdGjw+2iH38c1ggefhhycmJXJSIVpJGBbJ3ly8MhQ3/9a2gyM2ECDBigIBCpYhIVBmZW3cxmmtmI2LXIb1i/PjSYadAAhg8P+wfefx+OPTZ2ZSKyBZI2TdQJmA/8LnYhsgkzZ4auY++9ByeeGDaP1a8fuyoR2QqJGRmY2Z7AqUC/2LVIOQoK4IYbIDcXPv00nC00ZoyCQCQDJCYMgO7ArUA5jRDBzDqYWZ6Z5eXn56evMoGXXw53BXXvDldeGbqOnX++uo6JZIhEhIGZtQa+dvfpm/o5d+/r7rnunpujBcr0+OwzOOMMOPNMqFMH3n037BuoUyd2ZSKSQokIA6AFcLqZfQK8CBxvZs/HLSnLrVsHjzwSRgPjxkHXrjB9Ohx5ZOzKRKQSJCIM3P0Od9/T3esB7YE33P3CyGVlrylTwrrAzTfDcceFDWS33ALbbBO7MhGpJIkIA0mIFSvCyaLNm8M338BLL8Err8A++8SuTEQqWdJuLcXdJwITI5eRXdzhxRfDnUL5+XD99fDAA6E5sYhkhcSFgaTZRx+FYyTGjIHDDoPXXgtfRSSraJooW61ZA126QOPGMHky/PvfMHWqgkAkS2lkkI3eeiscMT1/Ppx9NvToAXXrxq5KRCLSyCCbfPMNXHYZHHNMOG56xAgYMkRBICIKg6zgHk4SbdAAnnsObr01NKA59dTYlYlIQmiaKNMtWBCmhN58M2wY69MHDjoodlUikjAaGWSqn36Cu++Ggw8OR0v36RPaUCoIRKQMGhlkorFjw+axjz6CCy+Ebt3gD3+IXZWIJJjCIIKCAhg8GBYtCqc/t2uXov1dX34JN94IgwaFJx47NvQbEBH5DQqDNJs0CVq1gqIiKCyEmjXD5/fIkdCy5RY+aVER9O0Lt98epofuuSc0pt9++5TWLiKZS2GQRgUFIQgKCkquFRaGr61awRdfQK1am/mks2eHrmNTpoSWk717wwEHpKpkEckSWkBOo8GDwy/xZSkqCo9XWGFhOEm0aVNYvBieeQbeeENBICJbRCODNFq0qGQkUFphYfhMr5BXX4WOHUPjmcsvh3/9C3bZJWV1ikj20cggjerXD2sEZalZE/bf/zeeYNkyOOssOP30MJ/01lvQr5+CQES2msIgjdq1g2rlvOPVqoXHy7RuXeg93LAhjBoFDz4IM2fCUUdVWq0ikl0UBmlUu3a4a6h27ZIRQs2aJdfLXDzOy4Mjjgi9Blq2DMdI3HEHbLttWmsXkcymNYM0a9ky3DU0eHBYI9h//zAi+H9BsHIl3HUXPPEE7LZb+B/OPRfMotQtIplNYRBBrVph3bdM7qHdZKdOsHx5aDzTpQvstFNaaxSR7KJpoiT55BNo3TqMAHbbLewd6NlTQSAilU5hkARr14bbQxs1CqeLPvoovPceHH547MpEJEtomii2d98NO4jnzoU2bUL7yb32il2ViGQZjQxi+f576NABWrSAH36A4cNh2DAFgYhEoTBIN3d4/vlwbMRTT8FNN8G8eWEjmYhIJJomSqcPPwx3B40fH/YOjBkDTZrErkpERCODtFi9Gu67L3QZy8uDJ5+Ed95REIhIYmhkUNkmTAg9iD/8ENq3h8ceg913j12ViMivKAwqS35+WA947jnYbz8YPRpOPjl2Vb9SaR3XRKTKSUQYmNlewLPA7kAR0Nfde8StagsVFYWF4VtvhVWroHPn8N8OO8Su7FcqpeOaiFRZiQgDYB1wk7vPMLPawHQzG+vu82IXtlk++CBMCU2aBEcfDb16hY1kCVMpHddEpEpLxAKyuy939xnFfy4A5gN141a1GX78MZwk2qQJzJ8fRgYTJyYyCCDFHddEJCMkZWTwCzOrBxwKTC3jsQ5AB4C99947rXWVa9QouPZa+PhjuOQSePhh2HXX2FVtUso6rolIxkjEyGADM6sFvAT83d1Xln7c3fu6e6675+bk5KS/wI198QW0bRvmVbbbLowEnn468UEAKei4JiIZJzFhYGbbEIJgoLsPjV1PudavDyeJNmwIr7wCDzwAs2bBMcfErqzCtrjjmohkrESEgZkZ0B+Y7+6Pxq6nXDNnwpFHwnXXhR3Ec+eGBjTbbRe7ss2yRR3XRCSjJWXNoAVwETDHzGYVX7vT3UdGrKlEQQHcc084UTQnB154IWwgq8JdxyrccU1EskIiwsDdJwHJ/GR9+eUwEvj883DU9D//CTvvHLuqlNhkxzURySqJmCZKpM8+gzPOgDPPhN//PvQd6NUrY4JARGRjCoPS1q6Fbt3CAvG4ceFW0bw8aNYsdmUiIpUmEdNEiTFlSpgKmj0bTjsNHn8c9tkndlUiIpVOIwOAFSvg6quheXP49lsYOjR0HlMQiEiWyO4wcIdBg6BBA+jbFzp1CsdJnHlmlb5TSERkc2XvNNFHH4WuY2PGQG5uuMG+adPYVYmIRJF9I4M1a6BLF2jcGCZPDusCU6YoCEQkq2XfyOCtt8Ku4XPPhe7dYY89YlckIhJd1oVBwREnMu6uPKasOYz6I9XdS0QEsiwMSrp7HabuXiIiG8maNYONu3ttOMu/sLDk+qpVcesTEYkpa8JA3b1ERMqXNWGg7l4iIuXLmjBQdy8RkfJlTRiou5eISPmyJgzU3UtEpHxZdWupunuJiJQtq8IA1N1LRKQsWTNNJCIi5VMYiIiIwkBERBQGIiICmLvHrmGLmFk+8OlWPMWuwDcpKqeq03tRQu/Fr+n9KJEp78U+7p5T+mKVDYOtZWZ57p4bu44k0HtRQu/Fr+n9KJHp74WmiURERGEgIiLZHQZ9YxeQIHovSui9+DW9HyUy+r3I2jUDEREpkc0jAxERKaYwEBGR7AsDMzvFzBaa2WIzuz12PTGZ2V5mNsHM5pvZB2bWKXZNsZlZdTObaWYjYtcSm5ntbGZDzGxB8d+RI2PXFIuZ3VD8b2SumQ0ys+1j15RqWRUGZlYdeAL4K9AIOM/MGsWtKqp1wE3u3hBoBlyb5e8HQCdgfuwiEqIHMNrdGwCHkKXvi5nVBa4Hct29MVAdaB+3qtTLqjAADgcWu/sSd18DvAicEbmmaNx9ubvPKP5zAeEfe924VcVjZnsCpwL9YtcSm5n9Djga6A/g7mvcfUXcqqKqAexgZjWAHYEvIteTctkWBnWBpRt9v4ws/vDbmJnVAw4FpsatJKruwK1AUexCEmA/IB94unjarJ+ZldNFPLO5++dAN+AzYDnwg7uPiVtV6mVbGFgZ17L+3lozqwW8BPzd3VfGricGM2sNfO3u02PXkhA1gKZAL3c/FCgEsnKNzczqEGYQ9gX2AGqa2YVxq0q9bAuDZcBeG32/Jxk43NscZrYNIQgGuvvQ2PVE1AI43cw+IUwfHm9mz8ctKaplwDJ33zBSHEIIh2x0IvCxu+e7+1pgKNA8ck0pl21h8B5Q38z2NbNtCYtAr0SuKRozM8Kc8Hx3fzR2PTG5+x3uvqe71yP8vXjD3TPut7+KcvcvgaVmdkDxpROAeRFLiukzoJmZ7Vj8b+YEMnAxPat6ILv7OjPrCLxOuCPgKXf/IHJZMbUALgLmmNms4mt3uvvIiDVJclwHDCz+xWkJcGnkeqJw96lmNgSYQbgDbyYZeDSFjqMQEZGsmyYSEZEyKAxERERhICIiCgMREUFhICIiKAxERASFgYiIoDAQEREUBiJbrbgJzDIze7bU9VfM7EMz2zFWbSIVpTAQ2UrF5/xfDlxkZm0AzOxSQm+ES9z9x5j1iVSEjqMQSREz6wO0AU4BJgB93P22uFWJVIzCQCRFivtCzCaceb8YOMzdV8etSqRiNE0kkiLuvgoYAWwH9FcQSFWikYFIiphZLjAZmAPsAxxY3BdAJPEUBiIpYGbbE867XwK0Bd4nNA06PWphIhWkaSKR1PgfYHfgyuK7hy4GTjWzS6JWJVJBGhmIbCUzawG8BVzk7i9sdP1h4Eqgsbsvi1WfSEUoDERERNNEIiKiMBARERQGIiKCwkBERFAYiIgICgMREUFhICIiKAxERAT4PxhdJibXIZh7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# estima os coeficientes\n",
    "b = estimate_coef(x, y) \n",
    "print(\"Estimated coefficients:\\nb_0 = {}  \\nb_1 = {}\".format(b[0], b[1])) \n",
    "  \n",
    "# mostra o ajuste linear\n",
    "plot_regression_line(x, y, b) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para quantificar a acurácia do modelo, usamos o erro padrão residual (residual standard error):\n",
    "$$\n",
    "RSE = \\sqrt{\\frac{1}{n-2}\\sum_{i=1}^n (y_i-\\hat{y}_i)^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RSE: 0.8384690232980003\n"
     ]
    }
   ],
   "source": [
    "#funcao que calcula o RSE\n",
    "def RSE(x,y,b):\n",
    "    n = len(y)\n",
    "    RSE = 0\n",
    "    for i in range(0,n):\n",
    "        y_pred = b[0]+ x[i]*b[1] # valor predito\n",
    "        RSE = RSE + (y[i]-y_pred)**2\n",
    "    RSE = math.sqrt(RSE/(n-2))\n",
    "    return RSE\n",
    "print('RSE:', RSE(x,y,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outra medida importante é o coeficiente R2, que mede a proporção da variabilidade em Y que pode ser explicada a partir de X.\n",
    "$$\n",
    "R^2 = 1-\\frac{\\sum_{i=1}^n (y_i-\\hat{y}_i)^2}{\\sum_{i=1}^n(y_i-\\bar{y})^2}, \\quad 0\\leq R^2\\leq 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R2: 0.952538038613988\n"
     ]
    }
   ],
   "source": [
    "def R2(x,y,b):\n",
    "    n = len(y)\n",
    "    c1 = 0\n",
    "    c2 = 0\n",
    "    ym = np.mean(y)\n",
    "    for i in range(0,n):\n",
    "        y_pred = b[0]+ x[i]*b[1] # valor predito\n",
    "        c1 = c1 + (y[i]-y_pred)**2\n",
    "        c2 = c2 + (y[i]-ym)**2\n",
    "    R2 = 1 - c1/c2\n",
    "    return R2\n",
    "\n",
    "print('R2:', R2(x,y,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quanto mais próximo de um, melhor é o ajuste da regressão linear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O código completo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coeficientes:\n",
      "b_0 = 1.2363636363636363  \n",
      "b_1 = 1.1696969696969697\n",
      "R2: 0.952538038613988\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAclUlEQVR4nO3deZRV1ZXH8e8GowE1Di1tO4J22yIYBEVxDoIxThGMsyVBIqIMShxiJJqFI4oKCoIMMoW2cAIUpZFRcELBYqYoQGUWCEVUEAsZd/9xsAWkpKrecN997/dZK6uqLlX19npLftmce+7Z5u6IiEj8VIq6ABERqRgFuIhITCnARURiSgEuIhJTCnARkZjaJ50vdthhh3mNGjXS+ZIiIrE3bdq0te5ebffraQ3wGjVqUFBQkM6XFBGJPTNbuqfrWkIREYkpBbiISEwpwEVEYkoBLiISUwpwEZGY2muAm9kAM1tjZnN3uva0mc03s9lm9oaZHZzaMkVE4ik/H2rUgEqVwsf8/OT97rJ04IOAi3e7Ng442d3rAAuBDskrSUQkO+TnQ6tWsHQpuIePrVolL8T3GuDu/j7w1W7Xxrr71h1ffgIcnZxyRESyxwMPQEnJrtdKSsL1ZEjGGvifgHdK+0Mza2VmBWZWUFxcnISXExGJh2XLyne9vBIKcDN7ANgKlPoPAnfv6+713b1+tWo/eRJURCRrHXts+a6XV4UD3MyaA5cDea6xPiIiP/H441C16q7XqlYN15OhQgFuZhcDfwWucPeSvX2/iEguysuDvn2cVtWGU5ltVK8OffuG68lQlm2ELwMfAyea2QozuwXoARwIjDOzmWbWOznliIhkkaIi8l5sSJ/iq9j68ussWZK88IYynEbo7jfs4XL/5JUgIpJlNm6ETp2gc2c44AB48UW49tqkv0xaj5MVEcl648ZB69bwxRdw003QpQv8+7+n5KX0KL2ISDKsXg033ggXXRQeuxw/Hv7nf1IW3qAAFxFJzPbt0Ls31KwJw4ZBx44wezY0bpzyl9YSiohIRc2eDbffDh9/DBdcAL16wYknpu3l1YGLiJTXd9/BfffBqafCZ5/B4MEwYUJawxvUgYuIlM/IkdCuXTiZqmXLsNPk0EMjKUUduIhIWaxYAVddBb//Pey/P7z/ftgeGFF4gwJcROTnbdsG3brBSSfBqFFhf/eMGXDeeVFXpiUUEZFSTZsWDvCePh0uvhh69oTjj4+6qv+nDlxEZHfr10P79nDGGbByJbz6aui+Myi8QR24iMiP3GH4cLjzTli1Ctq0CUcHHnRQ1JXtkTpwERGAJUvCDcqrrw5PT37yCfTokbHhDQpwEcl1W7bAU09B7dowaVI4u+TTT8PySYbTEoqI5K6PP4bbboM5c6BJE+jePXnjctJAHbiI5J6vvw7BffbZ4fM33oA334xVeIMCXERyiTsMGRIOnurXD+66C+bNg6ZNo66sQrSEIiK54fPPw66ScePg9NNh9GioVy/qqhKiDlxEstumTfDoo3DyyTBlSthZ8vHHsQ9vUAcuItnsvffCca/z54eRZs8+C0ceGXVVSaMOXESyz9q10KIFNGwYOvB33glPU2ZReIMCXESyiTsMHBhuUr70Etx/P8ydG84xyUJaQhGR7FBUFJZL3n8fzjknjDk7+eSoq0opdeAiEm8bN8KDD8Ipp4QHcl58MYR4loc3qAMXkTgbNw5at4YvvoBmzeCZZ1I6BT7TqAMXkfhZvRpuvBEuuggqVw7zKAcPzqnwBgW4iMTJ9u1hbbtmTRg2DB56CGbNgkaNoq4sElpCEZF4mD07nF/yySdwwQXQq1fap8Bnmr124GY2wMzWmNncna4dambjzOyzHR8PSW2ZIpKzvvsO/vIXOPXU8Dj84MFhyWQv4Z2fDzVqQKVK4WN+flqqTauyLKEMAnbfRHk/MMHdTwAm7PhaRCS53n4batUKNydbtIAFC8LNSrOf/bH8/DDKcunSsDV86dLwdbaF+F4D3N3fB77a7XIT4B87Pv8HEM+jvEQkM61YAVddBVdcAQceCB98ELYHHnpomX78gQegpGTXayUl4Xo2qehNzMPdfRXAjo+l3vo1s1ZmVmBmBcXFxRV8ORHJCdu2QbducNJJYYhwp05hIvy555br1yxbVr7rcZXyXSju3tfd67t7/WrVqqX65UQkrgoKwhizP/85BHZhIXToAPvuW+5fVdpchpjNa9irigb4P83sCIAdH9ckryQRySnr14cp8A0awMqV4dCpUaPg+OMr/CsffxyqVt31WtWq4Xo2qWiAvwU03/F5c2BEcsoRkZzhDkOHhuWSHj3CE5U/HPu6l5uUe5OXB337QvXq4VdVrx6+zstLUu0ZYq/7wM3sZaAhcJiZrQA6Ak8Cr5nZLcAy4JpUFikiWWbJEmjbNnTadeuGmZRJngKfl5d9gb27vQa4u99Qyh81TnItIpLttmwJQxUeeihs0O7SJSyf7KNnCitC75qIpMfHH4cnKefMgSZNoHv37LurmGY6C0VEUuvrr0Nwn302fPMNvPlm+J/CO2EKcBFJDXcYMiQcPNW/P9x9N8ybF7pvSQotoYhI8n3+edhVMn58uDk5Zky4WSlJpQ5cRJJn0yZ49NEwDWfq1LA9cPJkhXeKqAMXkeSYNCnMpFywIOzlfvbZrJsCn2nUgYtIYtauhZtvDmd0b94M77wTnqZUeKecAlxEKsYdBgwI53Ln54dzS+bOhYt3P31aUkVLKCJSfkVFYbnk/ffDwVO9e0Pt2lFXlXPUgYtI2W3cCA8+CKecEh7I6dcP3ntvl/DOhUk4mUIduIiUzdix0KYNfPFFmIrzzDM/mQL/wyScH4Yp/DAJB7L/XJIoqAMXkZ+3ejXccAP87ndQuXKYRzl48E/CG3JnEk6mUICLyJ5t3x4mv9esCcOHhwOoZs+GRo1K/ZFcmYSTKRTgIvJTs2aFs0vatIHTTgvr3R07wn77/eyP5coknEyhABeRH23YAPfeG0J70aKwVDJ+PPz3f5fpx3NlEk6mUICLSPD222E3SZcu0KJFmI7TrFm5puPkyiScTKFdKCK5bsWKMFThjTdCgH/wQbmnwO8sFybhZAp14CK5autWeO65MJNy9Gh44gmYPj2h8Jb0UgcukosKCsIG7Rkz4JJLoGdPOO64qKuSclIHLpJL1q8PyyUNGoT93a+9Bv/7vwrvmFIHLpIL3GHYMGjfHlatChPhH3sMDjoo6sokAerARbLdkiVw+eVwzTVw+OEwZQo8/7zCOwsowEWy1ZYt0Lkz1KoVDpx69tkwJef006OuTJJESygi2Wjy5DAJfu5caNoUuneHY46JuipJMnXgItnk66/D7pJzzoF162DEiLC/W+GdlRTgItnAHV56KUzHGTAA7rkH5s2DK66IujJJIS2hiMTdwoXh0KkJE+CMM8K53ZoCnxMS6sDN7C4zKzSzuWb2spn9MlmFicRFZBNoNm2CRx6BOnXg00+Z2rwnx6+eTKVT62oSTo6ocAduZkcBdwK13H2jmb0GXA8MSlJtIhkvsgk0kyaFmZQLFsB11zHs3Gf541+P0CScHJPoGvg+QBUz2weoCqxMvCSR+Ej7BJriYmjeHC64IGwTHD0aXnmFe545QpNwclCFA9zdvwSeAZYBq4B17j529+8zs1ZmVmBmBcXFxRWvVCQDpW0Czfbt0L9/mI7z8svwt7+FLYK/+11665CMUuEAN7NDgCbAccCRwP5mdtPu3+fufd29vrvXr1atWsUrFclAaZlAM28eNGwILVuGh3JmzgwTEqpUSW8dknESWUK5EFjs7sXuvgUYDpydnLJE4iGlE2g2bgxrIHXrQmEh9OsXnqisVSu9dUjGSiTAlwFnmllVMzOgMVCUnLJE4iFlE2jGjIGTT4ZOncJE+Pnz4ZZbwlaXdNYhGc3cveI/bPYwcB2wFZgBtHT3TaV9f/369b2goKDCryeS9VatgrvugldfDXMoe/cONywlp5nZNHevv/v1hB7kcfeOQMdEfoeIEG5S9ukDHTrA99/Dww/DX/+61ynwktv0JKZI1GbNCgdPTZkCjRtDr15wwglRVyUxoLNQRKKyYQPcey+cdhosXhzOMhk3TuEtZaYOXCQKb70F7drB8uVw663w5JNw6KFRVyUxow5cJJ2WL4crr4QmTcJEnA8/DNtFFN5SAQpwkXTYujVMxKlVK2wRfPJJmD49nNstUkFaQhFJtU8/DTcpZ8yASy+FHj00BV6SQh24SKqsWwd33AENGsDq1fD66zBypMJbkkYduEiyucPQodC+fQjutm3hscc0BV6STgEukkyLF4fAfucdqFcvzKTUFHhJES2hSGxFNglnT7Zsgc6doXZt+OCDcMNy6lSFt6SUOnCJpcgm4ezJRx+Fm5SFhdC0KXTvrinwkhbqwCWW0j4JZ0+++ir8v8a558L69WG55I03FN6SNgpwiaVIJ9C4h8fea9aEAQPgnnvC0IUrrkjDi4v8SAEusRTZBJqFC+G3v4VmzeD442HaNHjmGTjggBS/sMhPKcAlltI+gWbTJnjkEahTBwoK4IUXwtr3Kaek6AVF9k4BLrGU1gk0EyeGoO7YMZxjMn8+tG4NlSun4MVEyk67UCS28vJSvOOkuDgc9zp4cFguGT36/6fAi2QCdeAiu9u+Hfr3DzcpX345bG2ZO1fhLRlHHbjIzgoL4fbbwzGv550XZlLuYQq8SCZQBy4CsHFj6LTr1g1bAvv3h0mTFN6S0dSBi4wZA23awKJF0Lw5PP00VKsWdVUie6UOXHLXqlVw/fVw8cXwi1+E3SaDBim8JTYU4JJ7tm0L+7hr1oQ33wz7u2fNgoYNo65MpFy0hCK5ZebMcPDU1KnQuDH06qUp8BJb6sAlN2zYEM4sqV8fliwJZ5mMG6fwllhTBy7Zb8SIMNps+fJweuCTT8Ihh0RdlUjC1IFL9lq+PDz63rRpGGf20UfQp4/CW7JGQgFuZgeb2VAzm29mRWZ2VrIKE6mwrVvDRJxatcIWwc6dYfp0OPvsqCsTSapEl1C6AaPd/Woz2xeourcfEEmpTz8NNylnzIBLL4WePcO8NZEsVOEO3Mx+BZwP9Adw983u/k2yChMpl3XroF07aNAA/vlPeP11GDlS4S1ZLZEllOOBYmCgmc0ws35mtv/u32RmrcyswMwKiouLE3g5kT1wD2F90klhb3e7dlBUBFdfHc6ZFcliiQT4PsCpQC93rwd8B9y/+ze5e193r+/u9avpCTdJpsWL4bLL4Npr4Ygjwt7u7t3hV7+KujKRtEgkwFcAK9x9yo6vhxICXSS1tmwJWwFr14YPPoDnnoMpU8Ieb5EcUuGbmO6+2syWm9mJ7r4AaAzMS15pInvw0UfhJmVhIfzhD9CtGxx9dNRViUQi0X3gdwD5ZjYbqAt0SrwkkT346iu49VY491z49lt4+20YNkzhLTktoW2E7j4T0L9bJXXcIT8f7r47hPi998JDD8H+P7lfLpJz9Ci9ZK6FC8M53RMmhO2B48ZpCrzITvQovWSeTZvg4Yfh17+GgoJwYuDkyQpvkd2oA5fMMnFimEm5cCHccAN07Qr/8R9RVyWSkdSBS2YoLoY//hEaNQpnmYwZA0OGKLxFfoYCXKK1fTv06wcnngivvBIGC8+dCxddFHVlIhlPSygSncLCsFzy4Ydw/vnQu3d4JF5EykQduKRfSQl06AB164ZzSwYMgEmTFN4i5aQOXNJr9OiwNXDxYrj5Znj6aTjssKirEokldeCSHqtWwXXXwSWXwL77ht0mAwcqvEUSoACXCsnPD0dtV6oUPubnl/KN27aFoQo1a4bZlI88ArNmQcOG6StWJEtpCUXKLT8/zAYuKQlfL10avgbIy9vpG2fODAdPTZ0KF14YzuvWFHiRpFEHLuX2wAM/hvcPSkrCdQA2bIB77gnHuy5ZEhJ/7FiFt0iSqQOXclu27GeujxgBd9wRJsLfdhs88YSmwIukiDpwKbdjj/3ptaNZzuhfNoWmTeHgg8PZJb17K7xFUkgBLuX2+ONQtWr4vDJbuYuuFHESjbaNhaeegmnT4Kyzoi1SJAdoCUXK7Ycbla/dO5WHVt9GPWbyZd3LOOCNHpoCL5JG6sCl/NatI29yW0b880zqHbkGhg7lqOlvK7xF0kwBLmXnDq+9FvZ09+4dblYWFcFVV4FZ1NWJ5BwtoUjZLFoEbduGR+FPPTXMpNQUeJFIqQOXn7d5c9gKWLt2ODXwuedgyhSFt0gGUAcupfvww3Dca2Eh/OEP0K2bpsCLZBB14PJTX30FLVvCeefBt9+G5ZJhwxTeIhlGAS4/cofBg8N0nEGD4C9/gXnz4PLLo65MRPZASygSLFgQzul+910480zo0wfq1Im6KhH5GerAc93338NDD4WwnjYNevWCjz5SeIvEgDrwXPbuu9C6NSxcCDfcAF27agq8SIyoA89Fa9ZAs2bQuHEYuDBmDAwZovAWiZmEA9zMKpvZDDMbmYyCJIW2b4d+/cKTlK++Cg8+CHPmwEUXRV2ZiFRAMpZQ2gNFwK+S8LskVQoLw/ncH30E558fHoXXFHiRWEuoAzezo4HLgH7JKUeSrqQEOnSAunVh/vwwSHjSJIW3SBZItAN/DrgPOLC0bzCzVkArgGP3NAlAUuedd8L5JYsXQ4sW4axuTYEXyRoV7sDN7HJgjbtP+7nvc/e+7l7f3etXq1atoi8n5bFyJVx7LVx6Key3X+i4BwxQeItkmUSWUM4BrjCzJcArQCMzeykpVUnFbNsGPXuG5ZG33oJHHw2T4X/zm6grE5EUqHCAu3sHdz/a3WsA1wPvuvtNSatMymfGjDDGrF07aNAA5s4Nu0z22y/qykQkRbQPPO42bIC77w7Huy5dGvZzjxkD//VfUVcmIimWlCcx3X0SMCkZv0vK4c03w1ScFSvCFsEnntAUeJEcog48jpYtgyZN4MorQ2BPnhz2dSu8RXKKAjxOtm6FLl2gVi0YPz5sC5w2Lax9i0jO0WFWcTFlSlgmmTUrnM/9/POaAi+S49SBZ7p168LDOGedBWvXhsk4b72l8BYRBXjGcg8HTtWsGda377wTiorCbEqzqKsTkQygJZRMtGhR6LpHj4bTToORI8NHEZGdqAPPJJs3Q6dOULt2ODWwe/ew9q3wFpE9UAeeKT78MNyknDcPrroKunWDo46KuioRyWDqwKP2r39By5Zw3nnw3XdhuWToUIW3iOyVAryM8vPDxo9KlcLH/PwEf6E7DB4cblIOGgT33ReGLlx2WeLFikhO0BJKGeTnQ6tWYTYChCNHWrUKn+flVeAXLlgQhglPnBi2B/bpA7/+ddLqFZHcoA68DB544Mfw/kFJSbheLt9/Dx07Qp064fTAPn3C2rfCW0QqQB14GSxbVr7rezRhQui6P/sMbrwRunaFww9PSn0ikpvUgZdBaZPgyjQhbs0aaNYMLrwwrHuPHRvWZBTeIpIgBXgZPP44VK2667WqVcP1Um3fDi++GG5Svvoq/P3vMHs2/Pa3Ka1VRHKHArwM8vKgb1+oXj08xV69evi61BuYc+fC+eeHO5116oQDqB55BKpUSWvdIpLdtAZeRnl5ZdhxUlISgrpLFzjoIBg4EJo319klIpISCvBkGTUqnF+yZAm0aBHO6tYUeBFJIS2hJGrlSrjmmvAATpUq8N57MGCAwltEUk4BXlHbtkGPHuEm5dtvw2OPwcyZYe1bRCQNtIRSEdOnh4OnCgrCrpIXXtAUeBFJO3Xg5fHtt3DXXXD66bB8OQwZAmPGKLxFJBLqwMvqzTfhjjvgyy9D9/3EE3DwwVFXJSI5TB343ixbBk2awJVXwqGHhkELvXopvEUkcgrw0mzdGvZz16oF48fD00+HNe+zzoq6MhERQEsoezZlSlgmmTULfv97eP758PiliEgGUQe+s2++gTZtQpe9di0MHw4jRii8RSQjVTjAzewYM5toZkVmVmhm7ZNZWFq5hwOnTjopnNF9551QVBTWvTPsMfikTwYSkdhKZAllK3CPu083swOBaWY2zt3nJam29Fi0KHTdY8aE6e8jR2bsFPikTwYSkVircAfu7qvcffqOz78FioD4TOLdvBk6dYLatWHyZOjePax9Z2h4QxInA4lIVkjKTUwzqwHUA6bs4c9aAa0Aji3TBIQ0+OADuP12mDcPrr4annsuFlPgkzIZSESyRsI3Mc3sAGAY8Gd3X7/7n7t7X3ev7+71q1WrlujLJeZf/4JbbgnnlXz3XVguef31WIQ3JDgZSESyTkIBbma/IIR3vrsPT05JKeAO//hHOHhq8GC47z4oLAwnCMZIhSYDiUjWSmQXigH9gSJ375q8kpJswQJo1AhuvhlOOCEcRNW5M+y/f9SVlVu5JwOJSFZLZA38HKAZMMfMZu649jd3H5V4WUnw/ffhvJInnwxtap8+0LJl2H8XY2WaDCQiOaHCAe7uHwKZtUn6BxMmQOvW8NlncOON0LWrpsCLSNaJdzu6uzVr4Kab4MILw7r32LFh87TCW0SyUHYE+PbtYTH4xBPhtdfg73+HOXPCsAURkSwV/8Os5s4NB09Nngy/+Q307h12m4iIZLn4duAlJXD//VCvXthpMmgQTJyo8BaRnBHPDnzUKGjbFpYsgT/9CZ56Cv7t36KuSkQkreLVga9cCddcEx7AqVIF3nsP+vdXeItITopHgG/bBj16hOWRkSPhscdg5szwSLyISI6KxxLKrbfCwIFw0UXwwgvwn/8ZdUUiIpGLR4C3bh3C+7rrMm7AgohIVDJ+CSU/H2pcczqVbryeGseZJtCIiOyQ0R24JtCIiJQuoztwTaARESldRge4JtCIiJQuowNcE2hEREqX0QGuCTQiIqXL6ADXBBoRkdJl9C4U0AQaEZHSZHQHLiIipVOAi4jElAJcRCSmFOAiIjGlABcRiSlz9/S9mFkxsLSCP34YsDaJ5cSd3o8f6b3Yld6PXWXD+1Hd3avtfjGtAZ4IMytw9/pR15Ep9H78SO/FrvR+7Cqb3w8toYiIxJQCXEQkpuIU4H2jLiDD6P34kd6LXen92FXWvh+xWQMXEZFdxakDFxGRnSjARURiKhYBbmYXm9kCM/vczO6Pup6omNkxZjbRzIrMrNDM2kddUyYws8pmNsPMRkZdS9TM7GAzG2pm83f8d3JW1DVFxczu2vH3ZK6ZvWxmv4y6pmTL+AA3s8pAT+ASoBZwg5nViraqyGwF7nH3k4AzgbY5/F7srD1QFHURGaIbMNrdawKnkKPvi5kdBdwJ1Hf3k4HKwPXRVpV8GR/gwBnA5+6+yN03A68ATSKuKRLuvsrdp+/4/FvCX86joq0qWmZ2NHAZ0C/qWqJmZr8Czgf6A7j7Znf/JtqqIrUPUMXM9gGqAisjrifp4hDgRwHLd/p6BTkeWgBmVgOoB0yJtpLIPQfcB2yPupAMcDxQDAzcsaTUz8z2j7qoKLj7l8AzwDJgFbDO3cdGW1XyxSHAbQ/Xcnrvo5kdAAwD/uzu66OuJypmdjmwxt2nRV1LhtgHOBXo5e71gO+AnLxnZGaHEP6lfhxwJLC/md0UbVXJF4cAXwEcs9PXR5OF/xQqKzP7BSG88919eNT1ROwc4AozW0JYWmtkZi9FW1KkVgAr3P2Hf5UNJQR6LroQWOzuxe6+BRgOnB1xTUkXhwD/FDjBzI4zs30JNyLeirimSJiZEdY3i9y9a9T1RM3dO7j70e5eg/DfxbvunnVdVlm5+2pguZmduONSY2BehCVFaRlwpplV3fH3pjFZeEM344cau/tWM2sHjCHcSR7g7oURlxWVc4BmwBwzm7nj2t/cfVSENUlmuQPI39HsLAJaRFxPJNx9ipkNBaYTdm/NIAsfqdej9CIiMRWHJRQREdkDBbiISEwpwEVEYkoBLiISUwpwEZGYUoCLiMSUAlxEJKb+Dw/HpzxHt39XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from statistics import variance \n",
    "import math \n",
    "\n",
    "# define os dados\n",
    "x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) \n",
    "y = np.array([1, 3, 2, 5, 7, 8, 8, 9, 10, 12]) \n",
    "\n",
    "# funcao que estima os coeficientes\n",
    "def estimate_coef(x, y): \n",
    "    # número de observações\n",
    "    n = np.size(x) \n",
    "    # estimadores\n",
    "    m_x, m_y = np.mean(x), np.mean(y) \n",
    "    SS_xy = np.sum(y*x) - n*m_y*m_x \n",
    "    SS_xx = np.sum(x*x) - n*m_x*m_x \n",
    "    # calcula os coeficientes de regressão\n",
    "    b_1 = SS_xy / SS_xx \n",
    "    b_0 = m_y - b_1*m_x \n",
    "    return(b_0, b_1) \n",
    "\n",
    "# função para mostrar os dados e o ajuste linear\n",
    "def plot_regression_line(x, y, b): \n",
    "    # mostra os dados\n",
    "    plt.scatter(x, y, color = \"b\", marker = \"o\", s = 50) \n",
    "    # prediz os valores\n",
    "    y_pred = b[0] + b[1]*x \n",
    "    # mostra a reta de regressão\n",
    "    plt.plot(x, y_pred, color = \"r\") \n",
    "    plt.xlabel('x', fontsize = 15) \n",
    "    plt.ylabel('y', fontsize = 15) \n",
    "    plt.show(True) \n",
    "    \n",
    "def R2(x,y,b):\n",
    "    n = len(y)\n",
    "    c1 = 0\n",
    "    c2 = 0\n",
    "    ym = np.mean(y)\n",
    "    for i in range(0,n):\n",
    "        y_pred = b[0]+ x[i]*b[1] # valor predito\n",
    "        c1 = c1 + (y[i]-y_pred)**2\n",
    "        c2 = c2 + (y[i]-ym)**2\n",
    "    R2 = 1 - c1/c2\n",
    "    return R2\n",
    "    \n",
    "# estima os coeficientes\n",
    "b = estimate_coef(x, y) \n",
    "print(\"Coeficientes:\\nb_0 = {}  \\nb_1 = {}\".format(b[0], b[1])) \n",
    "print('R2:', R2(x,y,b))\n",
    "\n",
    "\n",
    "# mostra o ajuste linear\n",
    "plt.plot(x,y,'bo')\n",
    "plt.plot(x, b[0] + b[1]*x, 'r-')\n",
    "plt.savefig('plot.eps')"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}