2ex-Exercicio2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Exercício 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Revisitando o ajuste linear"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "No exercício 1 fizemos uma regressão linear simples para resolver o problema de ajustar de reta para um dado conjunto de dados e uma teoria, \n",
    "dadoum conjunto de dados $\\{(x_1,y_1),(x_2,y_2),...,(x_n,y_n)\\}$ para o qual a teoria diz que a relação entre $x$ e $y$ é dada por:\n",
    "\n",
    "$$\n",
    "y_i = a + bx_i\n",
    "$$\n",
    "\n",
    "Sabendo que as medidas estão sujeitas a várias formas de imprecisão e que a teoria não conta com todas tais contribuições,\n",
    "consideremos o modelo pára os dados\n",
    "\n",
    "$$\n",
    "y_i = a + bx_i + \\xi_i\n",
    "$$\n",
    "\n",
    "onde $\\xi_i$ são valores de ruído que assumimos ser Gausiano.\n",
    "\n",
    "Os valores de $a$ e $b$ que melhor descrevem essa teoria devem ser aqueles que maximizam a verossimilhança de obsevar\n",
    "um certo valor de ruído, ou seja, minimizam $\\sum_i\\xi_i^2 = \\sum_i(y_i - a - bx_i)^2$.\n",
    "\n",
    "A solução desse problema é dada por:\n",
    "$$\n",
    "\\hat{a} = \\bar{y} - \\hat{b}\\bar{x}\\\\\n",
    "\\hat{b} = \\frac{\\frac{1}{n}\\sum_i(x_i-\\bar{x})(y_i-\\bar{y})}{\\frac{1}{n}\\sum_i(x_i-\\bar{x})^2}\n",
    "$$\n",
    "\n",
    "onde\n",
    "$$\n",
    "\\bar{x} = \\frac{1}{n}\\sum_i x_i\\\\\n",
    "\\bar{y} = \\frac{1}{n}\\sum_i y_i\n",
    "$$\n",
    "são as médias amostrais de $x$ e $y$\n",
    "\n",
    "Agora que temos e scipy à disposição, refaça o exercíocio usando quaiquer ferramentas disponíveis nesses módulos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-20T18:18:36.613220Z",
     "start_time": "2017-07-20T18:18:36.599927Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "from numpy.random import randn, seed\n",
    "seed(42)\n",
    "a_real, b_real = 1.75, 3.5\n",
    "x = np.arange(-10,11)\n",
    "y = a_real + b_real * x + randn(x.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a_real = 1.75    a_chapeu = 1.6566513114326173\n",
      "b_real = 3.5    b_chapeu = 3.4352713112352835\n"
     ]
    }
   ],
   "source": [
    "b_chapeu = ((x - x.mean())*(y - y.mean())).mean()/x.var()\n",
    "a_chapeu = y.mean() - x.mean()*b_chapeu\n",
    "print(f\"a_real = {a_real}\",f\"a_chapeu = {a_chapeu}\", sep='    ')\n",
    "print(f\"b_real = {b_real}\", f\"b_chapeu = { b_chapeu}\", sep='    ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dessa vez os dados já estão nas variáveis adequadas, `x` e `y`.\n",
    "Tente encontrar o modo mais fácil de resover esse problema."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Encontrando os zeros de uma função"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def calcula(func, *x):\n",
    "    return func(*x)\n",
    "\n",
    "calcula(lambda x,y: x+y,2,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para funções bem comportadas de uma variável, podemos encontrar os valores onde ela se anula usando o método de Newton.\n",
    "Sejam $f:\\mathbb{R}\\rightarrow\\mathbb{R}$ uma função com pelo menos a primeira derivada contínua e $x_0$ um chute inicial de onde a função se anula, podemos usar a tangente à curva da função no ponto $x_0$ para chegar mais próximo do ponto onde a função se anula fazendo \n",
    "\n",
    "$$\n",
    "x_1 = x_0 - \\frac{f(x_0)}{f'(x_0)}\n",
    "$$\n",
    "\n",
    "Caso o chute inicial seja bom, iterando esse algoritmo podemos nos aproximar cada vez mais do zero da função fazendo\n",
    "\n",
    "$$\n",
    "x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_n)}\n",
    "$$\n",
    "\n",
    "Usando a seus conhecimentos de Python e numpy, encontre os zeros das seguintes funções:\n",
    "$f_1(x) = \\exp(x) - 1$, $f_2(x) = x\\ln(x) - x$ e $f_3(x) = \\sin(x)$ (cosidere $x$ no intervalo $(0,2\\pi)$ no caso de $f_3$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton(func, dfunc, x0, *, convergiu=0.01, nmax=100):\n",
    "    xn = x0 - func(x0)/dfunc(x0)\n",
    "    while (abs(xn - x0) > convergiu) and nmax > 0:\n",
    "        xn -= func(xn)/dfunc(xn)\n",
    "        nmax -= 1\n",
    "    return xn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.908684158558477e-17"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newton(lambda x: np.exp(x) - 1, lambda x: np.exp(x), -2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.7182818284590455"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newton(lambda x: x*np.log(x) - x, lambda x: np.log(x), 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589793"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newton(lambda x: np.sin(x), lambda x: np.cos(x), 1.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.337411025717124e-18"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.optimize as spo\n",
    "\n",
    "spo.newton(lambda x: np.exp(x) - 1, -2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Calculando o valor de pi usando Monte Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O exemplo mais simples de Monte Carlo é o cálculo do valor de $\\pi$, usando pontos aleatórios dentro do quadrado $x \\in [-1,1], y\\in[-1,1]$ e calculando a razão entre todos os números sorteados e aqueles que caem dentro do círculo.\n",
    "Calcule o valor de $\\pi$ usando o menor número de linhas de código possível com numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y = 2*np.random.rand(2,5) - 1\n",
    "dentro = x**2 + y**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.32650532, 1.03266702])"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dentro[dentro > 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4*dentro[dentro <= 1].shape[0]/x.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.1417648"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y = 2*np.random.rand(2,10000000) - 1\n",
    "\n",
    "dentro = (x**2 + y**2)\n",
    "4*dentro[dentro <= 1].shape[0]/x.shape[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "48px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}