{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# SME0822 Análise Multivariada e Aprendizado não-Supervisionado\n", "\n", "Por Cibele Russo - ICMC USP\n", "\n", "## Aula 5b: Aplicações em Python de testes de hipóteses para o vetor de médias\n", "\n", "- Testes de hipóteses para a média populacional (multidimensional)\n", "- Testes para a comparação de médias populacionais em amostras independentes\n", "- Testes para a comparação de médias populacionais em amostras correlacionadas\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testes de hipóteses para a média (multidimensional)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Seja $\\underline{X}_1,\\ldots,\\underline{X}_n$ uma amostra aleatória de uma distribuição normal p-variada com vetor de médias $\\underline{\\mu}$ e matriz de variâncias e covariâncias $\\Sigma$. Sejam $\\overline{\\underline{X}}$ e $S$ o vetor de médias amostrais e a matriz de variâncias e covariâncias amostrais.\n", "\t\t\t\n", "Queremos avaliar se\n", "\t\t\t\n", "$$\\begin{array}{l}H_0:\\underline{\\mu}=\\underline{\\mu}_0\\mbox{ contra }\\\\H_1:\\underline{\\mu}\\neq\\underline{\\mu}_0.\\end{array}$$ \n", "\t\t\n", "\n", "Temos que\n", "\n", "\n", "1. $\\overline{\\underline{X}}\\sim N_p\\left(\\underline{\\mu},\\displaystyle\\frac{\\Sigma}{n}\\right)$.\n", "2. $(n-1)S\\sim Wishart(n-1)$.\n", "3. $\\overline{\\underline{X}}$ e $S$ são independentes. \n", "\n", "\n", "\n", "Sob $H_0$, a estatística $T^2$ de Hotelling\n", "$$T^2 = \\sqrt{n}(\\overline{\\underline{X}}-\\underline{\\mu}_0)^\\top \\displaystyle\\left(\\frac{(n-1)S}{n-1}^{-1}\\right)\\sqrt{n}(\\overline{\\underline{X}}-\\underline{\\mu}_0)\\sim \\displaystyle\\frac{(n-1)p}{n-p}F_{p,n-p}$$ \n", " \n", "\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exemplo: dados banco.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "O arquivo banco.csv, disponível no E-disciplinas, mostra dados de todos os clientes de um banco de acordo com a descrição a seguir.\n", "\n", "- Sexo (1: Feminino, 0: Masculino);\n", "- Idade (em anos);\n", "- CartaodeCredito (1: Sim, 0: Não);\n", "- ChequeEspecial (1: Sim, 0: Não);\n", "- Renda (mensal, em reais);\n", "- LimiteCartaodeCredito (em reais);\n", "- LimiteChequeEspecial (em reais);\n", "- Devedor (1: Sim, 0: Não);\n", "- SaldoDevedor (em reais).\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercício 1**\n", "\n", "1. Considere uma amostra de tamanho 100. \n", "2. Teste se o vetor de médias pode ser considerado $\\underline{\\mu}_0=(0.3, 30, 0.48, 0.43, 3000, 2100, 1240, 0.28, 2150)^T$ usando a amostra.\n", "3. Verifique se $\\underline{\\mu}_0$ está na região de confiança de $\\underline{\\mu}$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import random\n", "\n", "random.seed(123)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "df = pd.read_csv('https://edisciplinas.usp.br/mod/resource/view.php?id=3179488', decimal=',',sep=',', index_col=0)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SexoIdadeCartaodeCreditoChequeEspecialRendaLimiteCartaodeCreditoLimiteChequeEspecialDevedorSaldoDevedor
ID
1026002471.510.000.000.00.00
2125112781.904172.852781.900.00.00
3029002567.380.000.000.00.00
4129103329.504994.240.000.00.00
5119012074.080.002074.081.05887.33
\n", "
" ], "text/plain": [ " Sexo Idade CartaodeCredito ChequeEspecial Renda \\\n", "ID \n", "1 0 26 0 0 2471.51 \n", "2 1 25 1 1 2781.90 \n", "3 0 29 0 0 2567.38 \n", "4 1 29 1 0 3329.50 \n", "5 1 19 0 1 2074.08 \n", "\n", " LimiteCartaodeCredito LimiteChequeEspecial Devedor SaldoDevedor \n", "ID \n", "1 0.00 0.00 0.0 0.00 \n", "2 4172.85 2781.90 0.0 0.00 \n", "3 0.00 0.00 0.0 0.00 \n", "4 4994.24 0.00 0.0 0.00 \n", "5 0.00 2074.08 1.0 5887.33 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Documentação: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f.html\n", "\n", "from scipy.stats import f\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def T2Hotelling(df, mu0, n, p):\n", " Xbarra=df.mean()\n", " S = df.cov()\n", " S_inv = np.linalg.inv(S)\n", " T2Hotelling = n*np.array(Xbarra-mu0).T.dot(S_inv).dot(np.array(Xbarra-mu0))\n", " qf = f.ppf(0.95, p , n-p, loc=0, scale=1)\n", " teste = T2Hotelling > (n-1) * p / (n-p) * qf\n", " pvalor = 1-f.cdf(T2Hotelling/((n-1) * p / (n-p) ), p, n-p)\n", " print('Rejeitamos H0') if teste else print('Não rejeitamos H0')\n", " print('Valor da estatística', T2Hotelling)\n", " print('valor p', pvalor)\n", " " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "df_amostra = df.sample(100)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Não rejeitamos H0\n", "Valor da estatística 14.271963352833856\n", "valor p 0.1758534221481306\n" ] } ], "source": [ "mu0 = [0.3, 30, 0.48, 0.43, 3000, 2100, 1240, 0.28, 2150]\n", "n=len(df_amostra)\n", "p=len(df_amostra.columns)\n", "\n", "T2Hotelling(df_amostra, mu0, n, p)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Sexo 0.3100\n", "Idade 30.2300\n", "CartaodeCredito 0.5100\n", "ChequeEspecial 0.4000\n", "Renda 2987.6827\n", "LimiteCartaodeCredito 2297.4788\n", "LimiteChequeEspecial 1190.2413\n", "Devedor 0.2800\n", "SaldoDevedor 2948.8496\n", "dtype: float64" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_amostra.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Região de confiança**" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Resultado: mu0 está na região de confiança de mu\n" ] } ], "source": [ "Xbarra = df_amostra.mean()\n", "mu0 = [0.3, 30, 0.48, 0.43, 3000, 2100, 1240, 0.28, 2150]\n", "S_inv = np.linalg.inv(df_amostra.cov())\n", "n = len(df_amostra)\n", "\n", "Teste = n*np.array(Xbarra-mu0).T.dot(S_inv).dot(np.array(Xbarra-mu0)) < (n-1) * p / (n-p) * f.ppf(0.95, p , n-p, loc=0, scale=1)\n", "\n", "print('Resultado: mu0 está na região de confiança de mu') if(Teste) else print('Resultado: mu0 não está na região de confiança de mu')\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testes de hipóteses para a comparação de médias em amostras independentes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sejam \n", "\n", "- $\\underline{X}_{11},\\ldots,\\underline{X}_{1n_1}$ vetores aleatórios $p\\times 1$ referentes a uma população com $E(\\underline{X}_{1j})=\\underline{\\mu}_1$ para $j=1,\\ldots,n_1$,\n", "- $\\underline{X}_{21},\\ldots,\\underline{X}_{2n_2}$ vetores aleatórios $p\\times 1$ referentes a uma população com $E(\\underline{X}_{2j})=\\underline{\\mu}_2$ para $j=1,\\ldots,n_2$,\n", "\n", "supondo que a população 1 é independente da população 2.\n", "\n", "\n", "Deseja-se avaliar as hipóteses \n", "\n", "$\\begin{array}{l}H_0:\\underline{\\mu}_1=\\underline{\\mu}_2\\mbox{ contra }\\\\H_1:\\underline{\\mu}_1\\neq\\underline{\\mu}_2\\end{array}$\n", "\n", "Primeiramente, consideramos um estimador para $\\Sigma$, por exemplo\n", "\t\t\t \n", "$S= S_{pooled} = \\displaystyle\\frac{\\displaystyle\\sum_{j=1}^{n_1} (\\underline{X}_{1j}-\\overline{\\underline{X}}_1)(\\underline{X}_{1j}-\\overline{\\underline{X}}_1)^\\top + \\displaystyle\\sum_{j=1}^{n_2}(\\underline{X}_{2j}-\\overline{\\underline{X}}_2)(\\underline{X}_{2j}-\\overline{\\underline{X}}_2)^\\top}{n_1+n_2-2}$\n", "\t\t\t\n", "ou seja\n", "\t\t\t\n", "$S_{pooled} = \\displaystyle\\frac{(n_1-1)S_1 + (n_2-1)S_2}{n_1+n_2-2}$\n", "\t\t\t\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\t\t\t\n", "Então, reescrevemos as hipóteses de interesse na forma mais geral\n", "\t\t\t\n", "$\\begin{array}{l}H_0:\\underline{\\mu}_1-\\underline{\\mu}_2=\\underline{\\delta}_0\\mbox{ contra } \\\\ \\ H_1:\\underline{\\mu}_1-\\underline{\\mu}_2\\neq\\underline{\\delta}_0.\\end{array}$\n", "\t\t\t \n", " \n", "e rejeitamos $H_0$ ao nível de significância $\\alpha$ se\n", "\n", "$T^2_{obs} = (\\overline{\\underline{X}}_1-\\overline{\\underline{X}}_2-\\underline{\\delta}_0)\\left[\\left(\\displaystyle\\frac{1}{n_1}+\\frac{1}{n_2}\\right)S \\right]^{-1}(\\overline{\\underline{X}}_1-\\overline{\\underline{X}}_2-\\underline{\\delta}_0) > c^2$\n", "\n", "com $c^2 = \\displaystyle\\frac{(n_1+n_2-2)p}{n_1+n_2-p-1} q_{F_{p,n_1+n_2-p-1,\\alpha}}.$\n", "\t\t\t\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "Se as hipóteses de interesse são\n", "\n", "\t\t\t\n", "$\\begin{array}{l}H_0:\\underline{\\mu}_1=\\underline{\\mu}_2\\mbox{ contra } \\ H_1:\\underline{\\mu}_1\\neq\\underline{\\mu}_2\\end{array}$\n", "\t\t\t \n", " \n", "então a estatística se simplifica, sob $H_0$, em\n", "\t\t\t\n", " \n", "$T^2 = ( \\overline{\\underline{X}}_1 - \\overline{\\underline{X}}_2 )^\\top \\displaystyle\\frac{S_{pooled}^{-1}}{n_1+n_2-2 } ( \\overline{\\underline{X}}_1 - \\overline{\\underline{X}}_2 ) \\sim \\displaystyle\\frac{(n_1+n_2-2)p}{n_1+n_2-p-1} F_{p,n_1+n_2-p-1}.$ \n", "\t\t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercício 2**\n", "\n", "1. Considere que a primeira população é composta por mulheres e a segunda população é composta por homens. Com base na amostra obtida no Exercício 1, verifique se os vetores de médias podem ser considerados iguais." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Médias amostrais**" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Sexo 0.000000\n", "Idade 30.579710\n", "CartaodeCredito 0.376812\n", "ChequeEspecial 0.188406\n", "Renda 3029.261449\n", "LimiteCartaodeCredito 1728.385797\n", "LimiteChequeEspecial 587.988841\n", "Devedor 0.362319\n", "SaldoDevedor 3241.563043\n", "dtype: float64" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_amostra[df_amostra['Sexo']==0].mean()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Sexo 1.000000\n", "Idade 29.451613\n", "CartaodeCredito 0.806452\n", "ChequeEspecial 0.870968\n", "Renda 2895.136452\n", "LimiteCartaodeCredito 3564.169677\n", "LimiteChequeEspecial 2530.738710\n", "Devedor 0.096774\n", "SaldoDevedor 2297.326129\n", "dtype: float64" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_amostra[df_amostra['Sexo']==1].mean()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABHQAAAI/CAYAAAAWQ44jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZhsVXnv8e+PGUEgCiIC5iABVFCOgAOKCEa9YBxQIYgkBE1EbsAxxpFLiNEkyo0xiKiHXERREQcUJEQwDEIQkPHAOSBChETERFBRcUA4/d4/9moo2q4+U5+u2vD9PE89vWvV2nu/e6jVVavevXaqCkmSJEmSJPXHaqMOQJIkSZIkScvHDh1JkiRJkqSesUNHkiRJkiSpZ+zQkSRJkiRJ6hk7dCRJkiRJknrGDh1JkiRJkqSesUNHYynJu5MsTnJNkquTPH3UMT3UJLlrSPmJSfZdjuXMS7Jo9iIbD0keneRzSf4jyRVJzkyy7XLM/65ZjOXgJMcuQ729k1ye5LokVyX5h5Vc713t72OSfLFNz0/ywpVZ7krEM90xOSTJGaOIp8V0fpIbWjt29eR+WsXrfEmSdyylzjKdM+MoyZK2Lxcl+WqSjWZpucvVts2F6drhJIcmOWg5l/PN9ndeklct4zzbtvfQjUmuTPL5JJsuxzo3SvLnyxPnUpZ3VJK3LkO9g9q5cW1r55Y6zwzLuu//V5JdkvymTe+R5JlteqTHY9Tv5SS3tH092cYdMwfrXOo+X9bzZZbimWyTFidZmOQvkqzS71jj2F7NtuX5LjBsf7T36oyfAdp76PbWXtyY5KzJ9/eqNOxztrS81hh1ANJUSXYFXgTsVFV3J9kYWGvEYUn3SRLgy8Anq+qVrWxHYFPgO8swb4B3AX+7ikMdXO8OwLHAH1TVt5OsDhwyTb01qure5Vl2Vd0GTH6Qmg/sApy5kiEvlxmOyUvmMo4hDqyqy+dqZVV1OnD6XK1vBH5VVfMBknwSOAx432hDmjtV9bEVmGfyy8k84FXAZ2eqn2Qd4F+At1TVV1vZHsAmwP8sbX1J1gA2Av4cOG55411RSfYG3gS8oKpuS7I28Ftf/Fewnbt8skMH2AO4C/jmGByPcbBnVd0xVytbkX2+ig22SY+iO54bAH810qgGJFm9qpaMOo5lNYLvAqdU1eFt3XsCpybZs6quX4XrXGYr0mbpocMMHY2jzYA7qupugKq6o30w2znJN9ov72cl2SzJhu3X7+0Akpyc5LXpHD3wK93+I92iHmv78ti2n/8NeNTAa0cmuazt5wXtSzXtWC1MspDuy9Zk/dXbcbms/eLyurnfolmxJ3DP4IfKqloIXJXknPbr6bVJXgr3/Qp7Q5JPAYuA/wes235x+kyr85V2bi9Ocl9HS5ID2rIWJXn/QPmrk3wnybeAZw2Ub5LkS20fX5Zk8rW3Ae+rqm+3eJdU1UfbPCcm+ViSS4EPJNk6yddaPBcmeXyrt1WSi1s87x1Y57wW31rAe4D927btn+QRbduuSXJJkifP4nEYNOyYXAisn+SLSb6d5DNTztMHtCkD5Qvb4+jc/wv9A34JT3JG+1JFkhe0fXNlki8kWX+mYJPs1/bZwiQXDCz/tHRZPTcm+auB+n+U5Fttv348XYccSfZq61yY5JypcSZ5cZJL0/3y+G9ZjgyLnrgY2BxghvP2xCTHJPlmku+m/YqbznK1beMgA5kH7Vz5x3SZd9cneWqSU9v5M/genfwl+O+BZ7fz6M0Z3ia/Crh4svMAoKrOr6pF7f1+YTvvrsz9mSp7tPLTgevaurZu6zo6yfqZpn1s876l7etFSd40UP7udO3cvwPbDZRPe6yBdwJvbZ3MVNXdVXX8wL76UJLLgTcuy/ufB/7/2gNYJ8k84FDgzW3bTkjyviTnJrkryX+1fTknx6M9fUzbHzcm+cDAeqZtl1q78e1WfkxaBkOmZLW04zGvTU/bBg2T5A3pskGvSfK5geWf1GK6MclrB+r/5cB2//VA+UGtbGGSk6bGme4z32Xt9S8ledhMca1qVfVDuh9LDk9n2mOaLpv0DybnS8swmaF+Mry9+v10bfy17Xxcu5XfkuT9Sa4E9pvL/TALhn0XWGrbPHh+Ay8fKF+mzyNVdR6wgPajV6Zpb9J9//jPtEysJOsl+V6SNaer3+pslek/QyXTfGfJb7ep0vSqyoePsXoA6wNX02U6HAc8B1gT+CawSauzP3BCm34+3Yf6VwJfa2WvAL4OrE6XNfFfwGaj3rY+PYC72t+XD+zLxwB3Avu21x4xUP8k4MVt+hpg9zZ9NLCoTR8CHNGm1wYuB7Ya9bauwL55A/CP05SvAWzQpjcGbqLLxpkHTADPmLp/B54/ov1dl67T55Ftf/8X3a+wawDnAvvQfdCZLF8LuAg4ts3/WWC3Nv1Y4Po2fSWw45DtORE4A1i9PT8H2KZNPx04t02fDhzUpg8bOEfmDRzjgydjac8/DPxVm34ucPUcH5M9gJ8CW9D9iHExsBsztynDzt+p23ZGW/7GwAXAeq387cCRbfp84Aa6Nu1q4OhWfi2weZveaGD5P2jHfvI82AV4AvBVYM1W7zi6rINNgO/R3kMD59B9cQK/A6RN/xnwD9NtS58eA+fd6sAXgL2Wct6e2OqtBjwRuKmVL3fbNqptnVJ2FF2nxeT59f42/UbgNrr2YW3gVuCRU/bZHsAZA8uatk0GPgi8cUhMDwPWadPbAJcPLPsXA+fjPNp7pz0f1j7u3N4P69H9/18MPGWg/GF02Q43DWz3sGP9Y2DDIXGfDxzXplfk/b8HcO80x+CoFvOftHWcCXxlDo/HwcB3gQ2BdYD/BLZkSLvU6nyvHbsAn5+MYXC72vNF7ThO2wa16VvacZps497cym8D1m7TGw0sfyFd+7Zxi+MxwAvovkCH7n16BrA7sD3dZ8GNB9+XU/b/IwfifS/w+um2ZQTv0zvpPn8OO6Yvo8sohe7/+PfafhlWf9r2auB4btvm+RTwpoFj87ZRtF2zsE9/67vA4DnQpgc/d544ZX9Md35P+3mEaf4f0n3W+tc2Pay9OY0uOw26NuSfl1J/2Geoab+zMKVN9eFj2MNLrjR2ququJDsDz6b71f0Uun/SOwBfb53xq9N98aGqvp5kP+AjwI5tMbsBJ1eXXvo/Sb4BPJUH92UIq8ru3L8vb0ty7sBreyZ5G90H7kcAi5NcSPfh7YJW5yRg7zb9AuDJuf865w3p/unevKo3Yo4E+Nsku9N14GxO988Z4D+r6pIZ5n1Dkpe16S3p9sumwPlVdTtAumye3VudwfJTgMnxe54HPHHgR6sNspRskeYLVbWk1X0m8IWBZazd/j6L7oMHdMf1/SzdbpPzVNW5SR6ZZIOq+tkyzDtbvlVVtwIkuZruC8qdTNOmpBuPZdj5O8wz6DoKLmrLWouu42jSdJdcXQScmOTzwKkD5V+vqh+1WE+l23/30n25vawtf13gh229F1TVzQBV9eNpYtsCOCVd9sFaPDjea+u247g5cD3dMZzpvAX4SlVNANfl/iyl5Wrb6L7QjqPJ/2vXAour6gcASb5L15b8aIZ5h7XJM1kTODbJfGAJ97c90L3Xhp1jw9rH3YAvV9UvWtyn0v3/X62V/7KVn97+Lu1Yz+SU9nc7Zu/9D12b8lng1XRZN6cCxzA3xwPgnKr6aVvPdcDv0l3yNl279Hjg5qq6sdX/NNNcfjvF7zN9GzRpukuurgE+k+QrdB1ck06rql8Bv0pyHvA0unPgBcBVrc76bbt3pPvfdAcMbeN2aNkOG7X5zlrKtsy1Ycf0X4F/atk0e9G15b9KMqz+sPZqO7rjOXm59+RlqB9qzyfP+V6Z7rtAurHhfr6Utnmm83vazyNDQpjM5J2pvTmFriPnPLoflY9bwc9Qw76z/IyZ21QJcAwdjanWqJ0PnJ/kWrp/TouratepdVu64xOAX9L9Gn3rHIb6kJXumv7jgF2q6ntJjqL7ZWTG2eh+PRu3D1zLazH3jxkz6EC6rImdq+qeJLdw/z75xbCFpUvlfx6wa1X9Msn5LH1fDrMaXSbQr6esYzHdB/KFQ+abjG814M5q4wFMo1YwrlVt2DEBuHtgegnd/74wTZuSmQfYvZcHXqo8eYxC1xFzwLIGW1WHphvg8Q+AK9oHV/jt/Vtt+Z+sqndOifXFy7CqDwMfrKrT23l21LLGOMZ+VVXz011acRbd/4cTmfm8HTwHZrx8agXbtlGa3LYJHridEyz9c960bXKSLemyY6fzZrpxdHakez8MtjVD2zlmbh+Xx0xt1GQ7d+40rw3GtyLv/2U118djcJ3wwDbut9ql1hE3zExt3G+1QUvxB3SdEC8G3p3kSa18WBv3d1X18Smxvn4Z1nMisE9VLUxyMF1Ww0gleRzdcfghM3zuaf/r/xddp8DnJounq58Vv9nATO/JsTbNd4HXAU9mbtrmp9D9YDBTe3M6XSf1I7i/3Vlvhvqw/J+henv8NHccQ0djJ8l2SQZ/kZpP16hukm6QNNo1qtu319/cXn8V8Ikka9KNm7F/umuRN6H7UPGtOduIB5cLuH9fbkb3Swnc/0/0jvaLxL4AVXUncGeS3drrBw4s6yzgf7djRLq7dqy3yrdg9p0LrJ0HjnXzZLpfRX/Yvqzs2Z4Pc8/kfqD7Be4nrTPn8XSZF9Cds89JsnG68QoOAL4BXNrKH9mWMXht/NnAfR+CBz68Hw28K+1OXElWS3Lo1KBa5szNLett8truycy3i+h+hYIHHtdBPwcePvD8wsm6rUPhjlWUnTPsmDx7SP0bmKZNWcr5ewswv+27Lel+WQa4BHhWkt9ry1ovS7njWZKtq+rSqjoSuJ3ul3uA56e7zn9dupTvi+jSt/dNN9jm5DgAv9vWu3uSrSbLp1nVhsD32/SfzBRT37TMjTcAf0HXoT/svB1mudq2B4mp789hbfJngWfmgWN87J5ucPUNgR+0jKc/pstuWZZ1bcj07eOFwD5JHtbW/bJWdkErXzfJw+k6BpbWRv0dcHSSR7fX1kryZ9PEtiLv/5m27Rbubxv3bvEvi9k4HsMMa5e+DcxLsnWrN9jhcwuwU6u/E92lPjC8DZpWuh/atqxuLJK30x37yUzRlyZZJ8kj6TpfLmvb/ZrcP8bP5m1d5wL7tbrD2riH02VXrcnw4zVn2mfOj9FdwlPM/LnnFLqsrmcDX2tlw+oPa69uoDuev9ee/zHd54ReG/Jd4IY2PVPbPNP5vUyfR5I8hy6r5/iZ2puquovu/P0nusu6lqzgZyi/s2ilmKGjcbQ+8OF0v5TdS3fd/CF011cfk2RDunP3Q0nupRsX4mlV9fN0g4seQfcr9K502QhFdw3xf8/5ljw4fJnuWuPr6K7rvRi6jpskx9NdY//fdP/UJr0aOCFJ0XUwTPpnutT0K5OE7ovsPqt6A2ZbVVW6y6M+lOTtdL9Q30J33h2T7peky+k+WAyzALgm3aB9rwEOTXI93QeWS9p6fpAuxfg8ul/t/qWqToNuUEi6Y3En3XXmk94AfCTJNXTvkwuAQ6vqmnSDjZ6cLrOh6MYpmM6BwEeTHEF3ecXn6N5LbwQ+27b5tCHznge8I90lMX/X9skJLZ5fsoo6FWY4Jl8ZUv836VLaH9Cm0P3CP+z8vYjukqXr6DqRr2zLuj3dL8Mnpw1GSdcOTabAfybJr9r0HVX1PLovnZPX+J9Dt3/n032I+xLdpVKfnrxUqx2Ls9sXpXuAw6rqknQdWKe28h/SjSk26Ci61O+f0H052ooHkaq6qp1bBzD8vB1mRdq2ufawJINZpx9cyeVdAyxJN+DviXRfROYxpU2uqp8meRHd++lDdOfcNXRtwHHAl9LdNvprDPkFuap+lOSidIOK/yvd5QVfndo+VtWVSU7k/i8w/1xVV8F9l5MupDu3B4/DtMe6qs5Md0ndv7XtKeCEaWJbkfc/wOrteKwBPKJ1Fp1Fdy69mm7Mq4fRfdF83HT7ZYrZOB7TGtYuVdV3WrvxL0l+SfdlcrJT6UvAQekyOi+ltWFVdd10bRDdeD0A5yWZvIPSNcCfAp9u+zbAMe19Nfn6eXRj6PxNdQNY35bkCcDFrc5dwB9V1eIk7wO+0ZZ/Fd14J4P+T4v19vb34cy9yctA16T73HoS979XZ/rcc3are1pV/WYp9Ye1V79O8mq6dn4NuvfJuN0FbEUM+y5wJzO0zW1/DDu/j2L455H903XkPozu//wr6v47XM30v+UUujHa9hhY1vJ+hvoy03xnyf2DvUszmhwoUZIkTZHuDi9nVNVMv4TP1roOpkslP3xVr0uS4L5MhbdW1YvmYF1H0Q0E+39X9bok6aHCS64kSZIkSZJ6xgwdSZIkSZKknjFDR5IkSZIkqWfs0JEkSZIkSeoZO3QkSZIkSZJWQpITkvyw3eFxuteT5JgkNyW5JslOK7tOO3TUa+3WhBoTHo/x4bEYHx6L8eGxGB8ei/Hi8RgfHovx4bHQCjgR2GuG1/cGtmmPQ4CPruwK7dBR39nQjhePx/jwWIwPj8X48FiMD4/FePF4jA+PxfjwWGi5VNUFwI9nqPJS4FPVuQTYKMlmK7NOO3QkSZIkSZJWrc2B7w08v7WVrbA1Vioc9do9d3y39/esP+4f3tv77Xj84/cddQizZuP1tmTrjXfq9fE4bYPHjDqEWXHkI5/Eose9qNfHYqsD1h11CLPimL134hfv3q/Xx+IpC24edQiz4lHrP5ZtN9ml18fis2tvMeoQZsU7N9qRy7fYp9fHYrt9fjPqEGbNP+2xPT8//IW9Ph6/s2DhqEOYFVl9Q9ZYa/NeH4tnP+qJow5hVmy70ePZc4vn9/pYAJx369cz6hhWlbn8HrjWJlu/jgdmbS2oqgVztf5h7NBRr/3ZQQeMOgQN2GCdjUcdgpr9NnjsqENQ85qnPG7UIajZcJ1NRh2CmpevN2/UIWjAq7f3f8a4WG219UYdgprHrPfg6EDX7GidNyvbgfN9YMuB51u0shXmJVeSJEmSJEmr1unAQe1uV88AflpVP1iZBZqhI0mSJEmS+mViyagjeIAkJwN7ABsnuRX4K2BNgKr6GHAm8ELgJuCXwKtXdp126EiSJEmSJK2EqppxPJCqKuCw2VynHTqSJEmSJKlfamLUEYycY+hIkiRJkiT1jBk6kiRJkiSpXybM0DFDR5IkSZIkqWfM0JEkSZIkSb1SjqFjho4kSZIkSVLfmKEjSZIkSZL6xTF0zNCRJEmSJEnqGzN0JEmSJElSvziGjhk6kiRJkiRJfWOHjiRJkiRJUs94yZUkSZIkSeqXiSWjjmDkzNCRJEmSJEnqGTN0JEmSJElSvzgoshk6kiRJkiRJfWOGjiRJkiRJ6pcJM3TM0JEkSZIkSeoZM3QkSZIkSVKvlGPomKEjSZIkSZLUN2boSJIkSZKkfnEMHTN0JEmSJEmS+sYMnTGQ5N3Aq4AlwATwuqq6dLRRSZIkSZI0phxDxw6dUUuyK/AiYKequjvJxsBaIw5LkiRJkiSNMS+5Gr3NgDuq6m6Aqrqjqm5LsnOSbyS5IslZSTZLsmGSG5JsB5Dk5CSvTefoJIuSXJtk/5FukSRJkiRJq9LEkrl7jCk7dEbvbGDLJN9JclyS5yRZE/gwsG9V7QycALyvqn4KHA6cmOSVwO9U1fHAy4H5wI7A84Cjk2w2kq2RJEmSJEmrnB06I1ZVdwE7A4cAtwOnAK8DdgC+nuRq4Ahgi1b/68C1wEeAP2uL2Q04uaqWVNX/AN8Anjrd+pIckuTyJJf/86dOXnUbJkmSJEnSqlITc/cYU46hMwaqaglwPnB+kmuBw4DFVbXr1LpJVgOeAPwS+B3g1uVc1wJgAcA9d3y3Vi5ySZIkSZI0CmbojFiS7ZJsM1A0H7ge2KQNmEySNZNs315/c3v9VcAn2uVZFwL7J1k9ySbA7sC35mwjJEmSJEnSnDJDZ/TWBz6cZCPgXuAmusuvFgDHJNmQ7jh9KMm9dJdZPa2qfp7kArrLsY4CdgUWAgW8rar+e863RJIkSZKkuTAxvpdCzRU7dEasqq4AnjnNS3fQZdpM9YSBed8yUP6X7SFJkiRJkh7k7NCRJEmSJEn9MsaDFc8Vx9CRJEmSJEnqGTN0JEmSJElSvziGjhk6kiRJkiRJfWOGjiRJkiRJ6pWqJaMOYeTM0JEkSZIkSeoZM3QkSZIkSVK/eJcrM3QkSZIkSZL6xgwdSZIkSZLUL97lygwdSZIkSZKkvjFDR5IkSZIk9Ytj6JihI0mSJEmS1Ddm6EiSJEmSpH6ZWDLqCEbODB1JkiRJkqSesUNHkiRJkiSpZ7zkSpIkSZIk9YuDIpuhI0mSJEmS1Ddm6EiSJEmSpH6ZMEPHDB1JkiRJkqSeMUNHkiRJkiT1i2PomKEjSZIkSZLUN2boSJIkSZKkfnEMHTN0JEmSJEmS+sYMHUmSJEmS1C9m6JihI0mSJEmS1Ddm6DyEPf7x+446BAHf/vYXRx2CBrxul7eNOgQ1P/v03aMOQc1Vb9l+1CGo+eDH7h11CGpOOmOdUYegAR/adM9Rh6DGjAHNlaolow5h5Hy/SZIkSZIk9YwZOpIkSZIkqV8cQ8cMHUmSJEmSpL4xQ0eSJEmSJPVLmaFjho4kSZIkSVLP2KEjSZIkSZLUM15yJUmSJEmS+sVBkc3QkSRJkiRJ6hszdCRJkiRJUr84KLIZOpIkSZIkSX1jho4kSZIkSeoXx9AxQ0eSJEmSJKlvzNCRJEmSJEn94hg6ZuhIkiRJkiT1jRk6kiRJkiSpXxxDxwwdSZIkSZKkvjFDR5IkSZIk9YsZOmboSJIkSZIk9Y0ZOpIkSZIkqV+8y5UZOpIkSZIkSX1jho4kSZIkSeoXx9AxQ0eSJEmSJGllJNkryQ1Jbkryjmlef2yS85JcleSaJC9c2XXaoSNJkiRJkrSCkqwOfATYG3gicECSJ06pdgTw+ap6CvBK4LiVXa+XXEmSJEmSpH4Zr0GRnwbcVFXfBUjyOeClwHUDdQrYoE1vCNy2siu1Q0eSJEmSJGnFbQ58b+D5rcDTp9Q5Cjg7yeuB9YDnrexKveRKkiRJkiT1y8TEnD2SHJLk8oHHISsQ8QHAiVW1BfBC4KQkK9UnY4fOCCS5a0j5iUn2XY7lzEuyaPYikyRJkiRJg6pqQVXtMvBYMKXK94EtB55v0coG/Snw+ba8i4F1gI1XJi47dCRJkiRJUr/UxNw9lu4yYJskWyVZi27Q49On1Pkv4PcBkjyBrkPn9pXZBXbojFA6x7Zbm/0b8KiB145MclmSRUkWJEkr3znJwiQLgcMG6q+e5Og2zzVJXjf3WyRJkiRJ0kNLVd0LHA6cBVxPdzerxUnek+QlrdpfAK9t3+VPBg6uqlqZ9Too8mi9DNiO7rZmm9KNgH1Ce+3YqnoPQJKTgBcBXwU+ARxeVRckOXpgWX8K/LSqnppkbeCiJGdX1c1ztC2SJEmSJM2NibG6yxVVdSZw5pSyIwemrwOeNZvrNENntHYHTq6qJVV1G3DuwGt7Jrk0ybXAc4Htk2wEbFRVF7Q6Jw3UfwFwUJKrgUuBRwLbTF3h4GBOP/v1HatimyRJkiRJ0ipmhs4YSrIOcBywS1V9L8lRdNfXzTgb8PqqOmumSm3wpgUAW2+800qld0mSJEmSNBJjlqEzCmbojNYFwP5t/JvNgD1b+WTnzR1J1gf2BaiqO4E7k+zWXj9wYFlnAf87yZoASbZNst4q3wJJkiRJkjTnzNAZrS/TXU51Hd2I1xdD13GT5HhgEfDfdCNmT3o1cEKSAs4eKP9nYB5wZRtA+XZgn1W9AZIkSZIkzbmVG0/4QcEOnRGoqvXb36IbCXu6OkcAR0xTfgWw40DR21r5BPCu9pAkSZIkSQ9iduhIkiRJkqR+cQwdx9CRJEmSJEnqGzN0JEmSJElSv5ihY4aOJEmSJElS35ihI0mSJEmS+qXM0DFDR5IkSZIkqWfs0JEkSZIkSeoZL7mSJEmSJEn94qDIZuhIkiRJkiT1jRk6kiRJkiSpX6pGHcHImaEjSZIkSZLUM2boSJIkSZKkfnEMHTN0JEmSJEmS+sYMHUmSJEmS1C9m6JihI0mSJEmS1Ddm6EiSJEmSpH4pM3TM0JEkSZIkSeoZM3QkSZIkSVKv1ESNOoSRM0NHkiRJkiSpZ8zQkSRJkiRJ/eJdrszQkSRJkiRJ6hszdCRJkiRJUr94lyszdCRJkiRJkvrGDh1JkiRJkqSe8ZIrSZIkSZLUL9623A6dh7LTNnjMqEMQ8Lpd3jbqEDTg45d/YNQhqDl/+3eOOgQ1f/7xu0YdgprjL37PqENQc+FTPBbj5L1r/s+oQ1DzrnseMeoQpIcMO3QkSZIkSVK/eNtyx9CRJEmSJEnqGzN0JEmSJElSv5ihY4aOJEmSJElS35ihI0mSJEmS+qW8y5UZOpIkSZIkST1jho4kSZIkSeoXx9AxQ0eSJEmSJKlvzNCRJEmSJEn9MuEYOmboSJIkSZIk9YwZOpIkSZIkqV/KMXTM0JEkSZIkSeoZM3QkSZIkSVK/OIaOGTqSJEmSJEl9Y4eOJEmSJElSz3jJlSRJkiRJ6pWacFBkM3QkSZIkSZJ6xgwdSZIkSZLULw6KbIaOJEmSJElS35ihI0mSJEmS+qUcQ8cMHUmSJEmSpJ4xQ0eSJEmSJPWLY+iYoSNJkiRJktQ3ZuhIkiRJkqR+mXAMHTN0JEmSJEmSembsO3SSPDrJ55L8R5IrkpyZZNvlmP9dsxjLwUmOXYZ6eye5PMl1Sa5K8g8rud672t/HJPlim56f5IUrs1xJkiRJknppoubuMabGukMnSYAvA+dX1dZVtTPwTmDTZZk3yWrArHXoLIskOwDHAn9UVU8EdgFumqbecl/uVlW3VdW+7el8wA4dSZIkSZIegsa6QwfYE7inqj42WVBVC4GrkpyT5Mok1yZ5KUCSeUluSPIpYBHw/4B1k1yd5DOtzldaps/iJIdMLjfJAW1Zi5K8f6D81Um+k+RbwLMGyjdJ8qUkl7XH5GtvA95XVd9u8S6pqo+2eU5M8rEklwIfSLJ1kq+1eC5M8vhWb6skF7d43juwznktvrWA9wD7t23bP8kj2gHBek4AACAASURBVLZdk+SSJE+exeMgSZIkSdL4qIm5e4ypcR8UeQfgimnKfw28rKp+lmRj4JIkp7fXtgH+pKouAUiyX1XNH5j3NVX14yTrApcl+RKwNvB+YGfgJ8DZSfYBLgX+upX/FDgPuKot55+Af6yqf0/yWOAs4Akt5pkusdoCeGZVLUlyDnBoVd2Y5OnAccBz27I/WlWfSnLY1AVU1W+SHAnsUlWHt+38MHBVVe2T5LnAp+iyeCRJkiRJ0oPMuHfoDBPgb5PsDkwAm3P/ZVj/OdmZM8QbkrysTW9J1wG0Kd1lXbcDtGye3VudwfJTgMnxe54HPLG7KgyADZKsvwyxf6F15qwPPBP4wsAy1m5/nwW8ok2fRNfZtDS7Tc5TVecmeWSSDarqZ4OVWlbSIQBHPvJJ7LfBY5dh0ZIkSZIkjZExHttmrox7h85iYN9pyg8ENgF2rqp7ktwCrNNe+8WwhSXZg64jZteq+mWS8wfmW16rAc+oql9PWcdiuoyehUPmm4xvNeDOKdlDg1bJ2VlVC4AFAIse9yLfAZIkSZIk9dC4j6FzLrD2lLFungz8LvDD1pmzZ3s+zD1J1mzTGwI/aZ05jwee0cq/BTwnycZJVgcOAL5Bd8nVc1q2y5rAfgPLPRt4/UBckx0zRwPvmrwTV5LVkhw6NaiWOXNzkv1avSTZsb18EfDKNn3gkO36OfDwgecXTtZtHVd3TM3OkSRJkiRJDw5j3aFTVQW8DHheu235YuDvgDOBXZJcCxwEfHuGxSwArmmXUX0NWCPJ9cDfA5e09fwAeAfdGDkLgSuq6rRWfhRwMV0ny/UDy31Di+GaJNcBh7ZlXQO8CTi5rWcR8LghsR0I/GmShXTZSC9t5W8EDmvbt/mQec+ju+Tr6iT7tzh3TnJN27Y/mWGfSJIkSZLUWzUxMWePcTXul1xRVbcBfzjNS7sOmWWHKfO/HXj7QNHeQ9ZzMnDyNOWfAD4xTfkdwP5DlnUGcMY05QdPeX4zsNc09W7mgdt3RCu/hbZ9VfVj4KlTZt1nungkSZIkSdKDy9h36EiSJEmSJD2AgyKP9yVXkiRJkiRJ+m1m6EiSJEmSpH4xQ8cMHUmSJEmSpL4xQ0eSJEmSJPVLje/dp+aKGTqSJEmSJEk9Y4aOJEmSJEnqF8fQMUNHkiRJkiSpb8zQkSRJkiRJvVJm6JihI0mSJEmS1Ddm6EiSJEmSpH4xQ8cMHUmSJEmSpJWRZK8kNyS5Kck7htT5wyTXJVmc5LMru04zdCRJkiRJUr9MTIw6gvskWR34CPB84FbgsiSnV9V1A3W2Ad4JPKuqfpLkUSu7XjN0JEmSJEmSVtzTgJuq6rtV9Rvgc8BLp9R5LfCRqvoJQFX9cGVXaoeOJEmSJEnSitsc+N7A81tb2aBtgW2TXJTkkiR7rexKveRKkiRJkiT1yxwOipzkEOCQgaIFVbVgORezBrANsAewBXBBkidV1Z0rGpcdOpIkSZIkSUO0zpuZOnC+D2w58HyLVjboVuDSqroHuDnJd+g6eC5b0bi85EqSJEmSJPXLRM3dY+kuA7ZJslWStYBXAqdPqfMVuuwckmxMdwnWd1dmF9ihI0mSJEmStIKq6l7gcOAs4Hrg81W1OMl7krykVTsL+FGS64DzgL+sqh+tzHq95EqSJEmSJPVK1dyNobMsqupM4MwpZUcOTBfwlvaYFWboSJIkSZIk9YwZOpIkSZIkqV/m8C5X48oMHUmSJEmSpJ4xQ0eSJEmSJPWLGTpm6EiSJEmSJPWNGToPYVsdsO6oQxDws0/fPeoQNOD87d856hDU7LH470Ydgpp/2un1ow5BzblP+ZtRh6Dm968+atQhaMD3n/zHow5BzX+ss+moQ9CA5486gFWozNAxQ0eSJEmSJKlvzNCRJEmSJEn9YoaOGTqSJEmSJEl9Y4aOJEmSJEnql4lRBzB6ZuhIkiRJkiT1jB06kiRJkiRJPeMlV5IkSZIkqVe8bbkZOpIkSZIkSb1jho4kSZIkSeoXM3TM0JEkSZIkSeobM3QkSZIkSVK/eNtyM3QkSZIkSZL6xgwdSZIkSZLUK97lygwdSZIkSZKk3jFDR5IkSZIk9Ytj6JihI0mSJEmS1Ddm6EiSJEmSpF5xDB0zdCRJkiRJknrHDB1JkiRJktQvjqFjho4kSZIkSVLfmKEjSZIkSZJ6pczQMUNHkiRJkiSpb+zQkSRJkiRJ6hkvuZIkSZIkSf3iJVdm6EiSJEmSJPWNGTqSJEmSJKlXHBS5Rxk6SR6d5HNJ/iPJFUnOTHJIkjNGGNP5SW5IcnV7fHEO1vmSJO9YSp2Dkxy7qmORJEmSJEmj0YsMnSQBvgx8sqpe2cp2BF4y0sA6B1bV5XO1sqo6HTh9rtYnSZIkSdLYMUOnNxk6ewL3VNXHJguqaiFwIbB+ki8m+XaSz7TOH5LsnOQbLZvnrCSbDZQvbI+jkyxq5Q/IaklyRpI92vQLklyc5MokX0iy/kzBJtkvyaK2jgsGln9ay+q5MclfDdT/oyTfalk+H0+yeivfq61zYZJzpsaZ5MVJLk1yVZJ/S7Lpyu9qSZIkSZI07vrSobMDcMWQ154CvAl4IvA44FlJ1gQ+DOxbVTsDJwDva/U/Aby+qnZclhUn2Rg4AnheVe0EXA68ZaDKZwYuuTq6lR0J/K+2jsEsoqcBrwCeDOyXZJckTwD2B55VVfOBJcCBSTYBjgde0Zaz3zTh/TvwjKp6CvA54G3Lsk2SJEmSJPVZTczdY1z14pKrpfhWVd0KkORqYB5wJ10n0Ndbws7qwA+SbARsVFUXtHlPAvZeyvKfQddZdFFb1lrAxQOvT3fJ1UXAiUk+D5w6UP71qvpRi/VUYDfgXmBn4LK2/HWBH7b1XlBVNwNU1Y+niW0L4JSWfbQWcPNStoUkhwCHAByz90685imPW9oskiRJkiRpzPSlQ2cxsO+Q1+4emF5Ct00BFlfVroMVW4fOMPfywIyldSZno+uIOWBZg62qQ5M8HfgD4IokO0++NLVqW/4nq+qdU2J98TKs6sPAB6vq9HZ52FHLENsCYAHAL96939R4JEmSJEkae+OcOTNX+nLJ1bnA2i27BIAkTwaePaT+DcAmSXZtdddMsn1V3QncmWS3Vu/AgXluAeYnWS3JlnSXRwFcQncZ1++1Za2XZNuZgk2ydVVdWlVHArcDW7aXnp/kEUnWBfahy+Q5B9g3yaPavI9I8rttvbsn2WqyfJpVbQh8v03/yUwxSZIkSZKkB49eZOhUVSV5GfChJG8Hfk3XAfOVIfV/k2Rf4JgkG9Jt54foMn1eDZyQpICzB2a7iO6SpeuA64Er27JuT3IwcHKStVvdI4DvtOnPJPlVm76jqp4HHJ1kG7rsm3OAhcB84FvAl+gulfr05KVaSY4Azk6yGnAPcFhVXdI6sE5t5T8Enj9lU48CvpDkJ3SdXlstdWdKkiRJktRzZuj0pEMHoKpuA/5wmpeOH6hz+MD01cDu0yznCmBHgCTzgBe28uKBGTuD85wLPHWa8j2G1H/51LI2Ps6tVbXPNPVPAU6ZpvxfgX+dUnYicGKbPg04bZr57qsjSZIkSZIefHrToSNJkiRJkgRAZdQRjNxDukOnqm6huxvWXKzrRMyakSRJkiRJs+Ah3aEjSZIkSZL6xzF0+nOXK0mSJEmSJDV26EiSJEmSJPWMl1xJkiRJkqReqQkHRTZDR5IkSZIkqWfM0JEkSZIkSb3ioMhm6EiSJEmSJPWOGTqSJEmSJKlXqhxDxwwdSZIkSZKknjFDR5IkSZIk9Ypj6JihI0mSJEmS1Dtm6EiSJEmSpF6pCcfQMUNHkiRJkiSpZ8zQkSRJkiRJvVI16ghGzwwdSZIkSZKknjFDR5IkSZIk9Ypj6JihI0mSJEmS1Dtm6EiSJEmSpF4xQ8cMHUmSJEmSpN6xQ0eSJEmSJKlnvORKkiRJkiT1irctN0NHkiRJkiSpd1J2az1kbbvJLh78MXDVW7YfdQga8Ocfv2vUIaj50cSvRx2Cmi9f+eFRh6DmhPlHjjoENTeusWTUIWjA9vesPuoQ1DhM7Xg5+PufftAeku8+6QVz9n32cdeePZb70QwdSZIkSZKknnEMHUmSJEmS1CtVY5k0M6fM0JEkSZIkSeoZM3QkSZIkSVKv1MSoIxg9M3QkSZIkSZJ6xgwdSZIkSZLUKxOOoWOGjiRJkiRJUt+YoSNJkiRJknrFu1yZoSNJkiRJktQ7ZuhIkiRJkqReqQkzdMzQkSRJkiRJWglJ9kpyQ5KbkrxjhnqvSFJJdlnZdZqhI0mSJEmSeqVq1BHcL8nqwEeA5wO3ApclOb2qrptS7+HAG4FLZ2O9ZuhIkiRJkiStuKcBN1XVd6vqN8DngJdOU+9vgPcDv56NldqhI0mSJEmStOI2B7438PzWVnafJDsBW1bVv8zWSr3kSpIkSZIk9cpcDoqc5BDgkIGiBVW1YDnmXw34IHDwbMZlh44kSZIkSdIQrfNmpg6c7wNbDjzfopVNejiwA3B+EoBHA6cneUlVXb6icdmhI0mSJEmSemWixuq25ZcB2yTZiq4j55XAqyZfrKqfAhtPPk9yPvDWlenMAcfQkSRJkiRJWmFVdS9wOHAWcD3w+apanOQ9SV6yqtZrho4kSZIkSeqVGq8MHarqTODMKWVHDqm7x2ys0wwdSZIkSZKknjFDR5IkSZIk9UrVqCMYPTN0JEmSJEmSesYMHUmSJEmS1CtjdperkTBDR5IkSZIkqWfM0JEkSZIkSb0ybne5GgUzdCRJkiRJknrGDp05kGRJkquTLEry1SQbzdJyT0yy72wsS5IkSZKkvqiau8e4skNnbvyqquZX1Q7Aj4HDRh2QJEmSJEnqL8fQmXsXA08GSLI18BFgE+CXwGur6ttJTgR+BuwCPBp4W1V9MUmADwPPB74H/GZyoUmOBF4MrAt8E3hd1Tj3JUqSJEmStGK8y5UZOnMqyerA7wOnt6IFwOuramfgrcBxA9U3A3YDXgT8fSt7GbAd8ETgIOCZA/WPraqntiygddt8kiRJkiTpQcgMnbmxbpKrgc2B64GvJ1mfrkPmC13iDQBrD8zzlaqaAK5Lsmkr2x04uaqWALclOXeg/p5J3gY8DHgEsBj46tRAkhwCHALwqPUfy4brbDJb2yhJkiRJ0pzwLldm6MyVX1XVfOB3gdCNobMacGcbW2fy8YSBee4emJ7xTE2yDl12z75V9STgeGCd6epW1YKq2qWqdrEzR5IkSZKkfrJDZw5V1S+BNwB/QTdmzs1J9gNIZ8elLOICYP8kqyfZDNizlU923tzRMn+885UkSZIkSQ9iXnI1x6rqqiTXAAcABwIfTXIEsCbwOWDhDLN/GXgucB3wX3QDLFNVdyY5HlgE/Ddw2arbAkmSJEmSRstBke3QmRNVtf6U5y8eeLrXNPUPnm7+dteqw4es4wjgiJWNVZIkSZIkjT87dCRJkiRJUq/UqAMYA46hI0mSJEmS1DNm6EiSJEmSpF5xDB0zdCRJkiRJknrHDB1JkiRJktQrZYaOGTqSJEmSJEl9Y4aOJEmSJEnqlYlRBzAGzNCRJEmSJEnqGTN0JEmSJElSrxSOoWOGjiRJkiRJUs+YoSNJkiRJknplokYdweiZoSNJkiRJktQzZuhIkiRJkqRemXAMHTN0JEmSJEmS+sYOHUmSJEmSpJ7xkitJkiRJktQr3rbcDB1JkiRJkqTeMUNHkiRJkiT1ysSoAxgDZuhIkiRJkiT1jBk6kiRJkiSpVxxDxwwdSZIkSZKk3jFDR5IkSZIk9Ypj6JihI0mSJEmS1Dtm6EiSJEmSpF4xQ8cMHUmSJEmSpN4xQ+ch7LNrbzHqEAR88GP3jjoEDTj+4veMOgQ15z7lb0YdgpoT5h856hDUvOZq26hxcfEObx91CBrwgbV+NuoQ1Pyfe9cedQh6iPAuV2boSJIkSZIk9Y4ZOpIkSZIkqVcmTNAxQ0eSJEmSJKlvzNCRJEmSJEm9MuEYOmboSJIkSZIk9Y0dOpIkSZIkST3jJVeSJEmSJKlXatQBjAEzdCRJkiRJknrGDB1JkiRJktQrE6MOYAyYoSNJkiRJktQzZuhIkiRJkqRemYi3LTdDR5IkSZIkqWfM0JEkSZIkSb3iXa7M0JEkSZIkSeodM3QkSZIkSVKveJcrM3QkSZIkSZJ6xwwdSZIkSZLUKxPe5MoMHUmSJEmSpL4xQ0eSJEmSJPXKBKbomKEjSZIkSZLUM2boSJIkSZKkXqlRBzAGzNCRJEmSJEnqGTt0JEmSJEmSesZLriRJkiRJUq9423IzdCRJkiRJknpnqR06Se6apuzQJActz4qSfLP9nZfkVcs4z7ZJzkxyY5Irk3w+yabLsc6Nkvz58sS5lOUdleSty1DvoCSLklyb5KplmWeGZc1LsqhN75LkmDa9R5JnruhyJUmSJEnqq4k5fIyrFcrQqaqPVdWnlnOeyc6HecBSO3SSrAP8C/DRqtqmqnYCjgM2WZb1JVkD2AiYtQ6dZVzv3sCbgBdU1ZOAZwA/HRLfcqmqy6vqDe3pHoAdOpIkSZIkPQStUIfOYKZKkvOT/GOSy5Ncn+SpSU5tWTXvHZhnMtPn74FnJ7k6yZuTrJ7k6CSXJbkmyetavVcBF1fVVyeXUVXnV9WilrVyYcvauXIyU6VlrVyY5HTguraurdu6jk6yfpJz2jzXJnnpQHxvaVk1i5K8aaD83Um+k+Tfge0GyrdO8rUkV7R1Pr699E7grVV1W4v57qo6fmBffSjJ5cAbk+yc5BttGWcl2azV2znJwiQLgcMG1rlHkjOSzAMOBd7ctu3ZbZ+c2/bhOUkeuyLHVpIkSZKkcVdz+BhXszUo8m+qapckbwROA3YGfgz8R5J/rKofDdR9B12Hx4sAkhwC/LSqnppkbeCiJGcDOwBXDFnfD4HnV9Wvk2wDnAzs0l7bCdihqm5uHR87VNX8tq41gJdV1c+SbAxc0jp/dgJeDTwdCHBpkm/QdXi9EphPt6+uHIhpAXBoVd2Y5Ol02UPPXUrcAGu1fbUm8A3gpVV1e5L9gfcBrwE+ARxeVRckOXrqAqrqliQfA+6qqv/btu2rwCer6pNJXgMcA+wzQxySJEmSJKmnZqtD5/T291pgcVX9ACDJd4EtgR8NmxF4AfDkJPu25xsC2yxlfWsCxyaZDywBth147VtVdfOQ+QL8bZLd6S6F2xzYFNgN+HJV/aLFfSrwbLoOnS9X1S9b+ent7/p0lzt9IblvaO21lxLzpFPa3+3oOn++3paxOvCDJBsBG1XVBa3eScDey7DcXYGXD8zzgekqtQ60QwDeudGOvHy9ecsYtiRJkiRJ48G7XM1eh87d7e/EwPTk86WtI8Drq+qsBxQmWwLPGTLPm4H/AXak63T59cBrv5hhXQfSjcGzc1Xdk+QWYJ2lxDed1YA7JzN/plhMl6F07pB5J+MLXefXroMvtg6dVaaqFtBlF/3/9u48yraqvhP49yeCKBAMTrHFFqKoQURU0BYnUNvGiEK6NTgixviMaW1Xu0w7oU0M3XFYS0GNxqdNg9gKcUZBMYooy6BCmMEBBQ2IEcQWZ6b69R/3PPtSVtV7vKFunfc+n7fuqnvP2efsfW7BH/Vb3713zt754JWcHgMAAAAWMYtty3+eZIepz6cmefEwBWnNzlbbJflgkn2r6slrGlbVY6pqj0xSPD/s7rkkz80k3bIufe2Y5OqhmLN/knsNx89IcnBV3WHo+0+GY18ejt++qnZI8pQk6e6fJbm8qp4+jKuq6kHDvf42yVuq6g+Gc9tU1Z8vMLZvJblLVT1iaLd1VT2gu3+a5KdV9aih3bPX8dn+KZPpYWuuOWOR6wAAAGDU7HK1bgmdO1TVlVOf37qBfV6Q5OZhwd9jkxydyc5X59Rk7tE1SQ7u7uuq6sAkR1XVUUluHK59WSbr1Xy0JlunfzaLpHK6+9qq+kpNtv3+TJI3JflUVV2Y5Owk3xzanVNVxyb5+nDp+7r73CSpqhOTnJ/Juj1nTd3+2UneXVWHZzIF7IQk53f3KTXZWv3zw/N0kmMWGNsNwzSzt1fVjpn8Lo7KJOHz/CTHVFUn+dwi3+OnknxkWNj5pcPrf1fVXw3f4fMXuQ4AAAAYueo262ZLZcrVyvCZrbab9RCY8sozXzPrITA47cF/M+shMPj+1htrhjYb6s/Oe8Osh8DgzD1eOeshMOXN2/xq1kNg8Lqb1nVpUZbDw6/62Ga70sx7dn7Osv09+6IrP7Aiv8dZTLkCAAAAYAMo6AAAAACj0rV8r3VRVQdU1beq6jtV9aoFzr+8qi6pqguq6gtVda+F7nNrKOgAAAAArKeq2irJ3yV5UpLdkzyzqnaf1+zcJHt3955JPpLkzRvar4IOAAAAMCorbJerhyX5Tndf1t03ZLJp0kHTDbr7i929ZsGvrybZ+VY/9DwKOgAAAADr7x5Jrpj6fOVwbDEvyGQn7g1i2woAAACARVTVqiSrpg6t7u7V63mv5yTZO8ljN3RcCjoAAADAqKzjVKiNYijeLFXA+UGSe0593nk4dgtV9YQkr03y2O6+fkPHZcoVAAAAwPo7K8luVbVrVW2T5BlJTppuUFUPTvKeJE/t7qs3RqcSOgAAAMCo9KwHMKW7b6qqlyQ5NclWSY7p7our6g1Jzu7uk5K8Jcn2ST5cVUnyL9391A3pV0EHAAAAYAN09ylJTpl37PVT75+wsftU0AEAAABGZa5mPYLZs4YOAAAAwMhI6AAAAACjspy7XK1UEjoAAAAAIyOhAwAAAIyKhI6EDgAAAMDoSOgAAAAAo9KzHsAKIKEDAAAAMDISOgAAAMCozNWsRzB7EjoAAAAAIyOhAwAAAIyKXa4kdAAAAABGR0EHAAAAYGRMuQIAAABGxbblEjoAAAAAoyOhswW738E3zHoIJDn+09vOeghMOePBb5j1EBg8/rwjZj0EBq9+hP8vVooz93jlrIfA4BEXvWnWQ2DKFXs8Z9ZDYHDZbXad9RCY8vBZD2ATmpPRkdABAAAAGBsJHQAAAGBUbFsuoQMAAAAwOhI6AAAAwKhYQUdCBwAAAGB0JHQAAACAUbGGjoQOAAAAwOhI6AAAAACjMlezHsHsSegAAAAAjIyEDgAAADAqc/a5ktABAAAAGBsJHQAAAGBU5HMkdAAAAABGR0EHAAAAYGRMuQIAAABGZW7WA1gBJHQAAAAARkZCBwAAABgV25ZL6AAAAACMjoQOAAAAMCryORI6AAAAAKMjoQMAAACMil2uJHQAAAAARkdCBwAAABgVu1xJ6AAAAACMjoQOAAAAMCryORI6AAAAAKOzzgWdqvrFAsf+oqoOvTUdVtU/DT93qapnreM1962qU6rq0qo6p6r+oaruVlWHVdU7b03/G1NVfa+qLqyq84bX25ehz7V+51V1RFW9YlOPBQAAAGZhbhlfK9UGTbnq7r9fj2v2Hd7ukuRZST64VPuq2jbJyUle3t2fGo7tl+Qut7bvTWT/7v7xcnW2Pt85AAAAsHnZoClX00mQqjq9qt5WVWdX1Teqap+q+tiQqjly6po1SZ83Jnn0kGz5r1W1VVW9parOqqoLqupFQ7tnJTlzTTEnSbr79O6+aPj4b6rqs0M/b57q54lVdeaQ6PlwVW0/HD+gqr45HH97VX16/rMMny+qql2G98+pqq8PY31PVW21lu/lv1TVJcNznDB1/+OHMV1aVS+cav9XU8/911PHDx2OnV9Vxy/wnb9wuO78qvpoVd1hHX5tAAAAMGq9jP9Wqo29KPIN3b13Vb0sySeTPDTJT5J8t6re1t3XTrV9VZJXdPeBSVJVq5Jc1937VNXtknylqj6XZI8k/7xEn3sleXCS65N8q6rekeTXSQ5P8oTu/mVVvTLJy4eCz3uTPC7Jd5KcuLYHqqo/SnJIkkd2941V9a4kz07y/qHJF6vq5uH9cd39tuHZdu3u66vqjlO32zPJv0uyXZJzq+rk4fl2S/KwJJXkpKp6TJJrh2fYt7t/XFU7LTC8j3X3e4dxHpnkBUnesbZnAgAAAMZtYxd0Thp+Xpjk4u7+YZJU1WVJ7plJkWIxT0yyZ1U9bfi8YyaFjrX5QndfN/RzSZJ7Jbljkt0zKQolyTZJzkxy/ySXd/elQ/sPJFm1lvs/PpPC1FnDvW6f5Oqp8wtNubogyf+pqk8k+cTU8U9296+T/LqqvphJEedRw7OfO7TZfnjuByX58Jp7d/dPFhjbHkMh547Ddaeu5VnWFM5WJcnR+z0gz3/Av13bJQAAAMAKs7ELOtcPP+em3q/5vLa+KslLu/sWRYmqumeSx65Dn0ly89BPJfnH7n7mvHvttcR9bsotp6BtOzWu47r71UsP/xaenOQxSZ6S5LVV9cDh+PysVg/3/9vufs+8sb50Hfo5NsnB3X1+VR2WZL+1XdDdq5OsTpKfv+SPV252DAAAABaxkhcrXi6z3Lb850l2mPp8apIXV9XWyW93ttouk0WT962qJ69pWFWPqao9lrj3V5M8sqruM7Tfrqrum+SbSXapqnsP7aYLPt9L8pCh/UOS7Doc/0KSp1XVXYdzO1XVvRbruKpuk+Se3f3FJK/MJGm0/XD6oKratqrulEnx5azhuf9sao2fewx9nZbk6UPbLDLlaockPxy+s2cv8X0AAAAAm5Fbk9C5Q1VdOfX5rRvY9wVJbq6q8zNJmhydyc5X59RkbtM1maRPrquqA5McVVVHJblxuPZli924u68ZEisfGtbjSZLDu/vbw5Sjk6vqV0nOyP8vKn00yaFVdXGSryX59nCvS6rq8CSfG4o1Nyb5z0m+P1w3vYbOBZmsY/OBqtoxk/TN27v7p8N0rQuSfDHJnZP8TXdfleSqYZ2eM4c2v0jynO6+uKr+R5IvDfc/N8lh8x71dcNYrxl+7hAAAADYzM2t4MWKl8s6F3S6e8k0T3fvN/X+OpeHuwAADgtJREFU9CSnL3Ju++HnjZksTjztNcNr/r2/meSABbo9dnitaXfg1PvTkuyzwL0+m8laOmu2P3/FcPzXmaxls9CznZgFFlDu7l0Wap/JujgLuaC7D13gPkdnUtCaf/y4JMfNO3bE1Pt3J3n3AtcdMf8YAAAAsPnY2GvoAAAAAGxS8jlbeEFnfpJoE/d1xHL0AwAAAGz+tuiCDgAAADA+1tCZ7S5XAAAAAKwHCR0AAABgVOZmPYAVQEIHAAAAYGQkdAAAAIBRaWvoSOgAAAAAjI2EDgAAADAq1tCR0AEAAAAYHQkdAAAAYFSsoSOhAwAAADA6CjoAAAAAI2PKFQAAADAqFkWW0AEAAAAYHQkdAAAAYFTm2qLIEjoAAAAAIyOhAwAAAIyKfI6EDgAAAMDoSOgAAAAAozInoyOhAwAAADA2EjoAAADAqLSEjoQOAAAAwNhI6AAAAACjMjfrAawAEjoAAAAAIyOhswX7/dXnz3oIJDnqbvvPeghMOXLrH816CAx+sOdzZz0EBq/e5v6zHgKDN2/zs1kPgcEVezxn1kNgytkXfWDWQ2Dwm8P/ctZDYAux0na5qqoDkhydZKsk7+vuN847f7sk70/y0CTXJjmku7+3IX1K6AAAAACsp6raKsnfJXlSkt2TPLOqdp/X7AVJ/m933yfJ25K8aUP7VdABAAAARqWX8d86eFiS73T3Zd19Q5ITkhw0r81BSY4b3n8kyeOrqjbkO1DQAQAAAFh/90hyxdTnK4djC7bp7puSXJfkThvSqYIOAAAAwCKqalVVnT31WjXrMSUWRQYAAABGZjm3Le/u1UlWL9HkB0nuOfV55+HYQm2urKrbJtkxk8WR15uEDgAAAMD6OyvJblW1a1Vtk+QZSU6a1+akJM8b3j8tyWndvUFbdUnoAAAAAKOygbWQjaq7b6qqlyQ5NZNty4/p7our6g1Jzu7uk5L8ryTHV9V3kvwkk6LPBlHQAQAAANgA3X1KklPmHXv91PvfJHn6xuxTQQcAAAAYlbl12058s2YNHQAAAICRkdABAAAARmU5d7laqSR0AAAAAEZGQgcAAAAYlbaGjoQOAAAAwNhI6AAAAACjYpcrCR0AAACA0ZHQAQAAAEalW0JHQgcAAABgZCR0AAAAgFGZm/UAVgAJHQAAAICRkdABAAAARqXtciWhAwAAADA2CjoAAAAAI2PKFQAAADAqc6ZcSegAAAAAjI2EznqqqpuTXJhk6yQ3JXl/krd19ybbPa2qjk3y6e7+yKbqAwAAAFa6bgkdBZ319+vu3itJququST6Y5PeS/PeZjmpKVW3V3TfPehwAAADAxmXK1UbQ3VcnWZXkJTWxVVW9parOqqoLqupFSVJVJ1TVk9dcV1XHVtXTlmhfVfXOqvpWVX0+yV2nrn18VZ1bVRdW1TFVdbvh+Peq6k1VdU6Spy/n9wAAAADLYS69bK+VSkFnI+nuy5JslUnR5QVJruvufZLsk+SFVbVrkhOT/GmSVNU2SR6f5OQl2v9Jkvsl2T3JoUn2Ha7dNsmxSQ7p7gdmkrR68dRwru3uh3T3CZv0oQEAAICZUNDZNJ6Y5NCqOi/J15LcKcluST6TZP8hTfOkJF/u7l8v0f4xST7U3Td391VJThvuf78kl3f3t4fPxw1t1zhxsYFV1aqqOruqzp6b++VGelwAAABYPr2M/1Yqa+hsJFX1h0luTnJ1kkry0u4+dYF2pyf5D0kOSbImQbNg+6r64/UczqKVmu5enWR1ktx2m3us3P8yAQAAgEVJ6GwEVXWXJH+f5J09WWr71CQvrqqth/P3rarthuYnJnl+kkcn+exwbLH2X05yyLDGzt2T7D+0/1aSXarqPsPn5yb50iZ9SAAAAFgh5rqX7bVSSeisv9sPU6TWbFt+fJK3Dufel2SXJOdUVSW5JsnBw7nPDW0/2d03rKX9x5M8LsklSf4lyZlJ0t2/qarnJ/lwVd02yVmZFJQAAACALYCCznrq7q2WODeX5DXDa/65G5PstK7tk7xkkT6+kOTBCxzfZalxAwAAwNit3NzM8jHlCgAAAGBkJHQAAACAUZmT0ZHQAQAAABgbCR0AAABgVCR0JHQAAAAARkdBBwAAAGBkTLkCAAAARqXblCsJHQAAAICRkdABAAAARsWiyBI6AAAAAKMjoQMAAACMSkvoSOgAAAAAjI2EDgAAADAqdrmS0AEAAAAYHQkdAAAAYFTsciWhAwAAADA6EjoAAADAqFhDR0IHAAAAYHQkdAAAAIBRsYaOhA4AAADA6EjoAAAAAKPSEjoSOgAAAABjo6ADAAAAMDKmXAEAAACjMmfbcgkdAAAAgLGR0NmCPfquu896CERVdaV5zY07zXoIDL677d1mPQQGNTfrEbDG62663ayHwOCy2+w66yEw5TeH/+Wsh8Bg2yPfNeshsIWwKLK/JQEAAABGR0IHAAAAGBVr6EjoAAAAAIyOhA4AAAAwKtbQkdABAAAAGB0JHQAAAGBUrKEjoQMAAAAwOhI6AAAAwKhYQ0dCBwAAAGB0JHQAAACAUbGGjoQOAAAAwOhI6AAAAACjYg0dCR0AAACA0VHQAQAAABgZU64AAACAUemem/UQZk5CBwAAAGBkJHQAAACAUZmzKLKEDgAAAMDYSOgAAAAAo9ItoSOhAwAAADAyEjoAAADAqFhDR0IHAAAAYHQkdAAAAIBRsYaOhA4AAADA6CjoAAAAAKMy171srw1VVTtV1T9W1aXDz99foM1eVXVmVV1cVRdU1SFru6+CDgAAAMCm86okX+ju3ZJ8Yfg836+SHNrdD0hyQJKjquqOS910syroVNVrp6pZ51XVw5doe2xVPW2B4/tV1afX0s9hVXVNVZ07VNhOrap9N8YzrKXfX2zqPgAAAGCl62X8txEclOS44f1xSQ7+nefp/nZ3Xzq8vyrJ1UnustRNN5tFkavqEUkOTPKQ7r6+qu6cZJtN2OWJ3f2Soe/9k3ysqvbv7m9swj7XWVXdtrtvmvU4AAAAYAt3t+7+4fD+X5PcbanGVfWwTOoZ312q3eaU0Ll7kh939/VJ0t0/7u6rqur1VXVWVV1UVaurquZfWFUHVNU3q+qcJP9x6vhOVfWJIfHz1arac6GOu/uLSVYnWTVcd++q+mxV/XNVnVFV96+qHavq+1V1m6HNdlV1RVVtvVD7oc2uwxy6C6vqyKlxVVW9ZXimC9fMrRvSRWdU1UlJLtlI3ysAAACsKN29bK+qWlVVZ0+9Vs0fT1V9fvgbff7roHnj7mTx2E9V3T3J8Ume391zS30Hm1NB53NJ7llV366qd1XVY4fj7+zufbp7jyS3zyTF81tVtW2S9yZ5SpKHJvmDqdN/neTc7t4zyWuSvH+J/s9Jcv/h/eokL+3uhyZ5RZJ3dfd1Sc5LsmZcByY5tbtvXKj90OboJO/u7gcmWVPNSyZFp72SPCjJE5K8ZfilJ8lDkrysu++7xFgBAACAddDdq7t776nX6gXaPKG791jg9ckkP1rzN/vw8+qF+qmq30tycpLXdvdX1zauzaag092/yKQgsyrJNUlOrKrDkuxfVV+rqguTPC7JA+Zdev8kl3f3pUOl7ANT5x6VSWUs3X1akjsNX/BCKkmqavsk+yb5cFWdl+Q9maSHkuTEJGtWqn7GMMal2j8yyYeG98fPG9eHuvvm7v5Rki8l2Wc49/XuvnyRMd6isnjVL69crBkAAACwcZyU5HnD++cl+eT8BlW1TZKPJ3l/d39kXW662ayhkyTdfXOS05OcPhRwXpRkzyR7d/cVVXVEkm03UfcPTvKNTIpkP+3uvRZoc1KS/1lVO2VSfDotyXZLtE+WiGIt4pdLnRwqiauTZP+d//1GWd0JAAAAltPcxlmseLm8Mck/VNULknw/yZ8mSVXtneQvuvvPh2OPySRIcthw3WHdfd5iN91sEjpVdb+q2m3q0F5JvjW8//GQhPmdXa2SfDPJLlV17+HzM6fOnZHk2cP998tkjZ6fLdD3YzNJBr13OH95VT19OFdV9aDktymiszKZSvXpIWGzaPskX8kkyZM145ga1yFVtVVV3SWTX/rXl/h6AAAAgBno7mu7+/HdvdswNesnw/Gzh2JOuvsD3b11d+819Vq0mJNsXgmd7ZO8oyb7tN+U5DuZFFl+muSiTFaSPmv+Rd39m2FBo5Or6leZFEt2GE4fkeSYqrogkz3hnzd16SFV9agkd0hyeZL/NLXD1bOTvLuqDk+ydZITkpw/nDsxyYeT7Dd1r8XavyzJB6vqlbllJOvjSR4xtOkk/627/3XNYsoAAACwOZusmLJlK1/ClsuUq5Xh6bXkjnUss91uuHHWQ2Dw3W22nvUQGGy75P4KLKc/6l/NeggMLsvtZz0Ephz41GtnPQQG2x75rrU3Ytlsfec//J1dnjcXd/69+y7b37M//tm3V+T3uDkldAAAAIAtwJxwyuazhg4AAADAlkJCBwAAABgVy8dI6AAAAACMjoQOAAAAMCpzkdCR0AEAAAAYGQkdAAAAYFSsoSOhAwAAADA6EjoAAADAqMxJ6EjoAAAAAIyNhA4AAAAwKm2XKwkdAAAAgLFR0AEAAAAYGVOuAAAAgFGxKLKEDgAAAMDoSOgAAAAAo9ISOhI6AAAAAGMjoQMAAACMim3LJXQAAAAARkdCBwAAABgVa+hI6AAAAACMjoQOAAAAMCoSOhI6AAAAAKMjoQMAAACMinyOhA4AAADA6JR5ZwAAAADjIqEDAAAAMDIKOgAAAAAjo6ADAAAAMDIKOgAAAAAjo6ADAAAAMDIKOgAAAAAj8/8AVAx5HLNc734AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import seaborn as sns\n", "\n", "corrmat = df_amostra.corr()\n", "corrmat\n", "\n", "fig, ax = plt.subplots(figsize=(20,10)) \n", "sns.heatmap(corrmat, vmax=1., square=False).xaxis.tick_top()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
IdadeCartaodeCreditoChequeEspecialRendaLimiteCartaodeCreditoLimiteChequeEspecialDevedor
Idade9.8559140.0903230.160215907.6309891.573383e+031.360068e+030.021505
CartaodeCredito0.0903230.161290-0.02580642.7032907.128339e+02-3.017626e+010.019355
ChequeEspecial0.160215-0.0258060.1161299.473860-8.959438e+013.374318e+02-0.020430
Renda907.63098942.7032909.473860211575.4193974.520597e+052.349786e+0517.880688
LimiteCartaodeCredito1573.382817712.833935-89.594376452059.7499223.545414e+061.940102e+05104.674366
LimiteChequeEspecial1360.068269-30.176258337.431828234978.6268751.940102e+051.187914e+06-44.347871
Devedor0.0215050.019355-0.02043017.8806881.046744e+02-4.434787e+010.090323
\n", "
" ], "text/plain": [ " Idade CartaodeCredito ChequeEspecial \\\n", "Idade 9.855914 0.090323 0.160215 \n", "CartaodeCredito 0.090323 0.161290 -0.025806 \n", "ChequeEspecial 0.160215 -0.025806 0.116129 \n", "Renda 907.630989 42.703290 9.473860 \n", "LimiteCartaodeCredito 1573.382817 712.833935 -89.594376 \n", "LimiteChequeEspecial 1360.068269 -30.176258 337.431828 \n", "Devedor 0.021505 0.019355 -0.020430 \n", "\n", " Renda LimiteCartaodeCredito \\\n", "Idade 907.630989 1.573383e+03 \n", "CartaodeCredito 42.703290 7.128339e+02 \n", "ChequeEspecial 9.473860 -8.959438e+01 \n", "Renda 211575.419397 4.520597e+05 \n", "LimiteCartaodeCredito 452059.749922 3.545414e+06 \n", "LimiteChequeEspecial 234978.626875 1.940102e+05 \n", "Devedor 17.880688 1.046744e+02 \n", "\n", " LimiteChequeEspecial Devedor \n", "Idade 1.360068e+03 0.021505 \n", "CartaodeCredito -3.017626e+01 0.019355 \n", "ChequeEspecial 3.374318e+02 -0.020430 \n", "Renda 2.349786e+05 17.880688 \n", "LimiteCartaodeCredito 1.940102e+05 104.674366 \n", "LimiteChequeEspecial 1.187914e+06 -44.347871 \n", "Devedor -4.434787e+01 0.090323 " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Dados mulheres\n", "\n", "S1 = df_amostra.iloc[:,1:8][df_amostra['Sexo']==1].cov()\n", "n1 = len(df_amostra[df_amostra['Sexo']==1])\n", "Xbarra1 = df_amostra[df_amostra['Sexo']==1].mean()\n", "S1" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "31" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n1" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
IdadeCartaodeCreditoChequeEspecialRendaLimiteCartaodeCreditoLimiteChequeEspecialDevedor
Idade15.503126-0.0247160.0285041530.5001557.437871e+023.557481e+020.041430
CartaodeCredito-0.0247160.2386680.003039-6.5271321.053693e+037.902214e+000.001473
ChequeEspecial0.0285040.0030390.1502182.2603851.498600e+014.464707e+020.004352
Renda1530.500155-6.5271322.260385240646.3496431.154767e+055.165076e+042.773754
LimiteCartaodeCredito743.7870711053.69313814.986003115476.6597274.868377e+067.295678e+047.124906
LimiteChequeEspecial355.7481027.902214446.47073551650.7584007.295678e+041.371908e+0616.766938
Devedor0.0414300.0014730.0043522.7737547.124906e+001.676694e+010.193529
\n", "
" ], "text/plain": [ " Idade CartaodeCredito ChequeEspecial \\\n", "Idade 15.503126 -0.024716 0.028504 \n", "CartaodeCredito -0.024716 0.238668 0.003039 \n", "ChequeEspecial 0.028504 0.003039 0.150218 \n", "Renda 1530.500155 -6.527132 2.260385 \n", "LimiteCartaodeCredito 743.787071 1053.693138 14.986003 \n", "LimiteChequeEspecial 355.748102 7.902214 446.470735 \n", "Devedor 0.041430 0.001473 0.004352 \n", "\n", " Renda LimiteCartaodeCredito \\\n", "Idade 1530.500155 7.437871e+02 \n", "CartaodeCredito -6.527132 1.053693e+03 \n", "ChequeEspecial 2.260385 1.498600e+01 \n", "Renda 240646.349643 1.154767e+05 \n", "LimiteCartaodeCredito 115476.659727 4.868377e+06 \n", "LimiteChequeEspecial 51650.758400 7.295678e+04 \n", "Devedor 2.773754 7.124906e+00 \n", "\n", " LimiteChequeEspecial Devedor \n", "Idade 3.557481e+02 0.041430 \n", "CartaodeCredito 7.902214e+00 0.001473 \n", "ChequeEspecial 4.464707e+02 0.004352 \n", "Renda 5.165076e+04 2.773754 \n", "LimiteCartaodeCredito 7.295678e+04 7.124906 \n", "LimiteChequeEspecial 1.371908e+06 16.766938 \n", "Devedor 1.676694e+01 0.193529 " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Dados homens\n", "S2 = df.iloc[:,1:8][df['Sexo']==0].cov()\n", "n2 = len(df_amostra[df_amostra['Sexo']==0])\n", "Xbarra2 = df_amostra[df_amostra['Sexo']==0].mean()\n", "S2" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "69" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n2" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
IdadeCartaodeCreditoChequeEspecialRendaLimiteCartaodeCreditoLimiteChequeEspecialDevedor
Idade13.7743870.0105000.0688241339.8259219.977450e+026.631931e+020.035331
CartaodeCredito0.0105000.214981-0.0057918.5434059.493485e+02-3.754461e+000.006947
ChequeEspecial0.068824-0.0057910.1397834.468591-1.702840e+014.130915e+02-0.003235
Renda1339.8259218.5434054.468591231747.0852822.185123e+051.077715e+057.398325
LimiteCartaodeCredito997.744952949.348484-17.028399218512.2995834.463388e+061.100140e+0536.986985
LimiteChequeEspecial663.193051-3.754461413.091478107771.5344641.100140e+051.315583e+06-1.941677
Devedor0.0353310.006947-0.0032357.3983253.698699e+01-1.941677e+000.161935
\n", "
" ], "text/plain": [ " Idade CartaodeCredito ChequeEspecial \\\n", "Idade 13.774387 0.010500 0.068824 \n", "CartaodeCredito 0.010500 0.214981 -0.005791 \n", "ChequeEspecial 0.068824 -0.005791 0.139783 \n", "Renda 1339.825921 8.543405 4.468591 \n", "LimiteCartaodeCredito 997.744952 949.348484 -17.028399 \n", "LimiteChequeEspecial 663.193051 -3.754461 413.091478 \n", "Devedor 0.035331 0.006947 -0.003235 \n", "\n", " Renda LimiteCartaodeCredito \\\n", "Idade 1339.825921 9.977450e+02 \n", "CartaodeCredito 8.543405 9.493485e+02 \n", "ChequeEspecial 4.468591 -1.702840e+01 \n", "Renda 231747.085282 2.185123e+05 \n", "LimiteCartaodeCredito 218512.299583 4.463388e+06 \n", "LimiteChequeEspecial 107771.534464 1.100140e+05 \n", "Devedor 7.398325 3.698699e+01 \n", "\n", " LimiteChequeEspecial Devedor \n", "Idade 6.631931e+02 0.035331 \n", "CartaodeCredito -3.754461e+00 0.006947 \n", "ChequeEspecial 4.130915e+02 -0.003235 \n", "Renda 1.077715e+05 7.398325 \n", "LimiteCartaodeCredito 1.100140e+05 36.986985 \n", "LimiteChequeEspecial 1.315583e+06 -1.941677 \n", "Devedor -1.941677e+00 0.161935 " ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S_pooled = ((n1-1)*S1 + (n2-1)*S2)/(n1+n2-2)\n", "S_pooled" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def T2Hotelling_duas_amostras(df1, df2, delta0):\n", " n1 = len(df1)\n", " n2 = len(df2)\n", " p = len(df1.columns)\n", " Xbarra1=df1.mean()\n", " Xbarra2=df2.mean()\n", " S1 = df1.cov()\n", " S2 = df2.cov()\n", " S_pooled = ((n1-1)*S1 + (n2-1)*S2)/(n1+n2-2)\n", " S_pooled_inv = np.linalg.inv(S_pooled)\n", " \n", " T2Hotelling_duas_amostras = np.array(Xbarra1-Xbarra2-delta0).T.dot(S_pooled_inv).dot(np.array(Xbarra1-Xbarra2-delta0)) / (n1+n2-2)\n", " qf = f.ppf(0.95, p , (n1+n2-2), loc=0, scale=1)\n", " teste = T2Hotelling_duas_amostras > (n1+n2-2) * p / (n1+n2-p-1) * qf\n", " pvalor = 1-f.cdf(T2Hotelling_duas_amostras/((n1+n2-2) * p / (n1+n2-p-1) ), p, (n1+n2-2))\n", " print('Rejeitamos H0') if teste else print('Não rejeitamos H0')\n", " print('Valor da estatística', T2Hotelling_duas_amostras)\n", " print('valor p', pvalor)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Não rejeitamos H0\n", "Valor da estatística 0.05027388876482565\n", "valor p 0.9999998142670615\n" ] } ], "source": [ "df1 = df_amostra.iloc[:,1:8][df_amostra['Sexo']==1]\n", "df2 = df_amostra.iloc[:,1:8][df_amostra['Sexo']==0]\n", "\n", "delta0 = [0,0,0,0,0,0,0]\n", "\n", "T2Hotelling_duas_amostras(df1,df2,delta0)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testes de hipóteses para a comparação de médias em amostras correlacionadas\n", "\n", "Sejam \n", "\n", "\n", "- $\\underline{X}_{11},\\ldots,\\underline{X}_{1n}$ vetores $p\\times 1$ que representam uma amostra aleatória de uma população normal multivariada **antes** de um tratamento com $E(\\underline{X}_{1j})=\\underline{\\mu}_1$ para $j=1,\\ldots,n$,\n", "\n", "\n", "- $\\underline{X}_{21},\\ldots,\\underline{X}_{2n}$ vetores $p\\times 1$ que representam uma amostra aleatória de uma população normal multivariada **após** de um tratamento com $E(\\underline{X}_{2j})=\\underline{\\mu}_2$ para $j=1,\\ldots,n$,\n", "\t\t\t \n", "sendo que $\\underline{X}_{11},\\ldots,\\underline{X}_{1n}$ e $\\underline{X}_{21},\\ldots,\\underline{X}_{2n}$ são amostras aleatórias de uma mesma população em diferentes situações, em que $\\underline{X}_{1j}$ e $\\underline{X}_{2j}$ são correlacionadas (por exemplo, vetores aleatórios de medições antes e após um tratamento).\n", "\n", "\n", "Sejam $\\underline{\\mu}_1$ e $\\underline{\\mu}_2$ os vetores de médias em situações 1 e 2, respectivamente. Essas situações podem não ser, necessariamente, antes e após um tratamento, mas são situções que indicam que as amostras estão correlacionadas. \n", "\n", "\n", "Deseja-se testar se não há diferença entre as situações 1 e 2 para verificar, por exemplo, que o tratamento não produz nenhum efeito, ou seja, se $\\underline{\\mu}_1 = \\underline{\\mu}_2$, ou equivalentemente se as médias são iguais em outra situação.\n", "\t\t\t \n", " \n", "Para avaliar as hipóteses \n", "\n", "\n", "$\\begin{array}{l}H_0:\\underline{\\mu}_1=\\underline{\\mu}_2\\mbox{ contra }\\\\H_1:\\underline{\\mu}_1\\neq\\underline{\\mu}_2\\end{array}$\n", "\t\t\t \n", "vamos considerar as diferenças\n", "\t\t\t\n", "$\\underline{D}_j = \\underline{X}_{1j}-\\underline{X}_{2j}.$\n", "\t\t\t\n", "Assim, $\\underline{D}_1,\\ldots,\\underline{D}_n$ são i.i.d e $\\underline{D}_j\\sim N(\\underline{\\mu}_D, \\Sigma_D)$.\n", "\n", "Então, avaliamos se \n", "\t\t\t\n", "\t\t\t\n", "$\\begin{array}{l}H_0:\\underline{\\mu}_D=\\underline{0}\\mbox{ contra }\\\\H_1:\\underline{\\mu}_D\\neq\\underline{0}\\end{array}$\n", "\t\t\t\n", "com a estatística $T^2$ de Hotelling:\n", "\t\t\t \n", "$T^2 = n(\\bar{\\underline{D}}-\\underline{0})^\\top S_D^{-1} (\\bar{\\underline{D}}-\\underline{0})\\stackrel{sob \\ H_0}{\\sim} \\displaystyle\\frac{(n-1)p}{n-p}F_{p,n-p},$\n", "\n", "em que $\\bar{\\underline{D}}$ e $S_D$ são o vetor de médias e a matriz de variâncias e covariâncias amostrais de $\\underline{D}$.\n", "\n", "\t\t\t\n", "\n", "Um teste análogo poderia ser desenvolvido para avaliar\n", "\t\t\t\n", "$\\begin{array}{l}H_0:\\underline{\\mu}_D=\\underline{\\mu}_{D0}\\mbox{ contra }\\\\H_1:\\underline{\\mu}_D\\neq\\underline{\\mu}_{D0}\\end{array}$\n", "\t\t\t \n", "com a estatística $T^2$ de Hotelling:\n", "\n", "$T^2 = n(\\bar{\\underline{D}}-\\underline{\\mu}_{D0})^\\top S_D^{-1} (\\bar{\\underline{D}}-\\underline{\\mu}_{D0})\\stackrel{ sob \\ H_0}{\\sim} \\displaystyle\\frac{(n-1)p}{n-p}F_{p,n-p},$\n", "\n", "\n", "em que $\\bar{\\underline{D}}$ e $S_D$ são o vetor de médias e a matriz de variâncias e covariâncias amostrais de $\\underline{D}$.\n", "\n", "A região de confiança, com nível de confiança $100(1-\\alpha)\\%$ nesse caso seria \n", "\n", "$\\{\\underline{\\mu}_D^\\star; n(\\bar{\\underline{D}}-\\underline{\\mu}_{D}^\\star)^\\top S_D^{-1} (\\bar{\\underline{D}}-\\underline{\\mu}_{D}^\\star) \\leq \\displaystyle\\frac{(n-1)p}{n-p}q_{F_{p,n-p,\\alpha} }\\}$\n", "\n", "\n", "em que $\\bar{\\underline{D}}$ e $S_D$ são o vetor de médias e a matriz de variâncias e covariâncias amostrais de $\\underline{D}$.\n", "\t\t\t \n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulação: Amostras pareadas" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "mean = [0, 0, 0]\n", "mean1 = [1, 0, 0]\n", "\n", "cov1 = [[2,1,0],[1,3,1],[0,1,4]] \n", "cov2 = [[0.01,0,0],[0,0.01,0],[0,0,0.01]]\n", "\n", "X1 = np.random.multivariate_normal(mean, cov1, 50)\n", "X2 = X1 + np.random.multivariate_normal(mean1, cov2, 50)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(50, 6)" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X = np.concatenate((X1,X2), axis=1)\n", "X.shape" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
012345
00.2618770.0464560.2616201.4162200.1679720.132674
10.4645192.3945490.8853941.5264912.3393890.870196
22.0216253.896733-1.1909462.9357793.718885-1.122018
3-2.4522830.3948801.681235-1.4477420.2003661.630899
40.1495040.6146911.6869461.1733630.6177071.710605
\n", "
" ], "text/plain": [ " 0 1 2 3 4 5\n", "0 0.261877 0.046456 0.261620 1.416220 0.167972 0.132674\n", "1 0.464519 2.394549 0.885394 1.526491 2.339389 0.870196\n", "2 2.021625 3.896733 -1.190946 2.935779 3.718885 -1.122018\n", "3 -2.452283 0.394880 1.681235 -1.447742 0.200366 1.630899\n", "4 0.149504 0.614691 1.686946 1.173363 0.617707 1.710605" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.DataFrame(X)\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 0.090345\n", "1 -0.025920\n", "2 0.425570\n", "3 1.103482\n", "4 -0.021276\n", "5 0.415850\n", "dtype: float64" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.mean()" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD4CAYAAABPLjVeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAR8UlEQVR4nO3df7Bkd1nn8fdnbhLBCSZWoRhmoknpsG5ECXhr1IpGVNAJUkTLXwkVMW7MtUqiuKgQ1Aoaf6PgrmVEBzZFZReTQl13p2Q0UDqEwhUyNxBiMhidjayZWXQIIhhQk3vv4x+3g811bv+Y6T6nz5n3a+rUdJ8+/e2nZrqe+9znfM/3pKqQJDVjR9sBSNKZxKQrSQ0y6UpSg0y6ktQgk64kNcikK0kN6mTSTbIvyYNJjia5se14xklya5ITSe5vO5ZJJLkwyaEkR5I8kORlbcc0TpInJbk7yfsHMf902zFNIslSkvcl+YO2Y5lEkg8m+fMk9yZZbTueLkrX5ukmWQL+Eng+cAw4DFxdVUdaDWyEJJcDjwK3VdUz245nnCQXABdU1XuTPAW4B/iWBf83DrCzqh5NcjbwLuBlVfXulkMbKcnLgWXgs6rqhW3HM06SDwLLVfVI27F0VRcr3b3A0ap6qKoeA+4Armw5ppGq6p3A37cdx6Sq6kNV9d7B438EPgDsajeq0WrTo4OnZw+2ha4okuwGvhl4Y9uxqDldTLq7gIeHnh9jwRNClyW5CHg28J52Ixlv8Kv6vcAJ4O1Vtegx/xfgFcBG24FMoYC3JbknyUrbwXRRF5OuGpLkXOD3gB+uqo+3Hc84VbVeVZcCu4G9SRa2lZPkhcCJqrqn7Vim9NVV9RzgCuClg9aZptDFpHscuHDo+e7BPs3QoC/6e8Cbq+p/th3PNKrqH4BDwL62YxnhMuBFgx7pHcDXJ/kf7YY0XlUdH/x9Avh9Ntt9mkIXk+5hYE+Si5OcA1wFHGg5pl4ZnJT6b8AHqup1bccziSSfk+T8weMns3mi9S/ajWp7VfWqqtpdVRex+R3+k6q6puWwRkqyc3BilSQ7gW8EOjEjZ5F0LulW1RpwA3Anmyd43lJVD7Qb1WhJbgf+DPgPSY4lua7tmMa4DPhuNquvewfbC9oOaowLgENJ7mPzB/Pbq6oT07A65GnAu5K8H7gbeGtV/VHLMXVO56aMSVKXda7SlaQuM+lKUoNMupLUIJOuJDWo00m3a1fEdC1e6F7MXYsXjHmRjVusKpt+bbD41n1JnjNuzE4nXaBr//Fdixe6F3PX4gVjXmRvYvRFNlcAewbbCvD6cQN2PelK0txMsFjVlWyuHliDFe3OH6zSt62zZhngyTz+yENzmwj8G6/92ZmPf/uzbprlcJ/m2vP2ctuua2b+73H9R+6a9ZCfsnTW+XzGky6cecyPHptPzPP4TgC86dL5fS+uOX8vb9g9++/FD5w4NOshP2XH0nmcfc6umcf8+GPHc9pjTPH/f87nfOH38+lV+/6q2j/Fx223ANeHtnvD3JPuPH3fS65uO4SpfN3OPW2HMLWlpXPbDmEqXftOAFzewe/Fjh072w5hJgYJdpoke9o6nXQl6d/ZWG/y06ZegMuerqR+WV+bfDt9B4CXDGYxfCXwsaratrUAVrqSeqZqdmvCDxarei7w1CTHgFezeVcSquo3gYPAC4CjwCeB7x03pklXUr9szC7pVtXIkwS1uWLYS6cZ06QrqV9mWOnOg0lXUr80eyJtaiZdSf1ipStJzanZzEqYG5OupH6Z4Ym0eTDpSuoX2wuS1CBPpElSg6x0JalBnkiTpAZ5Ik2SmlPV8Z5uki9mc3X0XYNdx4EDVfWBeQYmSadkwXu6I5d2TPJK4A4gwN2DLcDtSW4c8b6VJKtJVt942+2zjFeSRtvYmHxrwbhK9zrgS6rq8eGdSV4HPAD84sneNLwa+zxv1yNJ/86CV7rjku4G8HTg/23Zf8HgNUlaLOuPjz+mReOS7g8Df5zkr/i3m699PvBFwA3zDEySTkmXZy9U1R8leQawl08/kXa4Fv0UoaQzU8fbC9TmvS/e3UAsknT6ulzpSlLnmHQlqTnV8RNpktQtXe/pSlKn2F6QpAZZ6UpSg6x0JalBVrqS1KA1FzGXpOZY6UpSg+zpSlKDrHQlqUFneqV7+7NumvdHzNTV77+57RCmdvfyq9oOYWpv7tj34toOfi/uX/7xtkNoh5WuJDXI2QuS1KBa7DuEmXQl9cuZ3tOVpEaZdCWpQZ5Ik6QGrS/27RtNupL6xfaCJDVowZPujrYDkKSZqo3JtzGS7EvyYJKjSW48yeufn+RQkvcluS/JC8aNadKV1Cu1URNvoyRZAm4BrgAuAa5OcsmWw34SeEtVPRu4CviNcfHZXpDUL7NrL+wFjlbVQwBJ7gCuBI4MHVPAZw0enwf8/3GDmnQl9csUsxeSrAArQ7v2V9X+weNdwMNDrx0DvmLLED8FvC3JDwI7geeN+0yTrqR+maLSHSTY/WMP3N7VwJuq6rVJvgr470meWbV9w9ikK6lfZtdeOA5cOPR892DfsOuAfQBV9WdJngQ8FTix3aCeSJPUL1WTb6MdBvYkuTjJOWyeKDuw5Zi/Ab4BIMl/BJ4EfHjUoFa6kvplRpVuVa0luQG4E1gCbq2qB5LcDKxW1QHgR4A3JPnPbJ5Uu7ZqdDY36UrqlzFTwaZRVQeBg1v23TT0+Ahw2TRjnnJ7Icn3nup7JWlu1tcn31pwOj3dn97uhSQrSVaTrB76xF+dxkdI0nRqY2PirQ0j2wtJ7tvuJeBp271veBrGbbuuWexl3CX1ywzbC/Mwrqf7NOCbgI9u2R/g/8wlIkk6HR1fT/cPgHOr6t6tLyR5x1wikqTT0eVKt6quG/Hai2cfjiSdpjUXMZek5nS8vSBJ3dLl9oIkdU1bU8EmZdKV1C9WupLUIJOuJDXIW7BLUnPG3fusbSZdSf1i0pWkBjl7QZIaZKUrSQ0y6UpSc2r9DG8vXP+Ru+b9ETN19/Kr2g5har+6+gtthzC1nbsubzuEqbxr+RVthzC116++pu0Q2mGlK0nNccqYJDXJpCtJDVrslq5JV1K/1NpiZ12TrqR+Weyca9KV1C+eSJOkJlnpSlJzrHQlqUlWupLUnFprO4LRTLqSemXB78Bu0pXUMyZdSWqOla4kNcikK0kNqvW0HcJIJl1JvbLole6OcQck+eIk35Dk3C37980vLEk6NbWRibc2jEy6SX4I+N/ADwL3J7ly6OWfn2dgknQqamPyrQ3j2gvXA19eVY8muQj43SQXVdV/Bbb9MZFkBVgBWDrrfJaWzt3uUEmaqarF7umOay/sqKpHAarqg8BzgSuSvI4RSbeq9lfVclUtm3AlNWmWlW6SfUkeTHI0yY3bHPOdSY4keSDJb48bc1yl+3dJLq2qewEGFe8LgVuBLx0fsiQ1a2NGsxeSLAG3AM8HjgGHkxyoqiNDx+wBXgVcVlUfTfK548YdV+m+BPjb4R1VtVZVLwG6dTtXSWeEGZ5I2wscraqHquox4A7gyi3HXA/cUlUfBaiqE+MGHZl0q+pYVf3tNq/96bjBJalp0yTdJCtJVoe2laGhdgEPDz0/Ntg37BnAM5L8aZJ3TzKry3m6knqlplhOt6r2A/tP4+POAvaweb5rN/DOJF9aVf8w6g2S1BsznH97HLhw6Pnuwb5hx4D3VNXjwF8n+Us2k/Dh7QYde3GEJHVJVSbexjgM7ElycZJzgKuAA1uO+V9sVrkkeSqb7YaHRg1qpSupV9ZnNHuhqtaS3ADcCSwBt1bVA0luBlar6sDgtW9McgRYB36sqj4yalyTrqRemeXFEVV1EDi4Zd9NQ48LePlgm4hJV1KvtLWmwqRMupJ6ZZrZC20w6UrqFStdSWrQ+sZiT8oy6UrqFdsLktSgjQVf2tGkK6lXFn09XZOupF4549sLjx67a94fMVNvftZN4w9aMDt3dW+VzU8cf2fbIUzl1ku797148tO/pu0Qprb22NalDaZne0GSGuTsBUlq0IJ3F0y6kvrF9oIkNcjZC5LUoAlu8tsqk66kXimsdCWpMWu2FySpOVa6ktQge7qS1CArXUlqkJWuJDVo3UpXkpqz4HfrMelK6pcNK11Jak7nF7xJsheoqjqc5BJgH/AXVXVw7tFJ0pQ6fSItyauBK4Czkrwd+ArgEHBjkmdX1c81EKMkTWwji91eGLfa77cDlwGXAy8FvqWqfgb4JuC7tntTkpUkq0lW33jb7TMLVpLGWZ9ia8O49sJaVa0Dn0zyf6vq4wBV9U9Jtq3iq2o/sB/g8UceWvQWi6Qe6frshceSfGZVfRL48id2JjmPxW+dSDoDdX32wuVV9S8AVTWcZM8GvmduUUnSKVr0X61HJt0nEu5J9j8CPDKXiCTpNHS9vSBJnbLofU+TrqReWbfSlaTmWOlKUoNMupLUoAW/RZpJV1K/WOlKUoPaurx3UuPWXpCkTtnI5Ns4SfYleTDJ0SQ3jjju25JUkuVxY5p0JfXKxhTbKEmWgFvYXGnxEuDqwfK2W497CvAy4D2TxGfSldQrs0q6wF7gaFU9VFWPAXcAV57kuJ8Bfgn450niM+lK6pWaYhtehnawrQwNtQt4eOj5scG+T0nyHODCqnrrpPF5Ik1Sr0yz9sLwMrTTSrIDeB1w7TTvM+lK6pUZzl44Dlw49Hz3YN8TngI8E3hHNu9W8XnAgSQvqqrV7Qade9J906U3zfsjZura99/cdghTe9fyK9oOYWq3dux78Z/u7d734r3Lr2w7hFZszG5xx8PAniQXs5lsrwJe/MSLVfUx4KlPPE/yDuBHRyVcsKcrqWdmdSKtqtaAG4A7gQ8Ab6mqB5LcnORFpxqf7QVJvTLLRcwHdz0/uGXfSX9Nq6rnTjKmSVdSr3gZsCQ1aC2LfcMek66kXlnslGvSldQzthckqUEznDI2FyZdSb2y2CnXpCupZ2wvSFKD1he81jXpSuoVK11JalBZ6UpSc6x0JalBThmTpAYtdso16UrqmbUFT7tTr6eb5LZ5BCJJs1BT/GnDyEo3yYGtu4CvS3I+QFWddCHfwc3dVgCuOX8vl+/cM4NQJWm8rp9I2w0cAd7I4OaZwDLw2lFvGr7Z2xt2X7PYtb6kXln0KWPj2gvLwD3ATwAfq6p3AP9UVXdV1V3zDk6SpjWr2/XMy8hKt6o2gF9N8juDv/9u3HskqU3rtdiV7kQJtKqOAd+R5JuBj883JEk6db2ap1tVbwXeOqdYJOm0LXpP11aBpF7p+uwFSeqUXrUXJGnR2V6QpAb1YvaCJHWF7QVJapAn0iSpQfZ0JalBthckqUHliTRJao63YJekBtlekKQGnfHthR84cWjeHzFT9y//eNshTO31q69pO4SpPfnpX9N2CFN57/Ir2w5har+++ktth9AKK11JapBTxiSpQV4GLEkNsr0gSQ0y6UpSgxZ99sK4uwFLUqdsUBNv4yTZl+TBJEeT3HiS11+e5EiS+5L8cZIvGDemSVdSr9QUf0ZJsgTcAlwBXAJcneSSLYe9D1iuqi8DfhcYO3/TpCupV9ZrY+JtjL3A0ap6qKoeA+4Arhw+oKoOVdUnB0/fDeweN6hJV1KvVNXEW5KVJKtD28rQULuAh4eeHxvs2851wB+Oi88TaZJ6ZZrZC1W1H9h/up+Z5BpgGfjacceadCX1ygyvSDsOXDj0fPdg36dJ8jzgJ4Cvrap/GTeoSVdSr2zMbsrYYWBPkovZTLZXAS8ePiDJs4HfAvZV1YlJBjXpSuqVWVW6VbWW5AbgTmAJuLWqHkhyM7BaVQeAXwbOBX4nCcDfVNWLRo1r0pXUKxPMSphYVR0EDm7Zd9PQ4+dNO6ZJV1KvzLC9MBcmXUm90qulHZN8NZsThu+vqrfNJyRJOnWLXumOvDgiyd1Dj68Hfh14CvDqk12HLEltm9VlwPMyrtI9e+jxCvD8qvpwkl9h85K3XzzZmwZXdawA7Fg6jx07ds4iVkkaa73W2w5hpHFJd0eSz2azIk5VfRigqj6RZG27Nw1f5XH2ObsWu9aX1CuLvrTjuKR7HnAPEKCSXFBVH0py7mCfJC2UTi9iXlUXbfPSBvCtM49Gkk5T1yvdkxosZfbXM45Fkk7bos9ecJ6upF7p1TxdSVp0s7wMeB5MupJ6pZc9XUlaVPZ0JalBVrqS1KBOz9OVpK6x0pWkBjl7QZIa5Ik0SWqQ7QVJapBXpElSg6x0JalBi97TzaL/VJCkPhl5jzRJ0myZdCWpQSZdSWqQSVeSGmTSlaQGmXQlqUH/CkR7k0KAF1YKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "corrmat = df.corr()\n", "corrmat\n", "\n", "sns.heatmap(corrmat, vmax=1., square=False).xaxis.tick_top()" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "df_diff = pd.DataFrame((X1-X2))" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
012
0-1.154342-0.1215160.128946
1-1.0619710.0551600.015198
2-0.9141530.177847-0.068929
3-1.0045410.1945130.050337
4-1.023858-0.003016-0.023659
\n", "
" ], "text/plain": [ " 0 1 2\n", "0 -1.154342 -0.121516 0.128946\n", "1 -1.061971 0.055160 0.015198\n", "2 -0.914153 0.177847 -0.068929\n", "3 -1.004541 0.194513 0.050337\n", "4 -1.023858 -0.003016 -0.023659" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_diff.head()" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Rejeitamos H0\n", "Valor da estatística 5666.660411153529\n", "valor p 1.1102230246251565e-16\n" ] } ], "source": [ "mu0=[0,0,0]\n", "n=len(df_diff)\n", "p=len(df_diff.columns)\n", "T2Hotelling(df_diff, mu0, n, p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }