"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def genera_líneas(numero_lineas, numero_puntos):\n",
" a = 5\n",
" b = 1\n",
" X = np.linspace(0, 100, numero_puntos)\n",
" Y = np.empty((len(X), numero_lineas))\n",
" for i in range(Y.shape[1]):\n",
" eps = np.random.normal(scale = 100, size = numero_puntos)\n",
" Y[:,i] = a * X + b + eps\n",
" ahat, bhat = np.polyfit(X, Y[:,i], deg = 1)\n",
" plt.plot(X, X * ahat + bhat, color = (0,0,0.9))\n",
" plt.plot(X, X * a + b, 'k', lw = 10)\n",
" ahat, bhat = np.polyfit(X, Y.mean(axis = 1), deg = 1)\n",
" plt.plot(X, X * ahat + bhat, 'y', lw = 5)\n",
" plt.text(20, 400, 'ahat = {:3.2f}'.format(ahat))\n",
" plt.text(20, 375, 'bhat = {:3.2f}'.format(bhat))\n",
" plt.xlim(0,100)\n",
" plt.ylim(-50,550)\n",
" \n",
"interact(genera_líneas, \n",
" numero_lineas = widgets.FloatSliderWidget(min = 3, \n",
" max = 100, \n",
" step = 1, \n",
" value = 3, \n",
" description = \"Número de líneas\"), \n",
" numero_puntos = widgets.FloatSliderWidget(min = 5, \n",
" max = 100, \n",
" step = 1, \n",
" value = 5, description=\"Puntos para cada línea\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Las diferencias que vemos en el anterior gráfico entre la línea negra y el resto de líneas es debido a que usamos una muestra (líneas azules) para intentar obtener características de una población más grande (línea negra). Si la muestra de datos es más grande o si usamos un mayor número de líneas vemos que, en promedio (línea amarilla), nos vamos acercando cada vez más al valor esperado. Como hemos comentado en el párrafo anterior, hay veces que las líneas pueden quedar por encima y a veces pueden quedar por debajo en mayor o menor proporción pero, en conjunto, vemos que a medida que crece el número de datos para cada línea o el número de líneas, nos vamos acercando a un valor promedio. Esto nos está indicando que en el caso de la estimación de los parámetros mediante el método de mínimos cuadrados es no sesgada ([unbiased](https://en.wikipedia.org/wiki/Bias_of_an_estimator))."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Muy bien, pero normalmente solo tenemos un conjunto de datos y solo podemos hacer una estimación con ese número de datos. En determinados rangos podemos estar subreestimando mientras que en otros rangos podemos estar infraestimando (o viceversa). \n",
"\n",
"¿Cómo podemos estimar el error de la estimación?\n",
"\n",
"Para ello podemos usar el **error estándar** ([standard error](https://en.wikipedia.org/wiki/Standard_error)) usando las siguientes fórmulas para el caso de una regresión lineal simple:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$SE(\\hat{a})^2 = \\frac{\\sigma^2}{\\sum_{i=1}^{n}(x_{i}-\\bar{x})^2}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$SE(\\hat{b})^2 = \\sigma^2[\\frac{1}{n} + \\frac{\\bar{x}^2}{\\sum_{i=1}^{n}(x_{i}-\\bar{x})^2}]$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"siendo $\\sigma^2=Var(\\epsilon)$. \n",
"\n",
"Para que las fórmulas anteriores sean estrictamente válidas necesitamos que los errores $\\epsilon_{i}$ para cada observación no estén correlacionados con la varianza $\\sigma^2$. Normalmente podemos encontrar cierta correlación pero las fórmulas siguen siendo una buena aproximación.\n",
"\n",
"En general, $\\sigma^2$ no se conoce pero lo podemos obtener a partir de los datos. Su estimación se conoce como error estándar residual (residual standard error) y se puede obtener a partir de $RSE (también\\ conocido\\ como\\ \\hat{\\sigma}) = \\sqrt{RSS / (n - 2)}$. Siendo rigurosos, cuando $\\sigma^2$ se obtiene a partir de la muestra deberiamos escribir el error estándar con una caperuza, e.g. $\\hat{SE}(\\hat{a})$, para indicar que es una estimación. Si en algún momento se nos olvida poner la caperuza tened en cuenta que nos referimos a una estimación."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Intervalos de confianza"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Los errores estándar se pueden usar para calcular intervalos de confianza. Un intervalo de confianza, por ejemplo al 95%, se define como el rango de valores que contendrá el valor verdadero desconocido del parámetro con una probabilidad del 95%. El rango se calcula a partir de los datos de la muestra y contiene un valor superior y un valor inferior. El intervalo de confianza al 95% para el parámetro $a$ en la regresión lineal se podría obtener, de forma aproximada (\\*), usando:"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"$$\\hat{a} \\pm 2 · SE(\\hat{a})$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(*) En realidad, el 2 de la fórmula anterior tiende al valor 1.96 cuando el número de grados de libertad aumenta, para el caso de un intervalo de confianza del 95%. Como nuestro número de grados de libertad es bastante inferior (n - 2 = 9) el valor de la fórmula correcto a usar, en este caso concreto, estaría alrededor de 2.262.
\n",
"\n",
"Ese valor lo podemos obtener de la siguiente forma:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"2.262157162740992"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 0.025 es 0.05 (5%) / 2 para tener en cuenta ambas colas\n",
"stats.t.isf(0.025, len(X)-2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Por tanto, para obtener la probabilidad de un 95% de que el valor de $a$ esté incluido dentro de determinado intervalo usaremos la formula $[\\hat{a} - 2.262 · SE(\\hat{a}), \\hat{a} + 2.262 · SE(\\hat{a})]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Por tanto, si cogemos los datos del cuarteto de Anscombe que estamos usando tendremos que:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Probabilidad del 95% de que a esté en el intervalo : 0.233388666744 0.766793151438\n",
"Probabilidad del 95% de que b esté en el intervalo : 0.455913668282 5.5442681499\n"
]
}
],
"source": [
"# Calculamos los parámetros del ajuste, a y b, ax+b\n",
"ahat, bhat = np.polyfit(X, Y, deg = 1)\n",
"# Calculamos los residuos\n",
"residuos = Y - (ahat * X + bhat)\n",
"# Calculamos la suma cuadrática de los residuos (RSS, residual sum of squares)\n",
"RSS = np.sum(residuos**2)\n",
"# Calculamos el error estándar de los residuos (RSE, residual standard error)\n",
"RSE = np.sqrt(RSS/(len(X)-2))\n",
"# Error estándar de la estimación de la pendiente\n",
"SEa = np.sqrt(RSE**2 / np.sum((X - X.mean())**2))\n",
"SEb = np.sqrt(RSE**2 * (1 / len(X) + X.mean()**2 / np.sum((X - X.mean())**2)))\n",
"\n",
"print(\"Probabilidad del 95% de que a esté en el intervalo :\", ahat - 2.262 * SEa, ahat + 2.262 * SEa)\n",
"print(\"Probabilidad del 95% de que b esté en el intervalo :\", bhat - 2.262 * SEb, bhat + 2.262 * SEb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En el anterior gráfico hemos mostrado las rectas teniendo en cuenta todos los extremos del rango tanto en $\\hat{a}$ como en $\\hat{b}$."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Test de hipótesis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Los errores estándar se pueden usar, también, para comprobar hipótesis en los coeficientes. La prueba más evidente que se nos puede ocurrir sería conocer si existe relación entre $X$ e $Y$.. Vamos a montar nuestra prueba.\n",
"\n",
"* La hipótesis nula, $H_{0}$, sería considerar que no hay relación entre $X$ e $Y$.\n",
"* Mientras que la hipótesis alternativa, $H_{a}$, sería considerar que existe esa relación.\n",
"\n",
"Matemáticamente, habría que comprobar si la pendiente de la recta es nula y por tanto la relación se queda en $Y = b + \\epsilon$. Por tanto, la prueba se podría representar como:\n",
"\n",
"$$H_{0}: a = 0$$\n",
"$$H_{a}: a \\neq 0$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Para comprobar la hipótesis nula tendríamos que determinar que nuestra estimación de la pendiente, $\\hat{a}$, está lo sucientemente lejos de cero de forma que podemos estar seguros que $a$ no es cero. Pero, ¿cuán lejos sería suficiente? Esto, por supuesto, depende de la precisión en el cálculo de $\\hat{a}$ que dependerá del error estándar, $SE(\\hat{a})$. Si el error estándar es pequeño, incluso valores pequeños de la estimación de la pendiente, $\\hat{a}$ mostrarían evidencias de que existe relación entre $X$ e $Y$. Sin embargo, si el error estándar es grande, necesitariamos valores grandes de $\\hat{a}$ para poder asegurar que existe tal relación entre $X$ e $Y$. En la práctica, se calcula un estadístico t de Student con la siguiente forma:\n",
"\n",
"$$t = \\frac{\\hat{a} - 0}{SE(\\hat{a})}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El valor cero en la anterior ecuación corresponde al valor que consideremos en la hipótesis, en este caso estamos comprobando que no haya relación entre las variables $X$ e $Y$ y eso se cumple considerando que la pendiente es igual a $0$, i.e. $a = 0$.\n",
"\n",
"La ecuación anterior indica el número de desviaciones estándar que nuestra estimación de la pendiente está alejada de, en este caso, $0$. El valor crítico para considerar que aceptamos o rechazamos la hipótesis vendrá dado por una distribución t de Student de $n - 2$ grados de libertad. Con todo esto, es muy simple calcular la probabilidad de observar valores iguales o superiores a $|t|$ asumiendo $a = 0$. Esta probabilidad se conoce como *p-value*, que se puede interpretar como:\n",
"\n",
"* un valor bajo de $p$ indica que es improbable que la asociación entre el predictor y la respuesta sea debida a la casualidad (siempre que no haya asociación real entre $X$ e $Y$).\n",
"\n",
"Por tanto, si observamos un valor de $p$ bajo podemos inferir que existe asociación entre el predictor y la respuesta y podemos declarar que existe relación entre $X$ e $Y$ si el valor de $p$ es lo suficientemente bajo y podemos, por tanto, rechazar la hipótesis nula, $H_{0}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Veamos un ejemplo de todo esto con código:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Valor de la pendiente estimada, a = 0.500090909091\n",
"RSS = 13.76269\n",
"RSE = 1.23660332273\n",
"Error estándar = 0.117905500596\n",
"Valor de t = 4.24145528889\n",
"Valor de p = 0.00216962887308\n"
]
}
],
"source": [
"# X e Y serán el primer conjunto del cuarteto de Anscombe \n",
"# (https://en.wikipedia.org/wiki/Anscombe%27s_quartet)\n",
"X = np.array([10.0, 8.0, 13.0, 9.0, 11.0, 14.0, 6.0, 4.0, 12.0, 7.0, 5.0])\n",
"Y = np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])\n",
"\n",
"# Calculamos los parámetros del ajuste, a y b, ax+b\n",
"ahat, bhat = np.polyfit(X, Y, deg = 1)\n",
"# Calculamos los residuos\n",
"residuos = Y - (ahat * X + bhat)\n",
"# Calculamos la suma cuadrática de los residuos (RSS, residual sum of squares)\n",
"RSS = np.sum(residuos**2)\n",
"# Calculamos el error estándar de los residuos (RSE, residual standard error)\n",
"RSE = np.sqrt(RSS/(len(X)-2))\n",
"# Error estándar de la estimación de la pendiente\n",
"SEa = np.sqrt(RSE**2 / np.sum((X - X.mean())**2))\n",
"# Calculamos el estadístico t\n",
"t = ahat / SEa\n",
"# Calculamos el p-valor\n",
"# Aquí estamos usando scipy.stats \n",
"# donde metemos los grados de libertad (len(X) - 2) \n",
"# y calculamos la CDF del valor t\n",
"# Como calculamos la CDF y solo obtenemos una de las \n",
"# colas lo multiplicamos por dos para considerar ambos extremos\n",
"p = 2 * (1 - stats.t(len(X) - 2).cdf(t))\n",
"\n",
"print('Valor de la pendiente estimada, a = ', ahat)\n",
"print('RSS = ', RSS)\n",
"print('RSE = ', RSE)\n",
"print('Error estándar = ', SEa)\n",
"print('Valor de t = ', t)\n",
"print('Valor de p = ', p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A partir de los valores anteriores podemos construir la siguiente tabla resumen:\n",
"\n",
"\n",
" \n",
" \n",
" | \n",
" Coeficiente | \n",
" Error Estándar | \n",
" estadístico t de Student | \n",
" valor - p | \n",
"
\n",
" \n",
" \n",
" \n",
" | Pendiente | \n",
" 0.500 | \n",
" 0.117 | \n",
" 4.241 | \n",
" 0.002 | \n",
"
\n",
" \n",
"
\n",
"\n",
"En la tabla anterior vemos que el valor-p está lejos de y es inferior a, por ejemplo, 0.05 que sería lo que indicaría un intervalo de confianza del 95%. El intervalo de confianza al 95% de la pendiente que calculamos unas cuantas líneas más arriba indicaba el siguiente rango $[0.233, 0.767]$ por lo que parece que podríamos rechazar la hipótesis nula y aceptar la hipótesis alternativa de que existe relación entre $X$ e $Y$ con cierta confianza.\n",
"\n",
"A pesar de lo dicho anteriormente, tened en cuenta que las pruebas de hipótesis deberían guiar nuestra toma de decisiones, pero no hay que tomarlos de forma estricta ni confiar ciegamente en los mismos si tenemos ciertas incertidumbres (pocos datos, modelo poco preciso,...).
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Estimando la precisión del modelo (bondad del ajuste o *goodness of fit*)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La calidad de un ajuste por regresión lineal se realiza usando, normalmente, dos cantidades, el $RSE$ (Residual Standard Error o error estándar de los residuos, visto brevemente más arriba) y el estadístico $R^2$.\n",
"\n",
"#### RSE\n",
"\n",
"El $RSE$ es una estimación de la desviación estándar del término de error, $\\epsilon$. En otras palabras, es la cantidad promedio que la respuesta se desviará de la línea real de regresión. Como vimos anteriormente, se calcula usando la siguiente fórmula:\n",
"\n",
"$$RSE = \\sqrt{RSS / (n - 2)}$$\n",
"\n",
"donde RSS (Residual Sum of Squares) se expresa como $RSS = \\sum_{i=1}^{n}{r_i}^2 = \\sum_{i=1}^{n}({y_{i} - \\hat{y_{i}}})^2$.\n",
"\n",
"En el ejemplo anterior hemos visto que el $RSE$ valía $1.24$, esto indica que cualquier predicción que hagamos de $Y$ basándonos en $X$ tendrá un error de $1.24$ unidades. Este error será aceptable o no dependiendo del problema en cuestión. Para conocer el error relativo podemos usar la siguiente fórmula:\n",
"\n",
"$$Error\\ relativo = (\\frac{RSE}{\\bar{y}}) · 100$$\n",
"\n",
"Y en el ejemplo que estamos usando, el valor sería $16.5\\%$, aproximadamente.\n",
"\n",
"Como el $RSE$ se mide en las unidades de $Y$ no siempre puede parecer obvio si el resultado nos indica que el ajuste es bueno. Por ello es interesante obtener también el error relativo.\n",
"\n",
"Debería parecer obvio que el valor $RSE$ debería ser siempre lo más bajo posible dentro de su contexto (bajo error relativo) para considerar el ajuste bueno a partir de este estadístico.\n",
"\n",
"#### R2\n",
"\n",
"$R^2$ proporciona una forma alternativa de calcular la bondad del ajuste. Toma forma de proporción, la proporción de varianza explicada, y su valor se restringe al intervalo [0,1] y es independiente de las unidades de $Y$.\n",
"\n",
"Para calcular $R^2$ usamos la fórmula:\n",
"\n",
"$$R^2 = \\frac{TSS-RSS}{TSS} = 1 - \\frac{RSS}{TSS}$$\n",
"\n",
"donde $TSS$, la suma total de los cuadrados o *Total Sum of Squares*, mide la varianza total en la respuesta y se calcula como $TSS = \\sum_{i=1}^{n}({y_{i} - \\bar{y}})^2$. \n",
"\n",
"La fórmula de $R^2$ se puede ver como ( 1 - (cantidad de variabilidad que queda sin explicar después del cálculo de la regresión) / (cantidad de variabilidad inherente en la respuesta antes del cálculo de la regresión)) o, dicho de otra forma, $R^2$ mide la proporción de variabilidad en $Y$ que puede ser explicada usando $X$.\n",
"\n",
"Si el valor de $R^2$ es alto, cercano a 1, indicará que una gran proporción de la variabilidad de la respuesta ha sido explicada por la regresión. Si, por otra parte, el valor es cercano a cero tenemos que la regresión no es capaz de explicar la variabilidad de la respuesta indicándonos que el modelo es erróneo, el error inherente, $\\sigma^2$, es alto o que se combinan ambos problemas.\n",
"\n",
"Calculemos el valor de $R^2$ con código:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.666542459509\n"
]
}
],
"source": [
"# Calculamos ls suma total de los cuadrados\n",
"TSS = np.sum((Y - Y.mean())**2)\n",
"# Calculamos los residuos\n",
"residuos = Y - (ahat * X + bhat)\n",
"# Calculamos la suma cuadrática de los residuos (RSS, residual sum of squares)\n",
"RSS = np.sum(residuos**2)\n",
"\n",
"R2 = 1 - RSS/TSS\n",
"\n",
"print(R2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vemos que $X$ sería capaz de explicar 2/3 de la variabilidad en $Y$.\n",
"\n",
"¿Cuándo será bueno el valor de $R^2$? Normalmente, querremos que el valor sea lo más cercano posible a 1 pero también dependerá del problema concreto para decidir si un valor más bajo lo consideramos como adecuado.\n",
"\n",
"$R^2$ es una medida de la relación lineal entre $X$ e $Y$. El coeficiente de correlación, definido como:\n",
"\n",
"$$Cor(X,Y) = \\frac{\\sum_{i=1}^{n}({x_{i} - \\bar{x}})({y_{i} - \\bar{y}})}{\\sqrt{\\sum_{i=1}^{n}({x_{i} - \\bar{x}})^2}\\sqrt{\\sum_{i=1}^{n}({y_{i} - \\bar{y}})^2}}$$\n",
"\n",
"también es una medida de la relación lineal entre $X$ e $Y$. Podríamos usar $r = Cor(X,Y)$ como medida de la bondad del ajuste en lugar de $R^2$. De hecho, en el caso concreto de la regresión lineal simple, $r^2 = R^2$. En el caso de la regresión lineal múltiple no podemos extender este hecho ya que la correlación mide la relación entre una variable y una respuesta pero no entre múltiples variables y una respuesta. El trabajo de obtener la relación entre varias variables corresponderá a $R^2$ en el caso de la regresión lineal múltiple.\n",
"\n",
"Comprobemos que $r^2 = R^2$:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.666542459509 0.666542459509\n",
"None\n"
]
}
],
"source": [
"r2 = np.corrcoef(X, Y)[0,1]**2\n",
"print(r2, R2)\n",
"print(np.testing.assert_almost_equal(r2, R2, decimal = 7))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Regresión lineal simple con statsmodels"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Statsmodels es la librería de facto para hacer modelos y tests estadísticos. Aumenta enormemte las capacidades estadísticas nu Numpy y Scipy y se complementa muy bien con Pandas y scikit-learn.\n",
"\n",
"Todo lo que hemos hecho hasta ahora se podría hacer de la siguiente forma:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0001 1.125 2.667 0.026 0.456 5.544\n",
"x1 0.5001 0.118 4.241 0.002 0.233 0.767\n",
"==============================================================================\n",
"R^2 = 0.666542459509\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"d:\\users\\X003621\\AppData\\Local\\Continuum\\Miniconda3\\lib\\site-packages\\scipy\\stats\\stats.py:1233: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=11\n",
" int(n))\n"
]
}
],
"source": [
"result = sm.OLS(Y, sm.add_constant(X)).fit()\n",
"print(result.summary().tables[1])\n",
"print('R^2 = ', result.rsquared)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ¿Es nuestro modelo bueno?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recordad la cita que hemos visto más arriba. No os enamoréis del modelo. \n",
"\n",
"Regresemos de nuevo al cuarteto de Anscombe donde cuatro grupos de datos presentan el mismo valor de $R^2$ y los mismos parámetros de la recta de ajuste."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"ans_x_I = np.array([10.0, 8.0, 13.0, 9.0, 11.0,\n",
" 14.0, 6.0, 4.0, 12.0, 7.0, 5.0])\n",
"ans_y_I = np.array([8.04, 6.95, 7.58, 8.81, 8.33,\n",
" 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])\n",
"ans_x_II = np.array([10.0, 8.0, 13.0, 9.0, 11.0,\n",
" 14.0, 6.0, 4.0, 12.0, 7.0, 5.0])\n",
"ans_y_II = np.array([9.14, 8.14, 8.74, 8.77, 9.26,\n",
" 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])\n",
"ans_x_III = np.array([10.0, 8.0, 13.0, 9.0, 11.0,\n",
" 14.0, 6.0, 4.0, 12.0, 7.0, 5.0])\n",
"ans_y_III = np.array([7.46, 6.77, 12.74, 7.11, 7.81,\n",
" 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])\n",
"ans_x_IV = np.array([8.0, 8.0, 8.0, 8.0, 8.0,\n",
" 8.0, 8.0, 19.0, 8.0, 8.0, 8.0])\n",
"ans_y_IV = np.array([6.58, 5.76, 7.71, 8.84, 8.47,\n",
" 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])\n",
"x = [ans_x_I, ans_x_II, ans_x_III, ans_x_IV]\n",
"y = [ans_y_I, ans_y_II, ans_y_III, ans_y_IV]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Y usemos statsmodels para obtener más información del modelo, a parte del valor de $R^2$"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0001 1.125 2.667 0.026 0.456 5.544\n",
"x1 0.5001 0.118 4.241 0.002 0.233 0.767\n",
"==============================================================================\n",
"0.666542459509\n",
"\n",
"\n",
"\n",
"\n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0009 1.125 2.667 0.026 0.455 5.547\n",
"x1 0.5000 0.118 4.239 0.002 0.233 0.767\n",
"==============================================================================\n",
"0.666242033727\n",
"\n",
"\n",
"\n",
"\n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0025 1.124 2.670 0.026 0.459 5.546\n",
"x1 0.4997 0.118 4.239 0.002 0.233 0.766\n",
"==============================================================================\n",
"0.666324041067\n",
"\n",
"\n",
"\n",
"\n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0017 1.124 2.671 0.026 0.459 5.544\n",
"x1 0.4999 0.118 4.243 0.002 0.233 0.766\n",
"==============================================================================\n",
"0.666707256898\n",
"\n",
"\n",
"\n",
"\n"
]
}
],
"source": [
"for X, Y in zip(x, y):\n",
" result = sm.OLS(Y, sm.add_constant(X)).fit()\n",
" print(result.summary().tables[1])\n",
" print(result.rsquared)\n",
" print('\\n'*3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vaya, parece que los cuatro conjuntos de datos dan resultados similares para diferentes estadísticos. Los coeficientes del modelo son similares, los errores estándar también e incluso los intervalos de confianza salen parecidos."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"¿Podemos afirmar categóricamente que el ajuste lineal es el más adecuado para todos los casos? Un primer paso debería ser el análisis exploratorio de los datos:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.667\n",
"Model: OLS Adj. R-squared: 0.629\n",
"Method: Least Squares F-statistic: 17.99\n",
"Date: Thu, 13 Aug 2015 Prob (F-statistic): 0.00217\n",
"Time: 17:07:12 Log-Likelihood: -16.841\n",
"No. Observations: 11 AIC: 37.68\n",
"Df Residuals: 9 BIC: 38.48\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0001 1.125 2.667 0.026 0.456 5.544\n",
"x1 0.5001 0.118 4.241 0.002 0.233 0.767\n",
"==============================================================================\n",
"Omnibus: 0.082 Durbin-Watson: 3.212\n",
"Prob(Omnibus): 0.960 Jarque-Bera (JB): 0.289\n",
"Skew: -0.122 Prob(JB): 0.865\n",
"Kurtosis: 2.244 Cond. No. 29.1\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.666\n",
"Model: OLS Adj. R-squared: 0.629\n",
"Method: Least Squares F-statistic: 17.97\n",
"Date: Thu, 13 Aug 2015 Prob (F-statistic): 0.00218\n",
"Time: 17:07:12 Log-Likelihood: -16.846\n",
"No. Observations: 11 AIC: 37.69\n",
"Df Residuals: 9 BIC: 38.49\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0009 1.125 2.667 0.026 0.455 5.547\n",
"x1 0.5000 0.118 4.239 0.002 0.233 0.767\n",
"==============================================================================\n",
"Omnibus: 1.594 Durbin-Watson: 2.188\n",
"Prob(Omnibus): 0.451 Jarque-Bera (JB): 1.108\n",
"Skew: -0.567 Prob(JB): 0.575\n",
"Kurtosis: 1.936 Cond. No. 29.1\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.666\n",
"Model: OLS Adj. R-squared: 0.629\n",
"Method: Least Squares F-statistic: 17.97\n",
"Date: Thu, 13 Aug 2015 Prob (F-statistic): 0.00218\n",
"Time: 17:07:12 Log-Likelihood: -16.838\n",
"No. Observations: 11 AIC: 37.68\n",
"Df Residuals: 9 BIC: 38.47\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0025 1.124 2.670 0.026 0.459 5.546\n",
"x1 0.4997 0.118 4.239 0.002 0.233 0.766\n",
"==============================================================================\n",
"Omnibus: 19.540 Durbin-Watson: 2.144\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 13.478\n",
"Skew: 2.041 Prob(JB): 0.00118\n",
"Kurtosis: 6.571 Cond. No. 29.1\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.667\n",
"Model: OLS Adj. R-squared: 0.630\n",
"Method: Least Squares F-statistic: 18.00\n",
"Date: Thu, 13 Aug 2015 Prob (F-statistic): 0.00216\n",
"Time: 17:07:12 Log-Likelihood: -16.833\n",
"No. Observations: 11 AIC: 37.67\n",
"Df Residuals: 9 BIC: 38.46\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"const 3.0017 1.124 2.671 0.026 0.459 5.544\n",
"x1 0.4999 0.118 4.243 0.002 0.233 0.766\n",
"==============================================================================\n",
"Omnibus: 0.555 Durbin-Watson: 1.662\n",
"Prob(Omnibus): 0.758 Jarque-Bera (JB): 0.524\n",
"Skew: 0.010 Prob(JB): 0.769\n",
"Kurtosis: 1.931 Cond. No. 29.1\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
},
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAlAAAAJPCAYAAACzYuoLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
"AAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+Q3XV97/HXS1JMQlIpUtGWTaGhqNelXBi0YMtl69VC\n",
"toq203rp2JbCzf5RnSR2equgU0nn3tpaS9tsps6dxoaiA1hLxavTEAm2q8yIiIo0ASlNWs1iNaAg\n",
"JDRRfrzvH98TObt7dnM+3/P9nu+P83zMMNnzPWfPeX/Ifj557ff7+X4+jggBAACgf8+rugAAAICm\n",
"IUABAAAkIkABAAAkIkABAAAkIkABAAAkIkABAAAkWjJA2d5u+4Dt3V3HftX2fbafsX1u+SUCQD6L\n",
"jGEn2d5l+0Hbt9k+scoaATTTsc5AXSfpknnHdkv6JUmfLaUiAChOrzHsKkm7IuJMSZ/uPAaAJEsG\n",
"qIi4Q9Jj8449EBEPlloVABSg1xgm6VJJ13e+vl7Sm4ZaFIBWYA4UgFFzSkQc6Hx9QNIpVRYDoJkI\n",
"UABGVmR7WbGfFYBky8p6Y9sMSsAIighXXcMxHLD94oj4lu2XSHp4/gsYv4DRlDJ+DRqglvygBgyk\n",
"fbG9OSI2V11HEdrSlra0Q2pdW5oQPD4h6XJJ7+v8+fFeL2L8qh/aUj9taYeUPn4daxmDmyR9TtJL\n",
"bc/avtL2m2zPSjpf0j/YvjV/uQBQnh5j2BWS/ljS62w/KOk1nccAkGTJM1AR8WuLPNXzNzYA/bPH\n",
"J6WxjdIZL7PXnS/NTkfs2VF1XW2yxBj22qEWAqB1SpsD1TIzVRdQoJmqCyjITNUFDCILTxdskbad\n",
"kTVl4iekqbX2uAhRKNhM1QUUaKbqAgo0U3UBBZmpuoCqOLsJpYQ3tqMtcwiAotnrdkq3Xrzwmcmd\n",
"ETvWDb+iYrSl37elHQD6l9rvWcYAqMTq5b2Pr1ox3DoAAHkQoIBKHDzS+/ihw8OtAwCQBwEKqMTs\n",
"tDS1d+6x9fuk/VurqQcAkII5UEBFsonkazZkl+0OHZb2b236BPK29Pu2tANA/1L7PQEKQGHa0u/b\n",
"0g4A/WMSOQAAQMkIUAAAAIkIUAAAAIkIUAAAAIkIUAAAAIkIUAAAAIkIUAAAAIkIUAAAAIkIUAAA\n",
"AIkIUAAAAIkIUAAAAIkIUAAAAIkIUAAAAImWVV0AAAAojj0+KY1tlFYvlw4ekWanI/bsqLqutiFA\n",
"AQDQEll4umCLtO2M545OrbXHRYgqFpfwAABojbGNc8OTlD1es6GaetqLAAUAQGusXt77+KoVw62j\n",
"/biEBwDAkJQ/P+ngkd7HDx0u7jOYZyURoAAAGIrhzE+anZam1s79jPX7pP1bi3l/5lkd5Ygo543t\n",
"iAiX8uYAaqkt/b4t7UC92Ot2SrdevPCZyZ0RO9YV9znjk9mcp1UrsjNP+7cWGWyG1Y5hS+33S56B\n",
"sr1d0i9KejgizuocO0nS30r6CUlfk/TmiPhu7ooBYMhsb5K0XpIlbYuILRWXhJEwnPlJnbBU4pkg\n",
"5llJx55Efp2kS+Ydu0rSrog4U9KnO48BoBFsjysLT6+UdLak19teW21VGA3DmZ9Uvra0YzBLBqiI\n",
"uEPSY/MOXyrp+s7X10t6Uwl1AUBZXibprog4EhHPSPqMpF+uuCaMhNlpaWrv3GPFzk8ajra0YzB5\n",
"JpGfEhEHOl8fkHRKgfUAaKDn7shphD2S/rAzHeGIsmkKX6i2JIyCiD07shOgk6XNTxqGtrRjUMec\n",
"RG77NEmf7JoD9VhE/EjX849GxEk9vo9JmMAImHtHjtWEfm/7SklvlfSkpPskfS8ifqfrecYvYMQU\n",
"Ool8EQdsvzgivmX7JZIeXqKYzV0PZyJiJsfnAagp2xPSmmnpx9dKm6sup28RsV3Sdkmy/V5J++e/\n",
"hvELaLds/NJE7u/PcQbqTyR9JyLeZ/sqSSdGxIKJ5PwGB4wG+80z0kcv6jxqyhmoF0XEw7bXSPqU\n",
"pJ+JiCe6nmf8AkZM0csY3CTpIkkn256V9B5Jfyzpo7b/pzrLGOQvF0DzLXZHTq3dbPuFkp6S9Nbu\n",
"8AQA/WAhTQADaeIcqGNh/AJGzzDmQAHAD8y9I2fBunEA0EqcgQJQmLb0+7a0A/1jc1xwBgoAgARs\n",
"jltfdQ62BCgAwIgb2zg3PEnZ48kNKnVPOSyl7sH2WHvhAQDQcmyOW0+LBds1G6qpZy4CFABgxLE5\n",
"bj3VO9gSoAAAI47Nceup3sGWOVAAgJHG5rh1NTstTa2dexmvPsGWZQwAFKYt/b4t7QCaLptIvmYo\n",
"wTa13xOgABSmLf2+Le0A0L/Ufs8cKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQE\n",
"KAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgETLqi4AAICl2OOT\n",
"0thGafVy6eARaXY6Ys+OquvCaCNAAQBqKwtPF2yRtp3x3NGptfa4CFGoEpfwAAA1NrZxbniSssdr\n",
"NlRTD9rGHp+01+1M/T7OQAEAamz18t7HV60Ybh1oo7lnOJ30vbnPQNneZHu37T22N+V9HwAYNttX\n",
"276vM4bdaPv5VdeExRw80vv4ocPDrQPt1OsMZ39yBSjb45LWS3qlpLMlvd722jzvBQDDZPs0SVOS\n",
"zo2IsyQdJ+myKmvCUmanpam9c4+t3yft31pNPWiXxc5wHlveS3gvk3RXRByRJNufkfTLkt6ftxAA\n",
"GJInJD0laaXtZyStlPSNakvCYiL27Mh+Z5/ckF22O3RY2r+VCeQoxmJnOI8tb4DaI+kPbZ8k6Yik\n",
"X5T0hbxFAMCwRMSjtq+VtF/SYUmfiojbKy4LS+iEJQITSjA7LU2tzXMZL1eAiogHbL9P0m2SnpR0\n",
"j6Rn57/O9uauhzMRMZPn8wDUk+0JSRMVl5GkM93g7ZJOk/S4pL+z/ZaIuGHe6zZ3PWT8AlqmM369\n",
"Svq3L2RxRkkhyhFRRBHvlbQ/Iv5v17GIiLQp7QAarQn93vb/kPS6iFjfefwbks6PiLd1vab27QBQ\n",
"rNR+P8hdeC/q/LlG0i9JujHvewHAED0g6XzbK2xb0msl3V9xTQAaZpB1oG62/UJlkzHfGhFPFFQT\n",
"AJQmIu61/SFJX1Q29eDLkv6q2qoANE0hl/B6vjGnwEcO+1WhLf2+Le0A0L/Ufs9K5CgE+1UBAEYJ\n",
"e+GhIOxXBQAYHQQoFIT9qgAAo4MAhYKwXxUAYHQQoFAQ9qsCAIwO7sJDYbKJ5GvYr2qEtaXft6Ud\n",
"APqX2u8JUAAK05Z+35Z2AOjf0FYiBwAAGFUEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQEKAAAgEQE\n",
"KAAAgEQEKAAAgEQEKAAAgEQEKAAAgETLqi4AqKNsX7+xjdLq5dLBI9LsNPv6AQCOIkAB82Th6YIt\n",
"0rYznjs6tdYeFyEKACCxmTCwgL1up3TrxQufmdwZsWPd8Ctqjrb0+7a0Yxg4W4u2SO33nIECFli9\n",
"vPfxVSuGWwdQb5ytxShjEjmwwMEjvY8fOjzcOoC6G9s4NzxJ2eM1G6qpBxgeAhSwwOy0NLV37rH1\n",
"+6T9W6upB6grztZidHEJD5gnYs8Oe1zS5IbsH4JDh6X9W7kkAczH2VqMLgIU0EMnLBGYgCXNTktT\n",
"a+dexuNsLUYDd+EBKExb+n1b2jEM2UTyNZytReOl9nsCFIDCtKXft6UdAPo3tGUMbF8t6dclPStp\n",
"t6QrIuJ7ed8P6AdrzqAItl8q6SNdh35S0u9HxHRFJQFomFwByvZpkqYkvTwivmf7byVdJun64koD\n",
"5mLNGRQlIv5F0jmSZPt5kr4h6ZZKiwLQKHmXMXhC0lOSVtpeJmmlsgEIKBFrzqAUr5W0LyJmqy4E\n",
"QHPkClAR8aikayXtl/Qfkr4bEbcXWRiwEGvOoBSXSbqx6iIANEuuAGV7raS3SzpN0o9JWmX7LQXW\n",
"BfTAmjMolu3jJb1B0t9VXQuAZsk7ifw8SZ+LiO9Iku2PSXq1pBu6X2R7c9fDmYiYyfl5gFhzpn5s\n",
"T0iaqLiMQayT9KWIeGT+E4xfQLsNOn7lWsbA9tnKwtIrJR2R9DeSvhARf9n1Gm4DRuFYc6bemtbv\n",
"bX9E0q0Rcf28441qB4DBDW0dKNvvkHS5smUMvixpfUQ8lbcQAM3XpH5v+wRJX5d0ekQcnPdcY9oB\n",
"oBgspAmgMm3p921pB4D+pfb7vMsYAAAAjCwCFAAAQCICFAAAQCICFAAAQCICFAAAQCICFAAAQCIC\n",
"FAAAQCICFAAAQKK8e+EBaIBs65uxjdLq5dlmzLPTbH0DAIMjQAEtlYWnC7bM3Xx5aq09LkIUAAyG\n",
"S3hAa41tnBuepOzxmg3V1AMA7UGAAlpr9fLex1etGG4dANA+BCigtQ4e6X380OHh1gEA7UOAAlpr\n",
"dlqa2jv32Pp90v6t1dQDAO3hiCjnje2ICJfy5gD6kk0kX7Mhu2x36LC0f2uZE8jb0u/b0g4A/Uvt\n",
"9wQoAIVpS79vSzsA9C+133MJDwAAIBEBCgAAIBEBCgAAIBEBCgAAIBFbuSyBfcQAAEAvBKhFsI8Y\n",
"AABYDJfwFsU+YgAAoDcC1KLYRwwAAPTGJbxFsY8YgOZiDidQLgLUomanpam1cy/jsY8YgPpjDidQ\n",
"PrZyWcKw9xEDmq4N/V5qfjvsdTulWy9e+Mzkzogd64ZfEVB/qf0+1xko2y+V9JGuQz8p6fcjYjrP\n",
"+9VVJywRmICWsX2ipA9KeoWkkHRlRHy+2qqKxBxOoGy5AlRE/IukcyTJ9vMkfUPSLQXWBQBl2iJp\n",
"R0T8iu1lkk6ouqBiMYcTKFsRd+G9VtK+iJgt4L0AoFS2XyDpwojYLkkR8XREPF5xWQWbnZam9s49\n",
"xhxOoEhFTCK/TNKNBbwPAAzD6ZIesX2dpLMlfUnSpoj4z2rLKk7Enh32uKRJ5nACJRloErnt45Vd\n",
"vvsvEfHIvOcaPQmzbbilGcPQhH5v+zxJd0p6dUTcbfsvJD0REe/pek3t2wGgWEOZRN5lnaQvzQ9P\n",
"XcVs7no4ExEzA34ecuCWZpTF9oSkiYrLSPWQpIci4u7O45slXTX/RYxfQLsNOn4NegbqI5JujYjr\n",
"ezzHb3A1wS3NGJam9Hvbn5W0PiIe7ASlFRHxzq7nG9EOAMUZ2hko2ycom0A+lfc9MCzc0gzMs0HS\n",
"DZ1pCPskXVFxPQAaJneAiognJZ1cYC0oDbc0A90i4l5Jr6y6DgDNxWbCI4FbmgEAKBJbuYwItqXB\n",
"MLSl37elHQD6l9rvCVAACtOWft+WdgDoX2q/5xIeAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABA\n",
"IgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIU\n",
"AABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABA\n",
"omV5v9H2iZI+KOkVkkLSlRHx+aIKA4Ay2f6apCckPSPpqYh4VbUVAWiS3AFK0hZJOyLiV2wvk3RC\n",
"QTUBwDCEpImIeLTqQgA0T64AZfsFki6MiMslKSKelvR4kYUBwBC46gIANFPeOVCnS3rE9nW2v2x7\n",
"m+2VRRYGACULSbfb/qLtqaqLAdAseQPUMknnSvpARJwr6UlJVxVWFQCU72cj4hxJ6yS9zfaFVRcE\n",
"oDnyzoF6SNJDEXF35/HN6hGgbG/uejgTETM5Pw9ADdmekDRRcRm5RMQ3O38+YvsWSa+SdMfR5xm/\n",
"gHYbdPxyROT94M9KWh8RD3YGmhUR8c6u5yMimF8AjJCm9PvOlIPjIuKg7RMk3SbpDyLits7zjWgH\n",
"gOKk9vtB7sLbIOkG28dL2ifpigHeCwCG6RRJt9iWsnHwhqPhCQD6kfsM1DHfmN/ggJHTln7flnYA\n",
"6F9qv2clcgAAgEQEKAAAgESDzIECAORgj09KYxul1culg0ek2emIPTuqrgtA/whQADBEWXi6YIu0\n",
"7Yznjk6ttcdFiAKag0t4ADBUYxvnhicpe7xmQzX1AMijkWegOP0NoLlWL+99fNWK4dYBYBCNC1Cc\n",
"/gbQbAeP9D5+6PBw6wAwiAZewuP0N4Amm52WpvbOPbZ+n7R/azX1AMijcWegOP0NoMki9uywxyVN\n",
"bsjGrUOHpf1bOYMONEsDAxSnvwE0WycsEZiABmvgJTxOfwMAgGo1ci+8bCL5Gk5/AzXTlj3k2tIO\n",
"AP1L7feNDFAA6qkt/b4t7QDQPzYTBgAAKBkBCgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEB\n",
"CgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEBCgAAIBEBCgAAINGyvN9o\n",
"+2uSnpD0jKSnIuJVRRUFAGWzfZykL0p6KCLeUHU9AJold4CSFJImIuLRoooBgCHaJOl+SaurLgRA\n",
"8wx6Cc+FVAEAQ2T7VEmTkj4oxjEAOQwSoELS7ba/aHuqqIIAYAj+XNLvSXq26kIANNMgl/B+NiK+\n",
"aftHJe2y/UBE3NH9Atubux7ORMTMAJ8HoGZsT0iaqLiMJLZfL+nhiLinU/9ir9vc9ZDxC2iZQccv\n",
"R0QRRVwj6VBEXNt1LCKCU+PACGlCv7f9Xkm/IelpScsl/bCkv4+I3+x6Te3bAaBYqf0+1yU82ytt\n",
"r+58fYKkX5C0O897AcAwRcS7ImIsIk6XdJmkf+wOTwDQj7yX8E6RdIvto+9xQ0TcVlhVADA8g5+G\n",
"BzByCrmE1/ONOQXeF3t8UhrbKK1eLh08Is1OR+zZUXVdQB5t6fdtaQeA/qX2+0EmkWNAWXi6YIu0\n",
"7Yznjk6ttcdFiAIAoL7YyqVSYxvnhicpe7xmQzX1AACAfhCgKrV6ee/jq1YMtw4AAJCCAFWpg0d6\n",
"Hz90eLh1AACAFASoSs1OS1N75x5bv0/av7WaegAAQD+4C69i2UTyNRuyy3aHDkv7tzKBHE3Vln7f\n",
"lnYA6F9qvydAAShMW/p9W9oBoH9DWYkcAABglBGgAAAAEhGgAAAAEhGgAAAAEhGgAAAAEhGgAAAA\n",
"EhGgAAAAEhGgAAAAEhGgAAAAEhGgAAAAEhGgAAAAEhGgAAAAEpUaoOx1O+3xyTI/AwAAYNhKPgN1\n",
"68XSBVsIUQAAoE2GcAlv2xnSmg3lfw4AAMBwDGkO1KoVw/kcAACA8g0pQB06PJzPAQAAKN8QAtT6\n",
"fdL+reV/DgAAwHCUHKAmd0qf3xixZ0e5nwMAxeIuYgBLcUSU88Z2RIRLeXMAtdSWfm87pJA0tVe6\n",
"cxO/BALtlzp+DXQGyvZxtu+x/clB3qfubE9UXUNR2tKWtrRDaldbmsL2ctt32f6K7ftt/1HvVzb/\n",
"LuI2/XzRlvppSzvyGPQS3iZJ9yv7Va3NJqouoEATVRdQkImqCyjQRNUFjJqIOCLp5yPiv0r6aUk/\n",
"b/vner+68XcRT1RdQIEmqi6gQBNVF1CQiaoLqEruAGX7VEmTkj4oqfGn7AGMloj4z86Xx0s6TtKj\n",
"vV/JXcQAFhrkDNSfS/o9Sc8WVAsADI3t59n+iqQDkv4pIu5f+CruIgbQW65J5LZfL2ldRLytc/3z\n",
"dyPiDfNe0/bLegB6aNokctsvkPQpSVdFxEznGOMXMIJSxq9lOT/j1ZIutT0pabmkH7b9oYj4zTxF\n",
"AEBVIuJx2/8g6TxJM51jjF8AljTwMga2L5L0v+afgQKAurJ9sqSnI+K7tlcoOwP1BxHx6YpLA9AQ\n",
"ec9AzcfpbgBN8hJJ19t+nrK5oB8mPAFIUdpCmgAAAG1VylYubVlg0/aJtm+2/dXOYnvnV11TXrav\n",
"tn2f7d22b7T9/Kpr6pft7bYP2N7ddewk27tsP2j7NtsnVlljvxZpy/s7P2P32v5YZ1JzrfVqR9dz\n",
"v2v7WdsnVVFbERjD6oXxqx7aMn5JxYxhZe2F15YFNrdI2hERL1e22N5XK64nF9unSZqSdG5EnKVs\n",
"zZvLqqwp0XWSLpl37CpJuyLiTEmf7jxugl5tuU3SKyLibEkPSrp66FWl69UO2R6T9DpJXx96RcVi\n",
"DKsJxq9aacv4JRUwhhUeoNqywGYnRV8YEdslKSKejojHKy4rryckPSVppe1lklZK+ka1JfUvIu6Q\n",
"9Ni8w5dKur7z9fWS3jTUonLq1ZaI2BURR9dTu0vSqUMvLNEifyeS9GeS3jHkcgrFGFY7jF810Zbx\n",
"SypmDCvjDFRbFtg8XdIjtq+z/WXb22yvrLqoPCLiUUnXStov6T8kfTcibq+2qoGdEhEHOl8fkHRK\n",
"lcUU6EpJjdy41vYbJT0UEf9cdS0DYgyrEcavRmns+CWlj2GFBqjOApsPR8Q9avBvbh3LJJ0r6QMR\n",
"ca6kJ9Wc06xz2F4r6e2STpP0Y5JW2X5LpUUVKLI7IZp+qUW23y3p+xFxY9W1pOr8w/wuSdd0H66o\n",
"nNwYw+qH8asZmjx+SfnGsKLPQB1dYPPfJd0k6TW2P1TwZwzLQ8qS6N2dxzcrG4ya6DxJn4uI70TE\n",
"05I+puzvqskO2H6xJNl+iaSHK65nILZ/S9llo6b+w7BW2T9w93b6/6mSvmT7RZVWlY4xrH4Yv2qu\n",
"BeOXlGMMKzRARcS7ImIsIk5XNsnvH7tXJ2+SiPiWpFnbZ3YOvVbSfRWWNIgHJJ1ve4VtK2tLj32/\n",
"GuUTki7vfH25pI9XWMtAbF+i7JLRGyPiSNX15BERuyPilIg4vdP/H1I26bdR/zAwhtUS41eNtWH8\n",
"kvKNYWXdhfeDmkp+/7JtkHSD7XuV3cHy3orrySUi7pX0IUlflHT02u5fVVdRGts3SfqcpJfanrV9\n",
"haQ/lvQ62w9Kek3nce31aMuVkrZKWiVpV+fW+Q9UWmQfutpxZtffSbem9/2jmt6Oxo9hjF/10Zbx\n",
"SypmDGMhTQAAgERln4ECAABoHQIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIU\n",
"AABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABA\n",
"IgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIU\n",
"AABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAoiUDlO3ttg/Y3t117H/bvtf2\n",
"V2x/2vZY+WUCQDrGMABlcUQs/qR9oaRDkj4UEWd1jq2OiIOdrzdIOjsi1g+jWABIwRgGoCxLnoGK\n",
"iDskPTbv2MGuh6skfbuEugBgYIxhAMqyLM832f5DSb8h6T8lnV9oRQBQMsYwAINa8hKeJNk+TdIn\n",
"j57+nvfcVZJeGhFX9Hhu6TcG0EoR4apr6JZnDGP8AkZTyviV6wxUlxsl7SiikDqzvTkiNlddRxHa\n",
"0pa2tENqXVuaFjwWHcOaMn415eeHOovXlFobVGfS+JW8jIHtn+p6+EZJ96S+BwBUhTEMQBGWPANl\n",
"+yZJF0k62faspGskTdp+qaRnJO2T9NulVwkAOTCGASjLkgEqIn6tx+HtJdVSZzNVF1CgmaoLKMhM\n",
"1QUUaKbqAtpqRMawmaoL6NNM1QX0aabqAhLMVF1An2aqLqAMx5xEnvuN7WjKHAKgrezxSWlso7R6\n",
"uXTwiDQ7HbFn0XmLg39eO/p9W9oBoH+p/X7QSeQAaioLTxdskbad8dzRqbX2uMoMUQAwCtgLD2it\n",
"sY1zw5OUPV6zoZp6AKA9CFBAa61e3vv4qhXDrQMA2ocABbTWwSO9jx86PNw6AKB9CFBAa81OS1N7\n",
"5x5bv0/av7WaegCgPbgLD2ixbCL5mg3ZZbtDh6X9W7kL79ja0g4A/Uvt9wQoAIVpS79vSzsA9I9l\n",
"DAAAQCsNe227pRCgAABA7dVtbTsmkQMAgAao19p2BCgAANAA9VrbjgAFAAAaoF5r2xGgAABAA9Rr\n",
"bTuWMQBQmLb0+7a0A2ibMte2Yx0oAJVpS79vSzsA9C+133MJDwAAIBEBCgAAIBEBCgAAIBEBCgAA\n",
"IBEBCgAAIBEBCgAAINGSAcr2dtsHbO/uOvZ+21+1fa/tj9l+QfllAkA6xjAAZTnWGajrJF0y79ht\n",
"kl4REWdLelDS1WUUBgAFYAwDUIolA1RE3CHpsXnHdkXEs52Hd0k6taTaAGAgjGEAyjLoHKgrJRWy\n",
"hDoAVIAxDEAuuQOU7XdL+n5E3FhgPQAwFIxhAAaxLM832f4tSZOS/vsxXre56+FMRMzk+TwA9WR7\n",
"QtJExWUk62cMY/wC2m3Q8euYmwnbPk3SJyPirM7jSyRdK+miiPj2Et/HZpzAiKljv88zhtWxHQDK\n",
"ldrvlwxQtm+SdJGkkyUdkHSNsjtWjpf0aOdld0bEWwctBEDz1a3f5x3D6tYOAOUrNEANsxAAzdeW\n",
"ft+WdgDoX2q/ZyVyAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACA\n",
"RAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACARAQo\n",
"AACARAQoAACARAQoAACARAQoAACARAQoAACARAQoAACAREsGKNvbbR+wvbvr2K/avs/2M7bPLb9E\n",
"AMiHMQxAWY51Buo6SZfMO7Zb0i9J+mwpFQFAcRjDAJRi2VJPRsQdtk+bd+wBSbJdXlUAUADGMABl\n",
"YQ4UAABAIgIUAABAoiUv4Q3K9uauhzMRMVPm5wEYLtsTkiYqLqMUjF9Auw06fg0aoJacRBARmwd8\n",
"fwA11gkVM0cf276msmLyWXQMY/wC2m3Q8csRsfiT9k2SLpJ0sqQDkq6R9KikrZ1jj0u6JyLW9fje\n",
"iAhmaQIjpG79Pu8YVrd2AChfar9fMkANsxAAzdeWft+WdgDoX2q/ZxI5AABAIgIUAABAIgIUAABA\n",
"IgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIU\n",
"AABAIgKHaqkQAAARy0lEQVQUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABAIgIUAABA\n",
"IgIUAABAIgIUAABAIgIUAABAoiUDlO3ttg/Y3t117CTbu2w/aPs22yeWXyYApGMMA1CWY52Buk7S\n",
"JfOOXSVpV0ScKenTnccAUEeMYQBKsWSAiog7JD027/Clkq7vfH29pDeVUBcADIwxDEBZ8syBOiUi\n",
"DnS+PiDplALrAYCyMYYBGNiyQb45IsJ2LPa87c1dD2ciYmaQzwNQL7YnJE1UXEZuS41hjF9Auw06\n",
"fjli0fxz9ANOk/TJiDir8/gBSRMR8S3bL5H0TxHxsh7fFxHhvIUBaJ469vs8Y1gd2wGgXKn9Ps8l\n",
"vE9Iurzz9eWSPp7jPQCgKoxhAAa25Bko2zdJukjSycrmCrxH0v+T9FFJayR9TdKbI+K7Pb6X3+CA\n",
"EVO3fp93DKtbOwCUL7XfH/MS3rAKAdB8ben3bWkHgP4N4xIeAADASCNAAQAAJCJAAQAAJCJAAQAA\n",
"JCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJA\n",
"AQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJCJAAQAAJModoGxv\n",
"sr3b9h7bm4osCgDKxPgFYFC5ApTtcUnrJb1S0tmSXm97bZGFAUAZGL8AFCHvGaiXSborIo5ExDOS\n",
"PiPpl4srCwBKw/gFYGB5A9QeSRfaPsn2Skm/KOnU4soCgNIwfgEY2LI83xQRD9h+n6TbJD0p6R5J\n",
"zxZZGACUgfELQBFyBShJiojtkrZLku33Sto//zW2N3c9nImImbyfB6B+bE9Imqi4jGSMXwAGHb8c\n",
"EXk/+EUR8bDtNZI+JelnIuKJrucjIpy3MADN05R+z/gFYL7Ufp/7DJSkm22/UNJTkt7aPfgAGC32\n",
"+KQ0trHqOhIcc/yy1+2UZqcj9uwYfnkA6i73GahjvjG/wQEjIQtPF2yRtp0hWW3o97ZDCklTe6U7\n",
"NxGigPZLzS2sRA5gQGMbs/DURtvOkNZsqLoKAPVDgAIwoNXLq66gXKtWVF0BgPohQAEY0MEjVVdQ\n",
"rkOHq64AQP0QoAAMaHY6myvURuv3Sfu3Vl0FgPoZ5C48AFDEnh3Z9nKTGyRdUnU9xZncKe3fygRy\n",
"AL1wFx6AwrSl37elHQD6x114AAAAJSNAAQAAJCJAAQAAJCJAAQAAJCJAAUAP9rqd2TY1ALAQAQoA\n",
"err1YumCLYQoAL0QoABgUeyFB6A3AhQALIm98AAsRIACgCWxFx6AhQhQALAo9sID0Bt74QFAT+yF\n",
"B2Bx7IUHoDBt6fdtaQeA/rEXHgAAQMkIUAAAAImYAwVUJFugcWyjtHq5dPCINDvNfBsAaAYCFFCB\n",
"LDxdsCVbqPGoqbX2uAhRAFB/XMIDKjG2cW54klj1GgCaI3eAsn217fts77Z9o+3nF1kY0G6rl/c+\n",
"zqrXw8D4BWBQuQKU7dMkTUk6NyLOknScpMuKKwtou4NHeh9n1euyMX4BKELeM1BPSHpK0krbyySt\n",
"lPSNwqoCWm92WpraO/cYq14PCeMXgIHlmkQeEY/avlbSfkmHJX0qIm4vtDKgxSL27LDHJU1uyC7b\n",
"HTrMqtfDwfgFoAi5ViK3vVbSJyVdKOlxSX8n6eaIuKHrNazkC4yYJvR7xi8AvaT2+7zLGJwn6XMR\n",
"8Z3Oh35M0qsl3dD9Itubux7ORMRMzs8DUEO2JyRNVFxGqj7Hr5/YKz1yl3T4X8X4BbTOoONX3jNQ\n",
"ZysbbF4p6Yikv5H0hYj4y67X8BscMGKa0O/7Hb+kUDZP7c5NXFoF2m8oe+FFxL2SPiTpi5L+uXP4\n",
"r/K8FwAMU9r4xdpcAHrLdQaqrzduwG+iwGLYZiWftvT7585ASdKbPxPx0Ykq6wFQvmHNgQJai21W\n",
"MBdrcwFYiK1cgAXYZgVHsTYXgN44AwUswDYrkKR135ZmP8xZRwC9cAYKWIBtViBJt54sXfDr2SVd\n",
"AJiLAAUswDYrOIpLtwB64xIeMA/brGAuLt0CWIgABfTQCUsEJohLtwB64RIeACyKS7cAeuMMFAD0\n",
"NLmTS7cAFsNK5AAK05Z+35Z2AOgfK5Gj1dhiBQBQBwQoNAZbrAAA6oJJ5GgQtlgBANQDAQoNwhYr\n",
"AIB6IEChQdhiBQBQDwQoNAhbrAAA6oFlDNAo2UTyNWyxUlNt6fdtaQeA/qX2ewIUgMK0pd+3pR0A\n",
"+pfa77mEBwAAkIgABQAAkIgABQAAkIiVyFEYtlkBAIwKAhQKwTYrAIBRkusSnu2X2r6n67/HbW8s\n",
"ujg0CdusoDn6GcPsdTuzXwwAYKFcZ6Ai4l8knSNJtp8n6RuSbimwLjQO26ygOfobw269mLOoABZT\n",
"xCTy10raFxGzBbwXGottVtBYS4xhnEUF0FsRAeoySTcW8D5oNLZZQWMdYwzjLCqAhQaaRG77eElv\n",
"kPTORZ7f3PVwJiJmBvk81FfEnh32uKRJtlkZIbYnJE1UXEZuS49hmzt/3vXjticYv4B2GXT8Gmgr\n",
"F9tvlPTbEXFJj+fYCgEYMU3r94uNYbZDerek+74p/et6fhEA2m/YW7n8mqSbBnwPAKjKEmPY/5H0\n",
"I8OsBUCD5D4DZfsESV+XdHpEHOzxfKN+EwUwuCb1+6XGsOwM1NGxcXJnxI51Qy8QwFCljl+550BF\n",
"xJOSTs77/QBQpf7HMCaRA1iIlchHBNusAHmxFAeAhQhQI4BtVoC8WIoDQG9FrAOF2mObFSDd5E7p\n",
"8xv5JQNAL5yBGglsswKkYuI4gKVwBmoksM0KAABFIkCNBLZZAQCgSAOtRL7kGzdoPZhRkE0kX8M2\n",
"KyhVW/p9tg7UJZ/iblVgdKSOXwQoAIVpS79/biHNqb3SnZsIUUD7DXsrFwBoMe5WBdAbAQoAlsTd\n",
"qgAWIkABwJK4WxXAQqwDVTG2WAHqjLtVAfRGgKoQW6wAdTa5k7tVASyGu/AqZK/bKd168cJnJney\n",
"CjKaqC39vi3tANA/7sJrFLZYAQCgiQhQlWKLFQAAmogAVSm2WAEAoImYA1UxtlhBm7Sl37elHQD6\n",
"x1YuACrTln7flnYA6B+TyAEAAEpGgAIAAEhEgAIAAEiUO0DZPtH2zba/avt+2+cXWVgd2OOT9rqd\n",
"9ptnsj/HJ6uuCcDgRmH8AlCuQbZy2SJpR0T8iu1lkk4oqKZaYJsVoNVaPX4BKF+uu/Bsv0DSPRHx\n",
"k0u8ptF3sbDNCpCuCf1+FMYvAOmGdRfe6ZIesX2d7S/b3mZ7Zc73qim2WQFaagTGLwBlyxuglkk6\n",
"V9IHIuJcSU9KuqqwqmqBbVaAlhqB8QtA2fLOgXpI0kMRcXfn8c3qMQDZ3tz1cCYiZnJ+XgVmp6Wp\n",
"tXPnQLHNCtDN9oSkiYrLSDUC4xeAYxl0/Mq9Erntz0paHxEPdgaaFRHxzq7nGz+HgG1WgDRN6fej\n",
"MH4BSDO0rVxsny3pg5KOl7RP0hUR8XjeQgA0X1P6PeMXgPnYCw9AZdrS79vSDgD9Yy88AACAkhGg\n",
"AKAHdh8AsJRBViKvTDaojW3M1mo6eESanWZyN4Bi3Xoxuw8AWEzjAhRbrAAYnm1nSJMbJDG2AJij\n",
"gZfwxjbODU9S9njNhmrqAdBu7D4AYKEGBii2WAEwTOw+AGChBgYotlgBMCzsPgCgtwYGqNlpaWrv\n",
"3GMMcgCKtu7b0uc/zNxKAL00ciFNtlgB6qktC1DaDimU/bJ25ybGF6D9WIkcQGXa0u+fC1CSNLkz\n",
"Yse6SgsCUDpWIgeAQnGDCoCFCFAAsCRuUAGwUKkBiq0QADQbN6gA6K3klcjZCgFAU03u5AYVAIsp\n",
"dRI5kzCB0dKmSeRtaAeA/tV0EjmTMAEAQHsMKUAxCRMAALTHEAIUkzABNA83wQBYSsmTyJmECaCp\n",
"uAkGwOJYiRxAYdrS77kJBhg9NZ1EDgBNxU0wABYiQAHAkrgJBsBCBCgAWBQ3wQDoLfckcttfk/SE\n",
"pGckPRURryqqKAAo27HHMG6CAbC4Qc5AhaSJiDin7eHJ9kTVNRSlLW1pSzukdrWlYZYcwyJ2rGtC\n",
"eGrKzw91Fq8ptTalzlSDXsJr/N02fZqouoACTVRdQEEmqi6gQBNVFzDC2jCGTVRdQJ8mqi6gTxNV\n",
"F5BgouoC+jRRdQFlGPQM1O22v2h7qqiCAGBIGMMA5DbIQpo/GxHftP2jknbZfiAi7iiqMAAoGWMY\n",
"gNwKWUjT9jWSDkXEtV3HylmhE0CtNXEhzfljGOMXMJpSxq9cZ6Bsr5R0XEQctH2CpF+Q9Ad5iwCA\n",
"YTrWGMb4BeBY8l7CO0XSLbaPvscNEXFbYVUBQLkYwwAMpLS98AAAANqqlJXIbR9n+x7bnyzj/YfF\n",
"9om2b7b9Vdv32z6/6prysn217fts77Z9o+3nV11Tv2xvt33A9u6uYyfZ3mX7Qdu32T6xyhr7tUhb\n",
"3t/5GbvX9sdsv6DKGvvRqx1dz/2u7Wdtn1RFbYOqa19pUj9oys95U36OF6vT9obO/9M9tt9XVX1d\n",
"9fT6e3+V7S90MsHdtl9ZZY2dmsZs/1Onn++xvbFzPKk/lbWVyyZJ9+u57cybaoukHRHxckk/Lemr\n",
"FdeTi+3TJE1JOjcizpJ0nKTLqqwp0XWSLpl37CpJuyLiTEmf7jxugl5tuU3SKyLibEkPSrp66FWl\n",
"69UO2R6T9DpJXx96RQWoeV9pUj9oys95U36OF9Rp++clXSrppyNiXNKfVlHYPL3+f/6JpN+PiHMk\n",
"vafzuGpPSfqdiHiFpPMlvc32y5XYnwoPULZPlTQp6YNq8CJ1nd+OLoyI7ZIUEU9HxOMVl5XXE8p+\n",
"YFbaXiZppaRvVFtS/zq3lj827/Clkq7vfH29pDcNtaicerUlInZFxLOdh3dJOnXohSVa5O9Ekv5M\n",
"0juGXE6RattXmtQPmvJz3pSf40Xq/G1JfxQRT3Ve88jQC5tnkTq/Keno2cYTVYP+FBHfioivdL4+\n",
"pOzkyI8rsT+VcQbqzyX9nqRnj/XCmjtd0iO2r7P9ZdvbOnfuNE5EPCrpWkn7Jf2HpO9GxO3VVjWw\n",
"UyLiQOfrA8omBbfBlZJqv31IL7bfKOmhiPjnqmvJq4F9pan9oLY/5w36Of4pSf/N9udtz9g+r+qC\n",
"FnGVpGtt75f0ftXjzOMPdM46n6Ms1Cf1p0IDlO3XS3o4Iu5Rg88+dSyTdK6kD0TEuZKeVH1Ojyex\n",
"vVbS2yWdJunHJK2y/ZZKiypQZHdCNP1ysWy/W9L3I+LGqmtJ1fnl4l2Sruk+XFE5uTW5rzSlH9T5\n",
"57xhP8fLJP1IRJyv7KTFRyuuZzF/LWljRKyR9DuStldczw/YXiXp7yVtioiD3c/105+KPgP1akmX\n",
"2v53STdJeo3tDxX8GcPykLLfQu7uPL5ZWaBqovMkfS4ivhMRT0v6mLK/qyY7YPvFkmT7JZIerrie\n",
"gdj+LWWXvhvxj3UPa5WFjns7/f9USV+y/aJKq0rXtL7SqH7QgJ/zJv0cP6Ts51Odf6eetf3Cakvq\n",
"6VURcUvn65slLdi4uwq2f0hZePpwRHy8czipPxUaoCLiXRExFhGnK5t4+Y8R8ZtFfsawRMS3JM3a\n",
"PrNz6LWS7quwpEE8IOl82ytsW1lb7q+4pkF9QtLlna8vl/TxJV5ba7YvUfYb5Bsj4kjV9eQREbsj\n",
"4pSIOL3T/x9SNhG71v+g99C0vtKYftCEn/OG/Rx/XNJrJKnz79TxEfGdakvqaa/tizpfv0bZDQSV\n",
"6vTtv5Z0f0T8RddTaf0pIkr5T9JFkj5R1vsP4z9JZ0u6W9K9ypL+C6quaYC2vENZANytbHLcD1Vd\n",
"U0LtNymbj/J9SbOSrpB0kqTblXXG2ySdWHWdOdtypaR/VXa3zz2d/z5QdZ0J7fje0b+Tec//m6ST\n",
"qq4zZ9tq2Vea1A+a8nPelJ/jXnVK+iFJH+78nH5J0kSN6uz+GT1P2fyir0i6U9I5Najz55TN0/5K\n",
"18/jJan9iYU0AQAAEpW1DhQAAEBrEaAAAAASEaAAAAASEaAAAAASEaAAAAASEaAAAAASEaAAAAAS\n",
"EaAAAAAS/X8AdFeS7qn6BQAAAABJRU5ErkJggg==\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize = (10,10))\n",
"i = 1\n",
"for X, Y in zip(x, y):\n",
" result = sm.OLS(Y, sm.add_constant(X)).fit()\n",
" print(result.summary())\n",
" plt.subplot(2, 2, i)\n",
" plt.plot(X, Y, 'o')\n",
" i += 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El primer paso, antes de liarnos la manta a la cabeza y ponernos a calcular como locos debería ser un análisis exploratorio inicial de los datos incluyendo visualización y 'limpieza' en el caso de que fuera necesario.\n",
"\n",
"Los gráficos de la derecha parece evidente que no tienen un comportamiento lineal. Otros modelos serían más adecuados. En el gráfico de abajo a la izquierda parece evidente que tenemos un *outlier* que habría que filtrar."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Apuntes finales"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hemos hecho un breve repaso a través de uno de los modelos más simples que podemos usar para regresión (aprendizaje supervisado). Pero como siempre, solo tenemos una parte de la película y a partir de esa parte, a veces, es difícil conocer el final.\n",
"\n",
"Cuando aplicamos este modelo podemos seguir una serie de pasos.\n",
"\n",
"* Visualizar los datos.\n",
"* Visualizar los datos.\n",
"* Visualizar los datos.\n",
"* Visualizar los datos.\n",
"* Filtrar si es necesario.\n",
"* Considerar si existe una relación lineal entre las entradas y la salida y es razonable aplicar este modelo.\n",
"* Visualizar los residuos (homocedasticidad/heterocedasticidad, lineal/no lineal,...).\n",
"* Realizar una serie de tests que nos permitan no descartar la presunción de que el modelo es lineal (otros además de los que hemos hecho, statsmodels al rescate) y aplicable a nuestro problema de forma satisfactoria.\n",
"* Si la cantidad de información lo permite realizar validaciones cruzadas y comparar con otros modelos posibles.\n",
"* No enamorarte del modelo."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Referencias"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* [An Introduction to Statistical Learning](http://www-bcf.usc.edu/~gareth/ISL/ISLR%20Fourth%20Printing.pdf). Gran parte de toda la explicación sigue parte del capítulo 3 de este libro. En algunos puntos he intentado extender un poco la información del libro para hacerlo más accesible.\n",
"\n",
"* [Entrada en el blog de Connor Johnson](http://connor-johnson.com/2014/02/18/linear-regression-with-python/). Una referencia completa a regresión lineal usando statsmodels explicando el significado de muchas cosas interesantes que ofrece la salida y que he obviado en este test por simplificarlo un poco.\n",
"\n",
"* [Documentación de statsmodels](http://www.statsmodels.org/stable/index.html)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Post-Data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Como siempre, si encontráis errores en el texto no dudéis en indicarlos:\n",
"\n",
"* Haced un fork de https://github.com/Pybonacci/notebooks, modificad lo que consideréis y mandad un *pull request*.\n",
"* Comentad en la entrada del blog que corresponde a este notebook. http://pybonacci.org/2015/08/13/regresion-line…mple-explicada/\n",
"\n",
"Sería interesante que os descargáseis el notebook y lo ejecutáseis ya que hay algunas animaciones interactivas que os pueden ayudar a entender mejor algún concepto."
]
}
],
"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.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}