6 Measuring agreement and Model evaluation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Measuring agreement and Model evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem:\n",
    "\n",
    "Consider the Grubbs model:\n",
    "\n",
    "(a) Show that the method of the moments provides the following estimators of its paramters: <br>\n",
    "\n",
    "   $$\n",
    "    \\left( \\begin{array}{cccc}\n",
    "    \\mu_b \\\\\n",
    "    \\beta_0+\\mu_b \\end{array} \\right) =\n",
    "    \\left( \\begin{array}{cccc}\n",
    "    \\bar{Y}_{.1} \\\\\n",
    "    \\bar{Y}_{.2} \\end{array} \\right)\n",
    "   $$<br>\n",
    "   $$\n",
    "    \\left( \\begin{array}{cccc}\n",
    "    \\sigma_b^2+\\sigma_{e1}^2 & \\sigma_b^2 \\\\\n",
    "    \\sigma_b^2 & \\sigma_b^2+\\sigma_{e2}^2 \\end{array} \\right) =\n",
    "    \\left( \\begin{array}{cccc}\n",
    "    {S_1}^2 & S_{12} \\\\\n",
    "    S_{12} & {S_2}^2 \\end{array} \\right)\n",
    "   $$ \n",
    "   \n",
    "   <br>\n",
    "   \n",
    "   Therefore, we will have:\n",
    "   $$ \\mu_b=\\bar{Y}_{.1} $$ <br>\n",
    "   $$ \\beta_0+\\mu_b=\\bar{Y}_{.2} $$ <br>\n",
    "   $$ \\sigma_b^2+\\sigma_{e1}^2={S_1}^2 $$ <br>\n",
    "   $$ \\sigma_b^2+\\sigma_{e2}^2={S_2}^2 $$ <br>\n",
    "   $$ \\sigma_b^2=S_{12} $$ <br>\n",
    "   \n",
    "   So the following estimators of its parameters: <br>\n",
    "   $ \\hat{\\mu}_b=\\bar{Y}_{.1}$, $\\hat{\\beta}_0=\\bar{Y}_{.2}-\\bar{Y}_{.1}$, ${\\hat{\\sigma}}_b^2=S_{12}$, ${\\hat{\\sigma}}_{e1}^2=S_1^2-S_{12}$, ${\\hat{\\sigma}}_{e2}^2=S_2^2-S_{12}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) When the error variance estimators are positive, we have:\n",
    "\n",
    "$${\\hat{\\sigma}}_{e1}^2=S_1^2-S_{12}>0$$\n",
    "\n",
    "$${\\hat{\\sigma}}_{e2}^2=S_2^2-S_{12}>0$$\n",
    "\n",
    "and apply to the definitions:\n",
    "\n",
    "$${\\hat{\\sigma}}_{e1}^2=S_1^2-S_{12}=\\frac{1}{n-1}\\sum_{i=1}^{n}(Y_{i1}-\\bar{Y}_{.1})^2-\\frac{1}{n-1}\\sum_{i=1}^{n}{(Y_{i1}-\\bar{Y}_{.1})(Y_{i2}-\\bar{Y}_{.2})}>0$$\n",
    "\n",
    "$${\\hat{\\sigma}}_{e2}^2=S_2^2-S_{12}=\\frac{1}{n-1}\\sum_{i=2}^{n}(Y_{i2}-\\bar{Y}_{.2})^2-\\frac{1}{n-1}\\sum_{i=1}^{n}{(Y_{i1}-\\bar{Y}_{.1})(Y_{i2}-\\bar{Y}_{.2})}>0$$\n",
    "\n",
    "which means the auto-covariance between the $S_1$ should be larger than the covariance between these two methods $S_{12}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IPI from 52 kidneys using CT (method 1) and urography (method 2).\n",
    "\n",
    "(a) Scatterplot and a Bland-Altman plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "A=np.array([97,77,74,59,79,85,78,78,68,96,74,64,76,60,78,71,67,103,95,78,70,80,78,102,102,77,45,60,50,94,91,66,63,65,58,75,105,65,80,90,58,75,83,78,85,65,90,76,100,65,40,53])\n",
    "B=np.array([100,58,95,55,79,95,60,88,68,94,60,64,88,57,66,67,76,95,85,105,80,85,82,102,100,75,40,70,63,103,95,80,72,68,48,70,90,60,80,96,54,80,88,70,90,79,100,85,108,53,58,49])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHbBJREFUeJzt3XuUHOV55/HvD3HxgO0IxMBBA0IiZgUYFgQTfNGaJeBjDCagyMbBm8TC4VibXW98SYIRu9mEZPGiLN7YSTbBkcFGto/FzVhg2AUTyV7HOOAdIQziooBBSBoJabgIY6zEAj/7R70DzdCj6enq7rrM73NOn+6qrq566p2eZ2qeeustRQRmZlZfexQdgJmZdZcTvZlZzTnRm5nVnBO9mVnNOdGbmdWcE72ZWc050VupSdog6d0dWlefpG9Jel7SDZ1YZzd0cp/zknSppK91aF0XSPp+J9Zlk+NEX0OS/o2kH6SE9qykuyT9Ss51vu6XVNI1ki7LF21nSDpV0uYJFvsAcDAwIyLO60FYldJiG1oF7Vl0ANZZkt4M3Ar8B+B6YG/gXcC/FBlXM5L2jIiXerjJw4F/amebBcRq1jE+oq+ffwUQESsi4uWI2BkR346I+0cXkPRRSQ9LekHSQ5JOTPOXSPpxw/xfT/OPBr4AvEPSTyXtkLQY+E3g02net9KyMyV9Q9KIpCckfbxhu5dKulHS1yT9BLigYd51abv3Sjq+2Y5J2kfS5yVtSY/Pp3n7Af8HmJli+amkmWM++6fAHwO/kd6/UNIekv5I0pOStkv6iqRfSsvPlhRpuY3A6nFiOlvSfalNfiDpXze817Q9J/o5JCdIuj/9V3adpDeMs/0L0n9sn0sxPC7pnWn+prRfi8a04WclbZS0TdIXUklrd224d2qbFyQ9KGmwYX1HS/pu2vaDks5peG+GpFsk/UTSD4FfbrYP1gMR4UeNHsCbgWeA5cCZwP5j3j8PGAZ+BRDwFuDwhvdmkh0A/AbwInBIeu8C4Ptj1nUNcFnD9B7AGrKEujdwBPA4cEZ6/1JgF7AgLdvXMO8DwF7AHwJPAHulz2wA3p1e/xlwN3AQ0A/8APhv6b1Tgc0TtM2lwNcapn8HeCzF+UbgJuCr6b3ZQABfAfYD+pqs70RgO/A2YBqwKMW7Twvtubufwwbgh+mzBwAPA787zj5dALwEfCTFcBmwEfgbYB/gPcALwBvT8p8HbknrfRPwLeDy8dowtdk/A2el9V8O3J3e2yu1339OP+/T0rbmpvevJfuvcj/g2LS/39/dz8iPLuWFogPwows/VDiaLAlvTkngFuDg9N4dwCdaXM99wLnp9QVjf0l5faJ/G7BxzDKXAF9Ory8Fvjfm/UtHE0ea3gPYCrwrTW/g1UT/Y+CshmXPADak169LUk3251Jem+hXAf+xYXou2R+dPXk10R+xm/VdSfpD0zBvPfBvW2jPcX8OaZ9/q2H6fwBfGGfZC4BHG6aPS3Ef3DDvGeAEsj8oLwK/3PDeO4AnxmvD1GZ/3zB9DLAzvX4X8BSwR8P7K9JnpqW2PKrhvf8+9jvkR28ertHXUEQ8TJYAkHQU8DWyI7kPAYeRJczXkfRh4PfJkhxkR7kHTmLTh5P967+jYd404B8apjc1+dwr8yLiF+mE4Mwmy80EnmyYfnKc5VrVbH17kp2wfV1sTRwOLJL0ew3z9h6NaYL2HPfnkDzV8Ppn7H4/tzW83gkQEWPnvZHsv6B9gTWSRt8T2c9od8bG8gZJe6aYNkXELxrefxIYSNvak9e2X2NbWw850ddcRDwi6Rrg36dZm2hSK5V0OPBF4HTgHyPiZUn3kSUCyI4SX7f6MdObyI4Oj9xdSE3mHdYQxx7AocCWJsttIUuuD6bpWQ3LtTMM6+j6Rs0i+w9oW4phovVuAj4TEZ8Z+0YL7dn059BlT5Ml/bdGxHCT9yfbhluAwyTt0ZDsZwH/BIyQteVhwCMN71kBfDK2ZiQdJekPJB2apg8jO5K/Oy1yFfCHkk5S5i0pKe1H9os+kj73EbK66qhtwKGS9h4z74iG6R8CP5F0cTrBN03SsZq4a+dJkhamo8RPkvUQurvJciuAP5LUL+lAsnMBo328twEzRk+mtmgF8ClJcyS9kay0cF203rvmi8DvSnpbasv9JL1P0puYuD3H+zl0TUrGXwQ+J+mgFNeApDPSIpNtw3vISkGflrSXpFOBXwOujYiXyc55XCppX0nHkJ3DsAI40dfPC2S18nskvUiWMNcBfwAQETcAnwG+npZdCRwQEQ8B/xP4R7Jf+OOAuxrWu5rsSPopSU+neVcDx6QeFyvTL/evkdWDnyA7grwKmChx3Ex2svI54LeBhRGxq8lylwFDwP3AA8C9aR4R8QhZ4n48xdNKSedLwFeB76V4/xn4vd1+okFEDAEfBf5Xiv0xUslsovYc7+fQ6rZzuDjFebeynk9/T3ZuYtJtGBE/B84hO+n/NPC3wIfTegD+E1nJ6Cmy8zlf7vjeWEuUTpKYFULSpcBbIuK3io7FrK58RG9mVnNO9GZmNefSjZlZzfmI3sys5krRj/7AAw+M2bNnFx2GmVmlrFmz5umI6J9ouQkTvaQvAWcD2yPi2DTvPLLLnI8GTk7dzEaXvwS4EHgZ+HhE3DHRNmbPns3Q0NBEi5mZWQNJLV1t3Erp5hrgvWPmrQMWkvU/btzoMcD5wFvTZ/5W0kSXV5uZWRdNmOgj4nvAs2PmPRwR65ssfi7ZVXH/EhFPkF2YcXJHIjUzs7Z0+mTsAK8dxGhzmvc6khZLGpI0NDIy0uEwzMxsVKcTvZrMa9p/MyKWRcRgRAz29094LsHMzNrU6US/mYaRCBl/FEIzM+uRTnevvAX4uqS/IBur+kiyEQ3NzCpl5dphrrhjPVt27GTm9D4uOmMuC+Y1rUSXXivdK1eQ3XnmwHRDiD8hOzn712Q3F7hN0n0RcUZEPCjpeuAhsrGoP5ZGNDQzq4yVa4e55KYH2LkrS1/DO3ZyyU0PAFQy2ZdiCITBwcFwP3ozK4v5S1czvGPn6+YPTO/jriWnFRBRc5LWRMTgRMt5CAQzszG2NEnyu5tfdk70ZmZjzJzeN6n5ZedEb2Y2xkVnzKVvr9de1N+31zQuOmNuQRHlU4pBzczMymT0hOuU6XVjZjYVLZg3UNnEPpZLN2ZmNecjerMC1emiHCsvJ3qzgtTtohwrL5duzApyxR3rX0nyo3buepkr7mg2ArhZ+5zozQpSt4tyrLyc6M0KUreLcqy8nOjNClK3i3KsvHwy1qwgdbsop9vcQ6l9TvRmBarTRTnd5B5K+bh0Y2al5x5K+fiI3sxKr8w9lKpQUvIRvZmVXll7KI2WlIZ37CR4taS0cu1woXGN5URvZqVX1h5KVSkpuXRjZqVX1h5KZS4pNXKiN7NKKGMPpZnT+5reW7boktJYE5ZuJH1J0nZJ6xrmHSDpTkmPpuf903xJ+itJj0m6X9KJ3QzezKxIZS0pjdVKjf4a4L1j5i0BVkXEkcCqNA1wJnBkeiwGruxMmGZWVSvXDjN/6WrmLLmN+UtX9/xEZTe3v2DeAO8/aYBpEgDTJN5/Uvn+85gw0UfE94Bnx8w+F1ieXi8HFjTM/0pk7gamSzqkU8GaWbUU3Sul29tfuXaYb6wZ5uUIAF6O4BtrhmvT6+bgiNgKkJ4PSvMHgE0Ny21O88xsCiq6V0q3t1/0/rWq090r1WReNF1QWixpSNLQyMhIh8MwszIouldKt7df9P61qt1Ev220JJOet6f5m4HDGpY7FNjSbAURsSwiBiNisL+/v80wzKzMir7QqdvbL3r/WtVuor8FWJReLwJubpj/4dT75u3A86MlHjObeoruldLt7Re9f62asB+9pBXAqcCBkjYDfwIsBa6XdCGwETgvLf6/gbOAx4CfAR/pQsxmVhFFX+jU7e0XvX+tUkTTEnpPDQ4OxtDQUNFhmJlViqQ1ETE40XIe68bMrOac6M3Mas6J3sys5jyomZl1VRVuzFF3TvRm1jW+12s5uHRjZl1TlSEC6s5H9GZtqFM5opv7Mt5QAMM7djJnyW2Vb7tmyvjdcKI3m6Q6lSO6vS/j3ZgDeM1okp3aXtHK+t1w6cZskupUjuj2vjQbImCsqrZdM2X9bviI3mySqjJiYSu6vS9jhwgY7zr8TrZdkaWTsn43fERvNklVGbGwFb3YlwXzBrhryWk8sfR9DHR5e0Xf6KSs3w0nerNJqsqIha3o9b50e3tFl07K+t1w6cZskqoyYmErer0v3d5e0aWTsn43PHqlmdXG/KWrm/byGZjex11LTisgou7y6JVmNuWUtXRSNJduzKw2ylo6KZoTvZnVyoJ5A1M+sY/l0o2ZWc050ZuZ1ZwTvZlZzeVK9JI+IWmdpAclfTLNO0DSnZIeTc/7dyZUMzNrR9snYyUdC3wUOBn4OXC7pNvSvFURsVTSEmAJcHEngjWzzivjsLpTQS/bPU+vm6OBuyPiZwCS/i/w68C5wKlpmeXAd3GiNyulsg6rW3e9bvc8pZt1wCmSZkjaFzgLOAw4OCK2AqTng/KHaWbdUPTYMFNVr9u97SP6iHhY0p8DdwI/BX4EvNTq5yUtBhYDzJo1q90wzCyHoseGmap63e65TsZGxNURcWJEnAI8CzwKbJN0CEB63j7OZ5dFxGBEDPb39+cJw8zaVNZhdeuu1+2et9fNQel5FrAQWAHcAixKiywCbs6zDbOqWrl2mPlLVzNnyW3MX7q66ZjorSzTzRg8Nkwxet3ueYdA+IakGcAu4GMR8ZykpcD1ki4ENgLn5Q3SrGpaOdnW7RNyrazfY8MUo9ft7mGKzbqgleFyuz2k7lQbsncq8jDFZgVq5WRbt0/I+USrjXKiN+uCVk62dfuEnE+02ignerMuaOVkW7dPyPlEq43yePRmXdDKybZun5DziVYb5ZOxZmYV5ZOxZmYGuHRj1jWdHJ3QI0w253ZpjRO9WRd08mIojzDZnNuldS7dmHVBJ0cn9AiTzbldWudEb9YFnbxYyRc+Ned2aZ0TvVkXdPJiJV/41JzbpXVO9GZd0MmLlXzhU3Nul9b5ZKxZF3TyYiVf+NSc26V1vmDKzKyifMGUmZkBLt2YFcoX/FgvONGbFcQX/FivuHRjVhBf8GO94kRvVhBf8GO94kRvVhBf8GO9kivRS/qUpAclrZO0QtIbJM2RdI+kRyVdJ2nvTgVrVie+4Md6pe2TsZIGgI8Dx0TETknXA+cDZwGfi4hrJX0BuBC4siPR2pRTlaF+21n3ZC74ce8cyyNvr5s9gT5Ju4B9ga3AacC/S+8vBy7Fid7aUJWhfvOse8G8gQmXce8cy6vt0k1EDAOfBTaSJfjngTXAjoh4KS22GfA30dpSlaF+u917xr1zLK88pZv9gXOBOcAO4AbgzCaLNh1jQdJiYDHArFmz2g3DamRseWK4IkP9dnrd3WwHm5rynIx9N/BERIxExC7gJuCdwHRJo39ADgW2NPtwRCyLiMGIGOzv788RhtXBaHlieMdOgqw8oXGWLdtQv51cd7fbwaamPIl+I/B2SftKEnA68BDwHeADaZlFwM35QrSpoFl5IuB1Sa6MQ/12ct3dbgebmvLU6O8BbgTuBR5I61oGXAz8vqTHgBnA1R2I02puvDJEAAPT+1B6vnzhcW0P9Xv5wuM6sq5urrvb7WBTk4cptlKYv3R101r0wPQ+7lpyWgERFcPtYJPhYYqtUnzxUMbtYN3g0SutFHy3oIzbwbrBpRszs4pqtXTjI3qrHA8HYDY5TvRWKR4OwGzyfDLWKsXDAZhNno/orVKqfrMOl52sCD6it0qp8s06mg1vcMlND7By7XDRoVnNOdFbpVS5n7nLTlYUl26sUqrcz7zqZSerLid6q5xWbtZRRuMNOVyFspNVm0s3Zj1S5bKTVZuP6M16pMplJ6s2J3qzHqpq2cmqzaUbM7Oa8xG99YQvFDIrjhO9dZ3HpzErlks31nW+UMisWE701nW+UMisWE701nVVHp/GrA6c6K3rfKGQWbHaPhkraS5wXcOsI4A/Br6S5s8GNgAfjIjn2g/RyqSd3jNT4UKhdnsVuTeS9UJH7hkraRowDLwN+BjwbEQslbQE2D8iLt7d533P2GoY23sGsiPzyxceN6WTU7vt4va0vFq9Z2ynSjenAz+OiCeBc4Hlaf5yYEGHtmEFc++Z5tptF7en9UqnEv35wIr0+uCI2AqQng9q9gFJiyUNSRoaGRnpUBjWTe4901y77eL2tF7Jnegl7Q2cA9wwmc9FxLKIGIyIwf7+/rxhWA+490xz7baL29N6pRNH9GcC90bEtjS9TdIhAOl5ewe2YSXg3jPNtdsubk/rlU4MgfAhXi3bANwCLAKWpuebO7ANK4Fu957J0wOlyN4r7bbLVOiNZOWQq9eNpH2BTcAREfF8mjcDuB6YBWwEzouIZ3e3Hve6sTw9UNx7xaaqVnvd5Dqij4ifATPGzHuGrBeOWct21wNlomSd57Pux25TgUevtFLI0wOl3c96VE2bKjwEgpVCnh4o7X7W/dhtqnCit1LI0wOl3c+WoR/7yrXDzF+6mjlLbmP+0tWsXDvcs23b1OHSjZVCnh4o7X525vQ+hpsk9V71Y3fpyHqlI2Pd5OVeN1aEonvrzF+6uukfmoHpfdy15LSub9+qrye9bsyqrOh+7GUoHdnU4ERvU9qCeQOFlUmKLh3Z1OGTsWYF8RAI1is+orfCTPWLlYouHdnU4URvhXCPk0yRpSObOly6sUL4YiWz3vERfQGmeskC3OPErJd8RN9joyWL4R07CV4tWUy1KyJ90w2z3nGi7zGXLDLucWLWOy7d9JhLFhn3ODHrHSf6HvNFMq9yjxOz3nDppsdcsjCzXvMRfY+5ZGFmveZEXwCXLMysl1y6MTOruVxH9JKmA1cBxwIB/A6wHrgOmA1sAD4YEc/litKmhCpdSNapWKu0z1ZdeY/o/xK4PSKOAo4HHgaWAKsi4khgVZo2260qXUjWqVirtM9WbW0neklvBk4BrgaIiJ9HxA7gXGB5Wmw5sCBvkFZ/VbqQrFOxVmmfrdryHNEfAYwAX5a0VtJVkvYDDo6IrQDp+aBmH5a0WNKQpKGRkZEcYVgdVOlCsk7FWqV9tmrLk+j3BE4EroyIecCLTKJMExHLImIwIgb7+/tzhGF1UKWxbzoVa5X22aotT6LfDGyOiHvS9I1kiX+bpEMA0vP2fCHaVFClC8k6FWuV9tmqre1eNxHxlKRNkuZGxHrgdOCh9FgELE3PN3ck0jHcW6FeqnQhWZ5Yx35v33/SAN95ZKT0+2zVpoho/8PSCWTdK/cGHgc+QvZfwvXALGAjcF5EPLu79QwODsbQ0FDL2x17dyLIjoQuX3icf0mstPy9tU6TtCYiBidaLlc/+oi4D2i2kdPzrHciu+ut4F+Y7vJ/Uu3z99aKUskhENxboRi+z2s+/t5aUSo5BIJ7KxTD/b7z8ffWilLJRO/eCsXwEWk+/t5aUSpZuqlSD4068U1T8vH31oqSq9dNp0y2140Vw71GzMql1V43lSzdWG+sXDvM/KWrmbPkNuYvXQ3A5QuPY2B6HwIGpvc5yZtVQCVLN9Z94/WwuXzhcdy15LSCozOzyfARvTXlHjZm9eEjentF48VQ4525cQ8bs+pxojeg+YnWZtzDxqx6XLoxoHmpZiz3+TarJh/RG7D7koygJ32+PY6OWXc40Rsw/sVQA9P7etLLxuPomHWPSzcGFH95vnv5mHWPj+gNKP7yfI+jY9Y9TvT2igXzBgork3gcHbPucenGSqHo0pFZnfmI3kqh6NKRWZ050VtpFFk6Mqszl27MzGou1xG9pA3AC8DLwEsRMSjpAOA6YDawAfhgRDyXL0wzM2tXJ47ofzUiTmgY/H4JsCoijgRWpWkzMytIN0o35wLL0+vlwIIubMPMzFqU92RsAN+WFMDfRcQy4OCI2AoQEVslHdTsg5IWA4sBZs2alTOMqcFjwZhZO/Im+vkRsSUl8zslPdLqB9MfhWWQ3TM2Zxy157FgzKxduUo3EbElPW8HvgmcDGyTdAhAet6eN0jzWDBm1r62E72k/SS9afQ18B5gHXALsCgttgi4OW+Q1puxYMbeDHzl2uGOrdvMipOndHMw8E1Jo+v5ekTcLun/AddLuhDYCJyXP0zr9lgwLg2Z1VfbR/QR8XhEHJ8eb42Iz6T5z0TE6RFxZHp+tnPhTl3dHgvGpSGz+vIQCBXR7bFgPEywWX050VdIN8eC8TDBZvXlsW4M8DDBZnXmI3oDPEywWZ050dsrPEywWT25dGNmVnNO9GZmNedEb2ZWc070ZmY150RvZlZzTvRmZjXnRG9mVnNO9GZmNedEb2ZWc070ZmY150RvZlZzTvRmZjXnRG9mVnNO9GZmNedEb2ZWc7nHo5c0DRgChiPibElzgGuBA4B7gd+OiJ/n3U6vrVw77JtwmFktdOKI/hPAww3Tfw58LiKOBJ4DLuzANnpq5dphLrnpAYZ37CSA4R07ueSmB1i5drjo0MzMJi1Xopd0KPA+4Ko0LeA04Ma0yHJgQZ5tFOGKO9azc9fLr5m3c9fLXHHH+oIiMjNrX94j+s8DnwZ+kaZnADsi4qU0vRloWu+QtFjSkKShkZGRnGF01pYdOyc138yszNpO9JLOBrZHxJrG2U0WjWafj4hlETEYEYP9/f3thtEVM6f3TWq+mVmZ5Tminw+cI2kD2cnX08iO8KdLGj3JeyiwJVeEBbjojLn07TXtNfP69prGRWfMLSgiM7P2tZ3oI+KSiDg0ImYD5wOrI+I3ge8AH0iLLQJuzh1ljy2YN8DlC49jYHofAgam93H5wuPc68bMKil398omLgaulXQZsBa4ugvb6LoF8wac2M2sFjqS6CPiu8B30+vHgZM7sV4zM8vPV8aamdWcE72ZWc050ZuZ1ZwTvZlZzSmi6fVMvQ1CGgGebPPjBwJPdzCcXnLsxXDsxahq7GWO+/CImPCK01Ik+jwkDUXEYNFxtMOxF8OxF6OqsVc17kYu3ZiZ1ZwTvZlZzdUh0S8rOoAcHHsxHHsxqhp7VeN+ReVr9GZmtnt1OKI3M7PdcKI3M6u5SiZ6SdMkrZV0a5qeI+keSY9Kuk7S3kXH2IykDZIekHSfpKE07wBJd6bY75S0f9FxNiNpuqQbJT0i6WFJ7yh77JLmprYeffxE0ifLHvcoSZ+S9KCkdZJWSHpDhb7rn0hxPyjpk2leKdtd0pckbZe0rmFe01iV+StJj0m6X9KJxUXeukomeqp9Q/JfjYgTGvrlLgFWpdhXpeky+kvg9og4CjierP1LHXtErE9tfQJwEvAz4JuUPG4ASQPAx4HBiDgWmEZ234fSf9clHQt8lGwU2+OBsyUdSXnb/RrgvWPmjRfrmcCR6bEYuLJHMeYTEZV6kN21ahXZHa1uJbt94dPAnun9dwB3FB3nOLFvAA4cM289cEh6fQiwvug4m8T9ZuAJ0sn7KsXeEOt7gLuqEjfZvZY3AQeQDSd+K3BGFb7rwHnAVQ3T/5Xs3tKlbXdgNrCuYbpprMDfAR9qtlyZH1U8om/7huQlEMC3Ja2RtDjNOzgitgKk54MKi258RwAjwJdTyewqSftRjdhHnQ+sSK9LH3dEDAOfBTYCW4HngTVU47u+DjhF0gxJ+wJnAYdRgXZvMF6so3+AR5X1Z/AalUr0eW9IXgLzI+JEsn//PibplKIDatGewInAlRExD3iR8vzbPaFUxz4HuKHoWFqVasLnAnOAmcB+ZN+bsUr3XY+Ih8lKTHcCtwM/Al7a7Yeqo0r55hWVSvRU/IbkEbElPW8nqxWfDGyTdAhAet5eXITj2gxsjoh70vSNZIm/CrFDliDvjYhtaboKcb8beCIiRiJiF3AT8E6q812/OiJOjIhTgGeBR6lGu48aL9bNZP+djCrtz6BRpRJ9VPiG5JL2k/Sm0ddkNeN1wC1kMUNJY4+Ip4BNkuamWacDD1GB2JMP8WrZBqoR90bg7ZL2lSRebfPSf9cBJB2UnmcBC8navwrtPmq8WG8BPpx637wdeH60xFNqRZ8kyHHy5FTg1vT6COCHwGNk/57vU3R8TeI9guxf2B8BDwL/Jc2fQXZy+dH0fEDRsY4T/wnAEHA/sBLYvwqxA/sCzwC/1DCv9HGnOP8UeITsgOCrwD5V+K6n2P+B7A/Tj4DTy9zuZH+EtgK7yI7YLxwvVrLSzd8APwYeIOsVVfg+TPTwEAhmZjVXqdKNmZlNnhO9mVnNOdGbmdWcE72ZWc050ZuZ1ZwTvZlZzTnRm5nV3P8HHp4CdyEEeJAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c1d90a470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X98VXed5/HXpzS1SVtIS1BKoAJuYS2goZtWAQccYSU6YIHHiK08tIzudH2su6O7bmfKtKt1p92y2x236uy49qEOjkZr7AAWqqkFFEZBJUhqwRas0AqBtBAFrGE0pZ/945yUm+QmueH+OL/ez8fjPnLv+d5zzvd7z7mfnPv5nvM95u6IiEj6XRR1BUREpDIU8EVEMkIBX0QkIxTwRUQyQgFfRCQjFPBFRDJCAT+BzGydmd1TpmU/a2aLopo/KpWst5ndY2YnzayzEuuLQin3UTO728y+WoplZZ0CfgyFweesmb1oZr8xs0fNbFLU9eplZpeFdfv2MO+L1RfVzNzMfhfWvcPMPmVmo0a4jLea2dEi6jAJ+BhwnbuPL/Xyo2Bmq83sB1HXQ4angB9fS939cuBq4HngsxHXJ9efAr8H3m5mV0ddmRF6Y/i5LgTeC/x5hdf/WqDL3V+o8HpFFPDjzt3/BXgYuC5fuZldaWabzexE+Gtgs5lNzCn/vpn9jZn90Mx+a2bfNbO6nPL3mdlzZtZlZncWWK1bgf8H/AxYNUi9moC/Bt4THlE/kVOfe8xsZzh9k5mNNbNmMztjZrvNbHLOcj5tZkfCsj1m9kc5ZXebWYuZ/WPYtv1m1lhIA9z9aeCfgZl56v4qM3vAzI6FjwfCaZcB3wEmhHV/0cwm5Jl/TFinE+Fne5eZXRSmjB7PmX9dv/nyLj/8tVcXvucuM3vJzEaHr+8xsweGWu8g2+duM/ummX01/OyeNLNpZrbGzF4IP/O392vTF83sePjr6B4zG2VmryfYF+aE9T2Vs5orw1+nvzWzH5vZ63KWNzfc1qfDv3NzyqaY2fZwvseB3P310rDOXWZ2Kpz3NYNsZulHAT/mzKwGeA/wo0HechHwDwRHjtcAZ4G/6/ee9wJ/BrwauAT4r+GyrwM+B7wPmACMBSYyBDO7Bngr0Bw+3p/vfe7eCvwP4Bvufrm7vzGn+OZwnfXA64BdYRuuAp4CPpHz3t1AQ1j2NeCbZnZpTvm7gIeAWuCRPG0frB3XAX8E7M1TfCfw5nC9bwRuBO5y998B7wCOhW263N2P5Zn/s8AYYCqwgOAz+jN339Jv/tW5Mw2x/N3hcgDmA88B83Jebx9qvUN8DEuBrwBXhp/DYwT7Uz3w34HP57z3y8BLwL8CZgNvB/6duz8FfAjYFda3NmeeW4BPhst/BrgXwMyuAh4FPkOwz30KeNTMxobzfQ3YQxDo/4bgAKPXrWEbJ4Xzfohgn5dCuLseMXsAzwIvAqcIvmTHgFk55euAewaZtwH4Tc7r7xMEq97X/wFoDZ9/HHgop+wy4A/AoiHqdhfQHj6fAJwDZver+6Lw+d3AV/vN/33gzpzXfwt8J+f10t7lD7L+3xCkZXqXvyWn7Drg7BDzOnAmXMYvgXuAi/LU+5fAO3PmWww8Gz5/K3B0iHWMIkh3XZcz7d8D3y9w/gHlBEHvM8DFQCfwEWAtcClBsKsbbr151nM38Hi/z/1FYFT4+orw86oFXhMuuzrn/bcA3wufrwZ+0G/564Av5Lx+J/B0+Px9wE/6vX9XuJxrCPb5y3LKvta7HwEfAHYCb4j6e5rEh47w42uZB0dLrwL+I7DdzPJ18tWY2efDn/BngB1ArfXtjMw9G6QbuDx8PgE40lvgwRFmV86yX8x5XBNOfj/BkT0eHH1up+8RWCGez3l+Ns/r3vphZh8zs6fCn/6nCI7u6nLe379tl5rZxUOs+3p3v9LdX+fud7n7y3neM4HgKLrXc+G0QtQR/IrqP399gfPns53gH8H1wJMEaaEFBL9CnnH3kxe43v6f+0l3P5fzGoJt8VqgCjgeplFOERz9v3qYeg+13z3X7729dZ1AcMDyu35lvb5C8EvkoTDd9r/MrGqYekhIAT/m3P2cu68nOJJ+S563fAyYDrzJ3UcT/MQHsAIWf5zgp3EwQ5A+6v1ZjZ9PK1zu7r8K86zXAmvMrNOC0wrfBNwySJAtaijWMF//V8BK4MrwH+BpCmtbMY4RBLle14TTYPg2nQR68szfUeC68y1/J8E2Xg5sd/efh8v8E86nc4pd71COEBzh17l7bfgY7e4zhqjzUPp/vnC+rscJcv+X9SsLVuTe4+6fdPfrgLnAEgZJK8pACvgxZ4GbCPKgT+V5yxUER2OnwtzoJ/K8ZzAPA0vM7C1mdglB3naofeJWgqPL6whSRw0EnZ41BLnn/p4HJg/WcViAKwh+3p8ALjazjwOjL3BZI/F14C4zGxd2ln4c6D299HlgrJmNyTdjeITcAtxrZleY2WuB/5Iz/3AGLN/duwly2h/mfIDfSZCy2V6i9Q7K3Y8D3wX+1sxGhx3QrzOz3n6F54GJ4T5UiG8D08zsvWZ2sZm9h2Cf2uzuzwFtwCfN7BIzewtBugkAM/tjM5sV/oI9Q/BP7lyedUgeCvjxtcnMXiTYqe8FbnX3/Xne9wBQTXCE9yOgtdAVhMv7MEGO9DhBbjvvOeBhR+lK4LPu3pnzOEzwMztfWueb4d8uM/tpofXK8RjBWSsHCX7W/ws5Kagyuocg6PyMIIXy03AaHpzd83XgUJjeyJfq+U/A74BDwA8IPt8vFbLiIZa/nSCt8pOc11cQpPCKXm8B3k+QMvo5wX7yMMEpwwDbgP1Ap5mdHG5B7t5FcGT+MYIU4l8CS8LUFAQnGbwJ+DXBAcw/5sw+Plz3GYIDoO2U4J9aVljYESIiIimnI3wRkYxQwBcRyQgFfBGRjFDAFxHJiKEuUKm4uro6nzx5ctTVEBFJlD179px093HDvS9WAX/y5Mm0tbVFXQ0RkUQxs/5XLuellI6ISEYUHfDNbJKZfS8c72S/mX0knH6VmT1uZr8I/15ZfHVFRORCleII/yXgY+7+eoLBnD4cDj17B7DV3a8FtoavRUQkIkUHfHc/7u4/DZ//luBy53rgJoIxtAn/Lit2XSIicuFKmsO34E5Fs4EfA68JB13qHXwp71CqZnabmbWZWduJEydKWR0REclRsrN0zOxy4J+Aj7r7GbPCRrB19weBBwEaGxs1sI8MaePeDu5/7ADHTp1lQm01ty+ezrLZxQw1L5IdJQn44Q0I/gloDsduB3jezK529+MW3OhaN22Womzc28Ga9U9yticYDbfj1FnWrH8SQEFfpAClOEvHgC8CT7n7p3KKHuH8kLm3At8qdl2Sbfc/duCVYN/rbM857n/sQEQ1EkmWUhzhzyO4R+WTZtYeTvtrgntutpjZB4FfAe8uwbokw46dyn+v6sGmi0hfRQd8d/8Bg99ybmGxyxfpNaG2mo48wX1CbXUEtRFJHl1pK4lx++LpVFeN6jOtumoUty+e/srrjXs7mLd2G1PueJR5a7excW8pbukqkg6xGktHZCi9HbODnaWjTl2RoSngS6Ism10/aPAeqlNXAV8kZgG/q6uLdevW9Zk2Y8YMbrjhBnp6emhubh4wT0NDAw0NDXR3d9PS0jKgvLGxkZkzZ3L69Gk2bNgwoHzOnDlMnz6dkydPsnnz5gHl8+fPZ+rUqXR2dtLaOvD+4AsXLmTSpEkcOXKErVu3Dihvampi/PjxHDp0iB07dgwoX7JkCXV1dRw4cIBdu3YNKF++fDljxoxh3759eUcSXblyJTU1NbS3t9Pe3j6gfNWqVVRVVbF792727x94D/TVq1cDsHPnTg4ePNinrKqqilWrVgGwfft2Dh8+3Ke8pqaGlStXArBlyxaOHu17//PRo0ezYsUKAFpbW+ns7OxTPnbsWJYuXQrApk2b6Orq6lM+fvx4mpqaAFi/fj1nzpzpUz5x4kQWLVoEQEtLC2/oPsKsS86XH395NE+8NIFjp87S3NxMT09Pn/mnTZvG3LlzAQbsd6B9T/te4fted3d3n/IpU6awYMECgEj2vcEohy+pccnFo/JOV6euSMDc43Nxa2Njo2s8fLlQ/XP4EHTq3rdillI6kmpmtsfdG4d7X6xSOiLFGK5TVyTrFPAlVYbq1BXJOuXwRUQyQgFfRCQjFPBFRDJCAV9EJCPUaSsimZalm+oo4ItIZmVt/CWldEQks7J2Ux0FfBHJrKzdVEcBX0Qya7BxltI6/pICviSCbmwi5VDITXXSRJ22EntZ61iTysna+EsK+BJ7urGJlFOWxl9SSkdiL2sdayLlooAvsZe1jjWRclHAl9jLWseaSLkohy+xl7WONZFyUcCXRMhSx5pIuSjgC5CtAaTyyXr7JRsU8CXz57lnvf2SHeq0lcwNINVf1tsv2aGAL5k/zz3r7ZfsUMCXzJ/nnvX2S3Yo4Evmz3PPevslO9RpK5k/zz3r7ZfsMHePug6vaGxs9La2tqirISKSKGa2x90bh3ufUjoiIhlRkpSOmX0JWAK84O4zw2lXAd8AJgPPAivd/TelWF+a6QIgSSvt29Er1RH+OqCp37Q7gK3ufi2wNXwtQ+i9AKjj1Fmc8xcA6e5OknTat+OhJAHf3XcAv+43+Sbgy+HzLwPLSrGuNNMFQJJW2rfjoZw5/Ne4+3GA8O+r873JzG4zszYzaztx4kQZqxN/ugBI0kr7djxE3mnr7g+6e6O7N44bNy7q6kRKFwBJWmnfjodyBvznzexqgPDvC2VcVyroAiBJK+3b8VDOgP8IcGv4/FbgW2VcVyosm13PfStmUV9bjQH1tdXct2KWzmSQxNO+HQ8lufDKzL4OvBWoA54HPgFsBFqAa4BfAe929/4du33owisRkZEr9MKrkpyH7+63DFK0sBTLFxGR4mksHZER0MVD2ZK27a2AL1Ig3RkrW9K4vSM/LVMkKXTxULakcXsr4IsUSBcPZUsat7cCvkiBdPFQtqRxeyvgixRIFw9lSxq3tzptRQqkO2NlSxq3t+54JSKScLrjlYiI9KGALyKSEQr4IiIZoYAvIpIRiT9LJ21jXSRZWrZFWtoh0l+iA34ax7pIqrRsi7S0QySfRKd00jjWRVKlZVukpR0i+SQ64KdxrIukSsu2SEs7RPJJdMBP41gXSZWWbZGWdojkk+gc/u2Lp/fJt0K8x7ood2dgvuVDZS4Nj8u2KPYzjks7RMoh0QE/SWNdlLszMN/yb3/4CXDoednLss5ccdgWpfiM49AOkXLRWDoVMm/tNjry5IHra6v54R1vK9vy8ynVOuOm3J+xSFxpLJ2YKXdn4EiWk9YOSHW4igwt0SmdJJlQW5336LNUnYGDLX+w9/ZK00VG5f6MRZJOR/gVUu6bKeRbftUoo+oiG3SdvTnvjlNncc7nvDfu7ShJnSotjTesECklHeFXSLk7Awdb/lDrHOoioyQe5avDVWRo6rTNsCl3PEq+rW/A4bV/UunqiMgFKrTTVkf4Gaacd/mkqW8kSaL63JOyvZXDzzDlvMsjbX0jSRHV556k7a2An2HLZtdz34pZ1NdWYwTnq9+3YlYsj0ySRAOwRSOqzz1J21spnYxbNrteAb7EdD1ANKL63JO0vXWEL1JiGoAtGlF97kna3gr4IiWmvpFoRPW5J2l7K6UjUmK6HiAaUX3uSdreOg9fRCThNHiaiIj0UfaUjpk1AZ8GRgFfcPe15V6nxFOpL05JysUuIxXljWykdOK4f5Y1pWNmo4CDwL8FjgK7gVvc/ef53q+UTnr1vzkJBB1bF3ref6mXFxf52lV1kYFBz7nz39U0tDXNKr1/xiWlcyPwjLsfcvc/AA8BN5V5nRJDpb44JUkXu4xEvnb1vOx9gj2ko61pFtf9s9wBvx44kvP6aDjtFWZ2m5m1mVnbiRMnylwdiUqpL05J0sUuI6Eb2aRDXPfPcgd8yzOtz6GKuz/o7o3u3jhu3LgyV0eiUuqLU5J0sctIjKT+SW9rmsV1/yx3wD8KTMp5PRE4VuZ1SgyV+uKUJF3sMhJ5b2RzkVE1avAb2Uj8xHX/LPdZOruBa81sCtAB3Ay8t8zrlBgq9cUpSbrYZSQu5EY2Ej9x3T/LfuGVmb0TeIDgtMwvufu9g71XZ+mIiIxcbG6A4u7fBr5d7vWIiMjQYjWWTldXF+vWreszbcaMGdxwww309PTQ3Nw8YJ6GhgYaGhro7u6mpaVlQHljYyMzZ87k9OnTbNiwYUD5nDlzmD59OidPnmTz5s0DyufPn8/UqVPp7OyktbV1QPnChQuZNGkSR44cYevWrQPKm5qaGD9+PIcOHWLHjh0DypcsWUJdXR0HDhxg165dA8qXL1/OmDFj2LdvH/l+/axcuZKamhra29tpb28fUL5q1SqqqqrYvXs3+/fvH1C+evVqAHbu3MnBgwf7lFVVVbFq1SoAtm/fzuHDh/uU19TUsHLlSgC2bNnC0aNH+5SPHj2aFStWANDa2kpnZ2ef8rFjx7J06VIANm3aRFdXV5/y8ePH09TUBMD69es5c+ZMn/KJEyeyaNEiAFpaWuju7u5TPmXKFBYsWABAc3MzPT09fcqnTZvG3LlzAQbsd6B9T/tecve9wWhoBRGRjNDgaSIiCReXK21FRCQmYpXDF8kVx8GnJBvSuu8p4Ess9R98quPUWdasfxIgFV88ia8073tK6UgsxXXwKUm/NO97CvgSS3EdfErSL837nlI6EksTaqvpyPMFi3rwqaikNadcqEq2P837no7wJZbiOvhUFHpzyh2nzuKczylv3NsRddUqotLtT/O+p4AvsbRsdj33rZhFfW01BtTXVmf2Dk9pzikXotLtT/O+p5SOxNay2fWp+JIVK8055UJE0f607nsK+CIRGElO+kJyysXkvAudt1J59TTn1CtNKR2RChtpTnqkOeVict6FzlvJvHqac+qVpoAvUmEjzUmPNKdcTM670HkrmVdPc0690pTSEamwC8lJjySnXEzOu9B5K51XT2tOvdJ0hC9SYeW+wXUxyy903rjepFuGpoAvUmHlzkkXs/xC5610Xn3j3g7mrd3GlDseZd7abZm5BqHUlNIRqbBy3+C6mOUXOm8lb9Kd5sHMKk03QBGRWJu3dlve0zLra6v54R1vi6BG8aMboIhIKmT9wrNSUsAXkVhTB3HpKOCLSKzpwqvSUaetiMRaJTuI004BX0RiTxdelYYCvohIhCp5cxcFfBGRiFT6GgN12oqIRKTSN3dRwBcRiUilrzFQwBcRiUilrzFQDl8SpdQdXJXsMCuFpNU3qSr1Od++eHqfHD6U9xoDBXxJjFJ3cCVtUK6k1TepKvk5V/oaAw2eJolR6kG0kjYoV9Lqm1RJ/Jw1eJqkTqk7uJI2KFfS6ptUaf6ci0rpmNm7gbuB1wM3untbTtka4IPAOeAv3P2xYtYl6TSSXOmE2uq8R14X2sFVquVVKt97ofWNS96/EvUoxTpKvZ/FSbFH+PuAFcCO3Ilmdh1wMzADaAL+3sxGDZxdsqw3V9px6izO+VzpYHczKvUgWqVY3kjbUIwLqW8l6zeUStSjVOtI82BtRQV8d3/K3fNdIXAT8JC7/97dDwPPADcWsy5Jn5FedLJsdj33rZhFfW01RpBTvW/FrAs+SizF8ip54cyF1LfSF/ZEWY9SraPU+1mclOssnXrgRzmvj4bTBjCz24DbAK655poyVUfi6EJypaUeRKvY5VU63zvS+sYlH12JepRyHWkdrG3YI3wz22Jm+/I8bhpqtjzT8p4O5O4PunujuzeOGzeu0HpLCqThxhZxb0Nc6leJesSlrXE2bMB390XuPjPP41tDzHYUmJTzeiJwrNjKSrqkIVca9zbEpX6VqEdc2hpn5UrpPAJ8zcw+BUwArgV+UqZ1SUKl4cYWcW9DXOpXiXrEpa1xVtSFV2a2HPgsMA44BbS7++Kw7E7gA8BLwEfd/TvDLU8XXomIjFyhF14VdYTv7huADYOU3QvcW8zyRWRk4nLOvcSTxtIRSQmNtSPD0dAKIikRl3PuJb4U8EVSIi7n3Et8KeCLpITOQ5fhKOCn1Ma9Hcxbu40pdzzKvLXbKj52ilSezkOX4ajTNoXUeZdNOg9dhqOAn0JDdd7py59uaR0DRkpDKZ0UUuediOSjI/wUSvMNHMol3wVLoPSIpIsCfgrdvnh6nxw+qPNuKPn6PG5/+Alw6HnZX5mmfhBJOqV0UijNN3Aoh3x9Hj3n/JVg30sXMUnS6Qg/pdR5V7iR9G2oH0SSTEf4knkj6dtQP4gkmQK+ZF6+C5aqRhlVF/W9cZv6QSTplNKRzBvsgqV805QmkyQr6gYopaYboIiIjFyhN0BRSkdEJCMU8EVEMkIBX0QkIxTwRUQyQgFfRCQjFPBFRDJCAV9EJCMU8EVEMkIBX0QkIzS0Qgnku3mGLsEXkbhRwC+SbhguIkmhlE6RhrphuIhInCjgF0k3DBeRpFDAL9JgN8TQjTJEJG4U8IuU7+YZulGGiMSROm2LNNjNM9RhKyJxo4BfArphuIgkgVI6IiIZUVTAN7P7zexpM/uZmW0ws9qcsjVm9oyZHTCzxcVXVUREilHsEf7jwEx3fwNwEFgDYGbXATcDM4Am4O/NbNSgSxERkbIrKuC7+3fd/aXw5Y+AieHzm4CH3P337n4YeAa4sZh1iYhIcUqZw/8A8J3weT1wJKfsaDhtADO7zczazKztxIkTJayOiIjkGvYsHTPbAozPU3Snu38rfM+dwEtAc+9sed7v+Zbv7g8CDwI0NjbmfY+IlIcG/suWYQO+uy8aqtzMbgWWAAvdvTdgHwUm5bxtInDsQispIqWngf+yp6jz8M2sCfgrYIG7d+cUPQJ8zcw+BUwArgV+Mtzyurq6WLduXZ9pM2bM4IYbbqCnp4fm5uYB8zQ0NNDQ0EB3dzctLS0DyhsbG5k5cyanT59mw4YNA8rnzJnD9OnTOXnyJJs3bx5QPn/+fKZOnUpnZyetra0DyhcuXMikSZM4cuQIW7duHVDe1NTE+PHjOXToEDt27BhQvmTJEurq6jhw4AC7du0aUL58+XLGjBnDvn37aGtrG1C+cuVKampqaG9vp729fUD5qlWrqKqqYvfu3ezfv39A+erVqwHYuXMnBw8e7FNWVVXFqlWrANi+fTuHDx/uU15TU8PKlSsB2LJlC0ePHu1TPnr0aFasWAFAa2srnZ2dfcrHjh3L0qVLAdi0aRNdXV19ysePH09TUxMA69ev58yZM33KJ06cyKJFwfFIS0sL3d3dfcqnTJnCggULAGhubqanp6dP+bRp05g7dy7AgP0O0r/vNR+r42yPMemiU8y4+Py2+d6mpzn1RK32vQTve4Mp9sKrvwNeBTxuZgA/cvcPuft+M2sBfk6Q6vmwu58bYjkiUmEv/Pb3wKUDpv/hJX1V08rOZ2Gi19jY6PmOJESk9Oat3UZHnlFd62ur+eEdb4ugRnKhzGyPuzcO9z5daSuSURr4L3s0lo5IRmngv+xRwBfJMA38ly1K6YiIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEYo4IuIZIQCvohIRijgi4hkhAK+iEhGKOCLiGSEAr6ISEboFociEgsb93bo/rplpoAvIpHbuLeDNeuf5GzPOQA6Tp1lzfonART0S0gpHRGJ3P2PHXgl2Pc623OO+x87EFGN0kkBX0Qid+zU2RFNlwujgC8ikZtQWz2i6XJhFPBFJHK3L55OddWoPtOqq0Zx++LpEdUondRpKyKR6+2Y1Vk65aWALyKxsGx2vQJ8mSmlIyKSEQr4IiIZoYAvIpIRCvgiIhmhgC8ikhHm7lHX4RVmdgJ4LnxZB5yMsDqlkoZ2pKENkI52qA3xEad2vNbdxw33plgF/Fxm1ubujVHXo1hpaEca2gDpaIfaEB9JbIdSOiIiGaGALyKSEXEO+A9GXYESSUM70tAGSEc71Ib4SFw7YpvDFxGR0orzEb6IiJSQAr6ISEbEJuCb2Sgz22tmm8PXU8zsx2b2CzP7hpldEnUdh2Nmz5rZk2bWbmZt4bSrzOzxsB2Pm9mVUddzKGZWa2YPm9nTZvaUmc1JYBumh9ug93HGzD6awHb8ZzPbb2b7zOzrZnZpQr8XHwnbsN/MPhpOi/22MLMvmdkLZrYvZ1reelvgM2b2jJn9zMyuj67mg4tNwAc+AjyV8/p/Av/H3a8FfgN8MJJajdwfu3tDzvm5dwBbw3ZsDV/H2aeBVnf/18AbCbZJotrg7gfCbdAA/BugG9hAgtphZvXAXwCN7j4TGAXcTMK+F2Y2E/hz4EaC/WmJmV1LMrbFOqCp37TB6v0O4NrwcRvwuQrVcWTcPfIHMJHgw3sbsBkwgivYLg7L5wCPRV3PAtrxLFDXb9oB4Orw+dXAgajrOUT9RwOHCTvzk9iGPG16O/DDpLUDqAeOAFcR3LdiM7A4ad8L4N3AF3Je/zfgL5OyLYDJwL6c13nrDXweuCXf++L0iMsR/gMEO8HL4euxwCl3fyl8fZTgCxB3DnzXzPaY2W3htNe4+3GA8O+rI6vd8KYCJ4B/CNNrXzCzy0hWG/q7Gfh6+Dwx7XD3DuB/A78CjgOngT0k73uxD5hvZmPNrAZ4JzCJBG2Lfgard+8/6F6x3DaRB3wzWwK84O57cifneWsSzh+d5+7XE/y8+7CZzY+6QiN0MXA98Dl3nw38jnj+1C5ImN9+F/DNqOsyUmFu+CZgCjABuIxgv+ov1t8Ld3+KIA31ONAKPAG8NORMyZSImBV5wAfmAe8ys2eBhwjSOg8AtWbWewvGicCxaKpXOHc/Fv59gSBnfCPwvJldDRD+fSG6Gg7rKHDU3X8cvn6Y4B9AktqQ6x3AT939+fB1ktqxCDjs7ifcvQdYD8wlmd+LL7r79e4+H/g18AuStS1yDVbvowS/XHrFcttEHvDdfY27T3T3yQQ/v7e5+yrge8Cfhm+7FfhWRFUsiJldZmZX9D4nyB3vAx4hqD/EvB3u3gkcMbPp4aSFwM9JUBv6uYXz6RxIVjt+BbzZzGrMzDijpRK1AAAA3UlEQVS/LRL1vQAws1eHf68BVhBskyRti1yD1fsR4P3h2TpvBk73pn5iJepOhH4dJG8FNofPpwI/AZ4h+En+qqjrN0zdpxL8XH0C2A/cGU4fS9Ah/Yvw71VR13WYdjQAbcDPgI3AlUlrQ9iOGqALGJMzLVHtAD4JPE1w4PAV4FVJ+16E7fhngn9WTwALk7ItCP4xHQd6CI7gPzhYvQlSOv8X+CXwJMHZVZG3of9DQyuIiGRE5CkdERGpDAV8EZGMUMAXEckIBXwRkYxQwBcRyQgFfBGRjFDAFxHJiP8PNbRMSkzWC9wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c1d93df98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.title('Scatterplot for each method')\n",
    "plt.scatter(A,B)\n",
    "plt.show()\n",
    "\n",
    "def bland_altman_plot(data1, data2, *args, **kwargs):\n",
    "    data1 = np.asarray(data1)\n",
    "    data2 = np.asarray(data2)\n",
    "    mean = np.mean([data1, data2],axis=0)\n",
    "    diff = data1 - data2\n",
    "    md = np.mean(diff)\n",
    "    sd = np.std(diff,axis=0)\n",
    "    \n",
    "    plt.scatter(mean, diff, *args, **kwargs)\n",
    "    plt.axhline(md, color='gray', linestyle='--')\n",
    "    plt.axhline(md+1.96*sd, color='gray', linestyle='--')\n",
    "    plt.axhline(md-1.96*sd, color='gray', linestyle='--')\n",
    "\n",
    "bland_altman_plot(A, B)\n",
    "plt.title('Bland-Altman Plot of two methods')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) Most of the measurements are located between $\\mu-1.96\\sigma$ and $\\mu+1.96\\sigma$, and randomly distributed in both side of the mean value, which indicates that these two methods are agree with each other.\n",
    "\n",
    "(c) The are nearly one samples respectively lies above $\\mu+1.96\\sigma$, and below $\\mu-1.96\\sigma$, and it have large difference between each method, which can be considered as an evidence of heteroscedasticity or outlier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/regression/mixed_linear_model.py:2001: ConvergenceWarning: Gradient optimization failed.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/regression/mixed_linear_model.py:2039: ConvergenceWarning: The Hessian matrix at the estimated parameter values is not positive definite.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:1029: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return np.sqrt(np.diag(self.cov_params()))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               Mixed Linear Model Regression Results\n",
      "====================================================================\n",
      "Model:                 MixedLM     Dependent Variable:     values   \n",
      "No. Observations:      198         Method:                 REML     \n",
      "No. Groups:            2           Scale:                  13.2831  \n",
      "Min. group size:       99          Likelihood:             -548.6755\n",
      "Max. group size:       99          Converged:              No       \n",
      "Mean group size:       99.0                                         \n",
      "--------------------------------------------------------------------\n",
      "                          Coef.  Std.Err.   z    P>|z| [0.025 0.975]\n",
      "--------------------------------------------------------------------\n",
      "Intercept                 69.864    5.877 11.888 0.000 58.345 81.383\n",
      "subject                    0.490    1.862  0.263 0.792 -3.159  4.139\n",
      "Intercept RE              68.549                                    \n",
      "Intercept RE x subject RE  0.625                                    \n",
      "subject RE                 6.931                                    \n",
      "====================================================================\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+U3HV97/HnK8kCG6xsAoFDFkKiTYMChcBWQa5eJZYUr5KIIHDqNViuue21rVCbEq49Eq2UWNqjvcdWblpUrBaDgCGUWxEJ1pY22I0JkEBSUCDJJiQRWX+QLWw27/vH9zthMjuzMzvf2Z3vzL4e5+yZmc/3x7z3O8l7P/P5fL6fjyICMzNrX5OaHYCZmY0tJ3ozszbnRG9m1uac6M3M2pwTvZlZm3OiNzNrc070lmuSnpX0zgadq1PSvZJ+KukbjTjnWGjk75yVpBWSvtqgc10l6V8acS4bHSf6NiTpv0j61zSh/UTSw5J+LeM5h/0nlfRlSZ/OFm1jSHq7pJ1VdrsUOAE4NiIuG4ewWkqN19Ba0JRmB2CNJem1wD8AvwPcARwBvBV4uZlxlSNpSkQcGMe3PAX4j3reswmxmjWMa/Tt51cAIuL2iBiKiIGI+HZEPFbYQdKHJT0p6eeSnpB0dlq+XNIPi8rfm5a/AbgFOE/SLyT1S1oK/CbwR2nZvem+MyXdJWmfpGck/X7R+66QdKekr0r6GXBVUdnq9H1/IOnMcr+YpCMlfU7SrvTnc2nZ0cA/AjPTWH4haWbJsZ8EPgFcnm6/WtIkSX8s6TlJeyV9RdIx6f6zJUW633ZgXYWY3i1pU3pN/lXSrxZtK3s9q30OqbMkPZZ+K1st6agK739V+o3ts2kMP5L0lrR8R/p7LSm5hn8uabukPZJuSZu0RrqGR6TX5ueStkjqKTrfGyR9N33vLZIuLtp2rKS1kn4m6fvA68v9DjYOIsI/bfQDvBZ4AbgNuAiYVrL9MqAP+DVAwC8DpxRtm0lSAbgceAk4Md12FfAvJef6MvDpoteTgA0kCfUI4HXAj4CF6fYVwCCwON23s6jsUqAD+EPgGaAjPeZZ4J3p808B64HjgRnAvwJ/km57O7CzyrVZAXy16PVvAU+ncb4GuBv4u3TbbCCArwBHA51lznc2sBd4MzAZWJLGe2QN13Okz+FZ4PvpsdOBJ4HfrvA7XQUcAD6UxvBpYDvwV8CRwIXAz4HXpPt/DlibnveXgHuBmypdw/Sa/SfwrvT8NwHr020d6fX73+nnfUH6XvPS7V8n+VZ5NHB6+vv+y0ifkX/GKC80OwD/jMGHCm8gScI70ySwFjgh3XY/8NEaz7MJWJQ+v6r0PynDE/2bge0l+1wPfCl9vgL4Xsn2FYXEkb6eBOwG3pq+fpZXE/0PgXcV7bsQeDZ9PixJlfl9VnB4on8Q+F9Fr+eR/NGZwquJ/nUjnO8LpH9oisq2Af+1hutZ8XNIf+cPFL3+M+CWCvteBTxV9PqMNO4TispeAM4i+YPyEvD6om3nAc9UuobpNftO0es3AgPp87cCzwOTirbfnh4zOb2WpxZt+9PSf0P+GZ8ft9G3oYh4kiQBIOlU4KskNbkrgZNJEuYwkj4I/AFJkoOklnvcKN76FJKv/v1FZZOBfy56vaPMcYfKIuJg2iE4s8x+M4Hnil4/V2G/WpU73xSSDtthsZVxCrBE0u8VlR1RiKnK9az4OaSeL3q+n5F/zz1FzwcAIqK07DUk34KmAhskFbaJ5DMaSWksR0maksa0IyIOFm1/DuhO32sKh1+/4mtt48iJvs1FxFZJXwb+Z1q0gzJtpZJOAf4GWAD8W0QMSdpEkgggqSUOO33J6x0ktcO5I4VUpuzkojgmAScBu8rst4skuW5JX88q2q+eaVgL5yuYRfINaE8aQ7Xz7gBujIgbSzfUcD3Lfg5j7MckSf+0iOgrs32013AXcLKkSUXJfhbwH8A+kmt5MrC1aJs1gTtj24ykUyV9TNJJ6euTSWry69Nd/hb4Q0nnKPHLaVI6muQ/+r70uA+RtKsW7AFOknRESdnril5/H/iZpOvSDr7Jkk5X9aGd50i6JK0lXkMyQmh9mf1uB/5Y0gxJx5H0BRTGeO8Bji10ptboduBaSXMkvYakaWF11D665m+A35b05vRaHi3pv0n6Japfz0qfw5hJk/HfAJ+VdHwaV7ekhekuo72Gj5A0Bf2RpA5JbwfeA3w9IoZI+jxWSJoq6Y0kfRjWBE707efnJG3lj0h6iSRhbgY+BhAR3wBuBP4+3XcNMD0ingD+Avg3kv/wZwAPF513HUlN+nlJP07LbgXemI64WJP+534PSXvwMyQ1yL8FqiWOe0g6K18E/jtwSUQMltnv00Av8BjwOPCDtIyI2EqSuH+UxlNLk84Xgb8DvpfG+5/A7414RJGI6AU+DHw+jf1p0iazatez0udQ63tncF0a53olI5++Q9I3MeprGBGvABeTdPr/GPhr4IPpeQB+l6TJ6HmS/pwvNfy3sZoo7SQxawpJK4BfjogPNDsWs3blGr2ZWZtzojcza3NuujEza3Ou0ZuZtblcjKM/7rjjYvbs2c0Ow8yspWzYsOHHETGj2n65SPSzZ8+mt7e32WGYmbUUSTXdbeymGzOzNudEb2bW5pzozczanBO9mVmbc6I3M2tzuRh1Y2bWDtZs7OPm+7exq3+AmV2dLFs4j8Xzu5sdlhO9mVkjrNnYx/V3P87A4BAAff0DXH/34wBNT/ZuujEza4Cb7992KMkXDAwOcfP925oU0auc6M3MGmBX/8CoyseTE72ZWQPM7OocVfl4qproJX1R0l5Jm4vKbpa0VdJjkr4pqato2/WSnpa0rWiJMjOztrZs4Tw6Ow5fZ72zYzLLFs5rUkSvqqVG/2XgN0rKHgBOj4hfJVkI+HqAdF3IK4DT0mP+WlK1FebNzFre4vnd3HTJGXR3dSKgu6uTmy45o+kdsVDDqJuI+J6k2SVl3y56uR64NH2+iGRh4JeBZyQ9DbyJZN1MM7O2tnh+dy4Se6lGtNH/FvCP6fNuYEfRtp1p2TCSlkrqldS7b9++BoRhZmblZEr0kj4OHAC+Vigqs1vZJawiYlVE9EREz4wZVadTNjOzOtV9w5SkJcC7gQXx6nqEO4GTi3Y7CdhVf3hmZvmQ17tea1FXjV7SbwDXARdHxP6iTWuBKyQdKWkOMBf4fvYwzcyap3DXa1//AMGrd72u2djX7NBqUsvwyttJOlPnSdop6Wrg88AvAQ9I2iTpFoCI2ALcATwBfAv4SEQMVTi1mVlLyPNdr7WoZdTNlWWKbx1h/xuBG7MEZWaWJ3m+67UWntTMzKyKmV2d9JVJ6sd0dnD+ynW5b7f3FAhmZlWUu+u1Y5J46ZUDLdFu70RvZlZFubteX3PUFAaHDh89ntd2ezfdmJnVoPSu1znL7yu7X63t9uM5XNOJ3sysgpGScaV2+1pmqxzvRUrcdGNmVka1sfNZZqsc7+GaTvRmZmVUS8ZZZqsc7+GabroxMyujlmRc72yVWZp96uEavZlZGWO5YtR4L1LiGr2ZWaq48/WYzg46JuuwIZSNSsaFbwEedWNmNo5KR8L0DwzSMUlMm9pB//7Bhifj8VykxInezIzyna+DB4OpR0xh4ycubFJUjeE2ejMzWn/ispG4Rm9mE0qlm6DGeyTMeHKiN7Nca+RUASPdkbps4bzDtsHYjoQZT266MbPcavTKTiPdBJXlBqi8c43ezHKrUmL+2B2Pcu3qTaOu4Vdrhx/PkTDjyTV6M8utSol5KKKuGv5Y3gSVZ070ZpZbtSTg0UwGNt53pOZFLYuDf1HSXkmbi8ouk7RF0kFJPSX7Xy/paUnbJC0ci6DNbGIol5jLqXUIZDu3w4+kljb6LwOfB75SVLYZuAT4v8U7SnojcAVwGjAT+I6kX4mIwxvZzGxCqXfkTOlUAZMkhiKG7Teappd2bYcfSdVEHxHfkzS7pOxJAEmluy8Cvh4RLwPPSHoaeBPwb40I1sxaz2gX2Sj3R+Hh5ReUPRdMjKaXrBrdRt8N7Ch6vTMtM7MJajSLbFQbTjlRm16yavTwymFVfGD49yxA0lJgKcCsWbMaHIaZ5cVophaoNs4dJmbTS1aNTvQ7gZOLXp8E7Cq3Y0SsAlYB9PT0lP1jYGatbzRTCzRivpnxXHS7VTS66WYtcIWkIyXNAeYC32/we5hZCxnNkMas49wbfSdtu6hleOXtJJ2p8yTtlHS1pPdK2gmcB9wn6X6AiNgC3AE8AXwL+IhH3JhNbKNpV886zn28F91uFbWMurmywqZvVtj/RuDGLEGZWXuptV0968pL7TzVcBae68bMciVLZ2s7TzWchadAMLO2MVGnOKjGNXozaxvjveh2q3CiN7O24nH2w7npxsyszblGb2Zl+caj9uFEb2bDjHYiMss3N92Y2TC+8ai9ONGb2TC+8ai9ONGb2TATdW3VduVEb2bDVFrCb/8rBzJPELZmYx/nr1zHnOX3cf7KdRN+wrHx4M5YMxum0OG6Yu0W+gcGD5W/uH8wU6esO3mbwzV6Mytr8fxujj5yeF0wS6esO3mbw4nezCpqdKesO3mbw4nezCpqdKesO3mbw4nezCpq9GyQnl2yOdwZa2YVNXo2SM8u2RyKaP663D09PdHb29vsMMzMWoqkDRHRU20/N92YmbU5J3ozszZXtY1e0heBdwN7I+L0tGw6sBqYDTwLvD8iXpQk4C+BdwH7gasi4gdjE7qZjYc1G/sOu3Fq2tQObnjPaW5XbyG1dMZ+Gfg88JWisuXAgxGxUtLy9PV1wEXA3PTnzcAX0kczy7ly888DLPvGowwefLUv78X9gyy781HAd7O2iqqJPiK+J2l2SfEi4O3p89uA75Ik+kXAVyLp4V0vqUvSiRGxu1EBm1ljFCf2Yzo7eOmVAwwOJQm9MDXBUR2TDkvyBYNDwc33b3OibxH1Dq88oZC8I2K3pOPT8m5gR9F+O9OyYYle0lJgKcCsWbPqDMPM6lE650zxfDYFA4NDw6YrKOa7WVtHo8fRq0xZ2fGbEbEKWAXJ8MoGx2FmIyg358xoHdPZwfkr13k8fAuod9TNHkknAqSPe9PyncDJRfudBOyqPzwzGwu11sa7OjvomDS8/jZJ8NIrB+jrHyB4tanHUw7nU72Jfi2wJH2+BLinqPyDSpwL/NTt82b5U8vcMp0dk1lx8WncfNmZdHV2HCqfNrWDYzo7DrXnF3gWyvyqZXjl7SQdr8dJ2gncAKwE7pB0NbAduCzd/f+RDK18mmR45YfGIGYzq1Ghw7Wvf4DJEkMRdHd18o5TZ3DXhr7Dmm86JonXHDWF/v2Dw5piSptk5iy/r+z7ud0+n2oZdXNlhU0LyuwbwEeyBmVm1ZUbDlmckNds7GPZnY8eqnkPxasjau7a0Mf7zunmoa376mpjn9nVSV+ZpO5ZKPPJk5qZtZjSG5ig/EpNn7x3y7DmlYKBwSEe2rqPh5dfUFcMyxbOO2zUDngWyjxzojfLoUq19dJhkcUKbeSFRP/i/uFDJotlaWbxLJStxYneLGdGWle12rDI0STvAM5fua7uBL14frcTe4twojfLmZHWVS3XLl5sksQfr3mch7buq+m9vDj3xODZK81yplKtvK9/oOwdicWGIvjq+u1V/yAU87DI9udEb5YzlUauTJbK32beAB4W2d6c6M1yptK6qkNjuBqch0W2Nyd6s5xZPL+bmy45g+6uTgR0d3Ueej0WPCyy/bkz1iyHKo1oqTS0sl5eRGRicKI3y6lyY+lvuuSMQ2X1NORMljgY4XHvE4wTvdkYqDY9QS3Hl46lv2b1psNq4OevXFdxdE1XyUIikDTR3HTJGU7uE5Db6M0arJCks0zhW+nGqBf3Dx46V6VO289dfhabbriQmy89c1g7v5P8xOQavVmDjXTDU2mirVTzH2kcfOFchXlqKn1zKJ2moDBW3sl+4nGiN2uwSmPSS8srTXXQ+9xPEBWWZis510jTEIw0lYKT/cTiphuzBuua2lFTeaWa/1fXb6/a0VrLuPeRvlnYxOJEb9ZAazb20V9h1sjS+53qvRu11nHvtX6zsPbnphuzBik0lVSqjfcPDB62mHbX1I6qUwmX6h7FCB4vDmIFTvRmDfLJe7eMeDOT4FDiHc2kY1Df0EgvDmIFTvRmDbBmY9+ItfNqnasjHVfvzU1eHMQKnOjNGmCkDs7CotyjNW1qBxs/cWGWsLw4iAEZO2MlfVTSZklbJF2Tlk2X9ICkp9LHaY0J1Sy/Rurg/Iv3nznqCck6Josb3nNa1rDMgAyJXtLpwIeBNwFnAu+WNBdYDjwYEXOBB9PXZm2tUgdnV2cHi+d3l72LtZJJgpsvPdM1cWuYLE03bwDWR8R+AEn/BLwXWAS8Pd3nNuC7wHUZ3scsNyrdyVqp4/PdZ554aE6aaqtDQVKTd5K3RsuS6DcDN0o6FhgA3gX0AidExG6AiNgt6fhyB0taCiwFmDVrVoYwzMZHLXeaFv8ReMepM1j97zsOTSxWrZV+spzkbWzUnegj4klJnwEeAH4BPAocGMXxq4BVAD09PWO3dI5Zg1S707S0pv/Je7ccNnvkSDyzpI2lTJ2xEXFrRJwdEW8DfgI8BeyRdCJA+rg3e5hmzTfSot3lZqus5WYozyxp4yHT8EpJx0fEXkmzgEuA84A5wBJgZfp4T+YozXKg0p2mQNmafjXdXZ2HZqA0G0tZx9HflbbRDwIfiYgXJa0E7pB0NbAduCxrkGbjaTQdrvXyHao2njIl+oh4a5myF4AFWc5r1iy1driOdgqDYl6n1cabZ680K1Ktw3Xx/G4eXn5BTUMly/nc5Wex8RMXOsnbuHKiNytS69S+9cwA2d3V6QRvTeFEb1akUgIvLV+2cB4dk4bX6zsmiw+cO6vsWq5uk7dm8aRm1vIqdZ7Wc+w7Tp3BXRv6hjXf7H/lAGs29g1bj3XF2i30DyTDKIvb3ntOme5ZIy03FHXMqtdoPT090dvb2+wwrAWVdp5C7TcfVTr2fed08w+P7j6UwEd7XrPxImlDRPRU289NN9bSsqyLWunYh7bu4+gjh3/ZHRgc4prVm5j/qW+zZmNftsDNxpGbbqylZVkXtd5jX9w/yLI7HwVw7d5agmv01tJq7Twd7bHVjh8cipq+NZjlgRO9tbRy87zXOsKl3LEdk8T+Vw7UdENULd8azPLATTfW0rKsi1p67DGdHbz0yoGaJiOD+sbSmzWDR92YpQoLhNTCC4RYHnjUjdko1doUIy/1Zy3GTTc2YVS7sWqkaYgLPJbeWpFr9DYhFG6OKl0cpHg8fLUFvKdN7XCSt5bkGr1NCCPdWFU6rUFhGuLJEkMRdHsKA2txTvQ2IdR6c9Ti+d1O6NZ2nOitJWSZuAwqt797iKRNBG6jt9yrpX29miw3Vpm1OtfoLXdKa+/7XzlQtX293LHHdHYgQf/+QWZ2dfK+c7p5aOs+Tx1sE44TveVKuTVbKyltXy89tnia4b7+Ae7a0OdRMzYhZWq6kXStpC2SNku6XdJRkuZIekTSU5JWSzqiUcFa+ys3OqaS0vb1asfWOn2xWbupO9FL6gZ+H+iJiNOBycAVwGeAz0bEXOBF4OpGBGoTQ613p5ZrX88yNbFZO8vaGTsF6JQ0BZgK7AYuAO5Mt98GLM74HjaBdE3tKF/e2UF3VydKnx/VMYlrV2/i/JXrDnXKZpma2Kyd1Z3oI6IP+HNgO0mC/ymwAeiPiAPpbjuBsg2ikpZK6pXUu2/fvnrDsDayZmMfv/jPA8PKOyaLFRefxsPLL+Czl5/FywcO8uL+wUMjcK5dvYnZy+/jpZcP0DF5+ILdBR5lYxNV3Z2xkqYBi4A5QD/wDeCiMruWnR4zIlYBqyCZvbLeOKx93Hz/NgYPDv+nMDgUXLN6E5+8dwvAsHb4whH9A4N0TBLTpnbQv39w2Kgbj7KxiSrLqJt3As9ExD4ASXcDbwG6JE1Ja/UnAbuyh2kTQS1L+FUzeDCYesQUNn7iwkaFZdbysrTRbwfOlTRVkoAFwBPAQ8Cl6T5LgHuyhWgTRaPaz93hana4LG30j5B0uv4AeDw91yrgOuAPJD0NHAvc2oA4bQJ4x6kzGnIed7iaHS7TDVMRcQNwQ0nxj4A3ZTmvTUwPba2tU76rs4Ojj5xCX/8A4vBOIHe4mg3nO2MtN2ppcimMwCl0qmad7MxsInCit9yotsLTtKkd3PCe0w5L5J5W2Kw6J3obVyPVwJctnMc1qzeVPW6y5JE0ZnXyNMU2bqpNN7x4fjfnv3562WOvfPPJ4xipWXtxordxU2k5v4/d8eihZP+1D5/HB86dxWQld7hOlvjAubP49OIzxj1es3bhphsbN5U6W4ciuP7ux4GkVv/pxWc4sZs1kBO9NUQto19G6mwdGBxixdotHkFjNgbcdGOZ1brUX7nl/Ir1DwxmWi7QzMpTRPPnE+vp6Yne3t5mh2F1On/lurI19c6OSUw/+sjDaugAH7vjUYZq/HfX3dXJw8svaGi8Zu1C0oaI6Km2n2v0VtaajX2cv3Idc5bfd9ic7+VUansfGDw4rIYO8BfvP3PEmn0t5zaz2rmN3oYpt25rcWdpYZ9Ce/okqaYaemEpv0INvXQB8HKzU3reGrPsnOhtmErDIG++fxuL53cP+0NQazMMvFpDL72jtfSc4HlrzBrFid6GqdRcUigfzQLepSrV0AtJ36NuzBrPid6GqTQMspCk6203r1ZD97w1ZmPDnbE2TLlhkMVJulq7uUhGy3zg3FmHFvTu7urkpkvOcCI3awLX6G2Y0maUwtqr167exM33b+Mdp87ga+u3l10M2MMhzfLHNXora/H8bh5efgGfvfwsXj5wkBf3Dx4aJnnXhj7e8vrpqOQYd56a5ZMTvY2o0gicZ18Y4LOXn+WmGbMW4KYbG9FII3DceWrWGpzobcQJyaqNwDGz/Ku76UbSPEmbin5+JukaSdMlPSDpqfRxWiMDtsaqNiFZtRE4ZpZ/dSf6iNgWEWdFxFnAOcB+4JvAcuDBiJgLPJi+tpwa6S5YSDplb7rkDLfFm7WwRjXdLAB+GBHPSVoEvD0tvw34LnBdg97HGqxSG3xxc43b4s1aW6NG3VwB3J4+PyEidgOkj8eXO0DSUkm9knr37dvXoDBsNNZs7GOSSgdJJpRuN7PWlznRSzoCuBj4xmiOi4hVEdETET0zZszIGoaNUqFtvtKEZAGHmm/MrLU1okZ/EfCDiNiTvt4j6USA9HFvA97DGqyWick8F7xZe2hEor+SV5ttANYCS9LnS4B7GvAe1mC1JHEPoTRrD5kSvaSpwK8DdxcVrwR+XdJT6baVWd7Dxka1JO4hlGbtI1Oij4j9EXFsRPy0qOyFiFgQEXPTx59kD9MabaSFuj2E0qy9eK6bCWrx/G7ed053xYnJnOTN2ocT/QT20NZ9w6YaLr5Zyszag+e6aQMjzVUzkmpLBppZe3CNvsVVm6tmJJU6ZD3axqy9ONG3uGpz1YzEE5aZTQxuumlxWZpfSpcMHE2zj5m1Dif6FldpvvhJEms29lVN2p6wzKz9uemmxVUaDz8UUXNbvZm1Nyf6FleYL35ymVkoPVTSzMCJvi0snt/NwQqzUHqopJk50bcJD5U0s0qc6BtgzcY+zl+5jjnL7+P8leua0i7uoZJmVolH3WRUuGGpMJa9cMMSMK6jWTxU0swqcaLPaKQblsY7yXqopJmV46abjDxfjJnlnRN9Ru4ENbO8c6LPyJ2gZpZ3bqPPyJ2gZpZ3TvQN4E5QM8uzrIuDd0m6U9JWSU9KOk/SdEkPSHoqfZzWqGDNzGz0srbR/yXwrYg4FTgTeBJYDjwYEXOBB9PXZmbWJHUnekmvBd4G3AoQEa9ERD+wCLgt3e02YHHWIM3MrH5Z2uhfB+wDviTpTGAD8FHghIjYDRARuyUdX+5gSUuBpQCzZs3KEMbYqnc9VjOzvMjSdDMFOBv4QkTMB15iFM00EbEqInoiomfGjBkZwhg7WdZjNTPLiyyJfiewMyIeSV/fSZL490g6ESB93JstxObJsh6rmVle1J3oI+J5YIekwp1BC4AngLXAkrRsCXBPpgibyNMbmFk7yDqO/veAr0k6AvgR8CGSPx53SLoa2A5clvE9mqbSeqye3sDMWkmmRB8Rm4CeMpsWZDlvXixbOO+wKYjB0xuYWevxnbEj8PQGZtYOnOir8PQGZtbqPHulmVmbc6I3M2tzTvRmZm3Oid7MrM050ZuZtTknejOzNudEb2bW5pzozczanBO9mVmbc6I3M2tzLTsFgld+MjOrTUsm+sLKT4VZJQsrPwFO9mZmJVqy6cYrP5mZ1a4lE71XfjIzq11LJvpKKzx55Sczs+FaMtEvWziPzo7Jh5V55Sczs/JasjPWKz+ZmdWuJRM9eOUnM7NaZUr0kp4Ffg4MAQciokfSdGA1MBt4Fnh/RLyYLUwzM6tXI9ro3xERZ0VET/p6OfBgRMwFHkxfm5lZk4xFZ+wi4Lb0+W3A4jF4DzMzq1HWRB/AtyVtkLQ0LTshInYDpI/HlztQ0lJJvZJ69+3blzEMMzOrJGtn7PkRsUvS8cADkrbWemBErAJWAfT09ETGOMzMrIJMiT4idqWPeyV9E3gTsEfSiRGxW9KJwN5q59mwYcOPJT2XJZZROg748Ti+32g5vvrlOTZwfFk5vsOdUstOiqivMi3paGBSRPw8ff4A8ClgAfBCRKyUtByYHhF/VNebjBFJvUWdx7nj+OqX59jA8WXl+OqTpUZ/AvBNSYXz/H1EfEvSvwN3SLoa2A5clj1MMzOrV92JPiJ+BJxZpvwFklq9mZnlQEvOddMAq5odQBWOr355jg0cX1aOrw51t9GbmVlrmKg1ejOzCcOJ3syszbV9opfUJelOSVslPSnpPEnTJT0g6an0cVqTYpsnaVPRz88kXZOX+NIYr5W0RdJmSbdLOkrSHEmPpPGtlnREE+P7aBrbFknXpGVNu36Svihpr6TNRWVl41Hi/0h6WtJjks5uUnyXpdfvoKSekv2vT+PbJmlhk+K7Of3/+5ikb0rqyll8f5LGtknStyXNTMvH/fOtKCLa+odkvp3/kT5lzWmuAAADgklEQVQ/AugC/gxYnpYtBz6TgzgnA8+T3ACRi/iAbuAZoDN9fQdwVfp4RVp2C/A7TYrvdGAzMJVkBNl3gLnNvH7A24Czgc1FZWXjAd4F/CMg4FzgkSbF9wZgHvBdoKeo/I3Ao8CRwBzgh8DkJsR3ITAlff6ZouuXl/heW/T894FbmvX5Vvpp6xq9pNeSfDC3AkTEKxHRTz4nXlsA/DAiniNf8U0BOiVNIUmou4ELgDvT7c2M7w3A+ojYHxEHgH8C3ksTr19EfA/4SUlxpXgWAV+JxHqgK72bfFzji4gnI2Jbmd0XAV+PiJcj4hngaZK738c7vm+nny/AeuCknMX3s6KXR5PMAVaIb1w/30raOtEDrwP2AV+StFHS36Z38dY08do4uwK4PX2ei/giog/4c5Ib33YDPwU2AP1F//F2ktT8m2Ez8DZJx0qaSlKDOpmcXL8ileLpBnYU7dfMa1lOHuP7LZJaMuQoPkk3StoB/CbwibQ4N/G1e6KfQvI16wsRMR94iRzOj5+2cV8MfKPZsRRL25IXkXwtnklSW7mozK5NGaMbEU+SfJV/APgWydf4AyMelC8qU5an8c65ik/Sx0k+368Visrs1qx/ix+PiJNJYvvdtDg38bV7ot8J7IyIR9LXd5Ik/j2Fr1C1Trw2xi4CfhARe9LXeYnvncAzEbEvIgaBu4G3kHwFLdxVfRKwq0nxERG3RsTZEfE2kq/UT5Gf61dQKZ6dJN9ACpp6LcvITXySlgDvBn4z0gZwchRfkb8H3pc+z018bZ3oI+J5YIekeWnRAuAJYC2wJC1bAtzThPCKXcmrzTaQn/i2A+dKmqpkUqPC9XsIuDQH8aFkimwkzQIuIbmOebl+BZXiWQt8MB2dcS7w00ITT06sBa6QdKSkOSQd3d8f7yAk/QZwHXBxROzPYXxzi15eDBSma8/P59usXuDx+gHOAnqBx4A1wDTgWJJlDp9KH6c3Mb6pwAvAMUVleYrvkyT/cDcDf0cywuF1JP+hniZpbjqyifH9M8kfn0eBBc2+fiR/aHYDgyQ1uqsrxUPy1f6vSEaLPE7RiJdxju+96fOXgT3A/UX7fzyNbxtwUZPie5qkrXtT+nNLzuK7K/3/8RhwL9DdrM+30o+nQDAza3Nt3XRjZmZO9GZmbc+J3syszTnRm5m1OSd6M7M250RvZtbmnOjNzNrc/weZJdndiRQe2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c1d8eca20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X2cVOV99/HPT1hwF4VVdg1hgQAmYAXjkqyJYooPkEASVORuiYY2atLQp7tNclsSrb6iaW3jXdImaXrfaX0llqShJuiN1IcGVDTQRkxEwQhRNEqUB1FAF6KLusDv/uOcgdmZOTOzM2d3zpz9vl+vfbFzrvNwzdnhN9f5Xde5jrk7IiJS/46rdQVERCQeCugiIimhgC4ikhIK6CIiKaGALiKSEgroIiIpoYCeMGa21Mxu6qN9/9rMZtVq+1rpz3qb2U1mttfMdvfH8Wohzs+omd1oZj+IY1+igN7vwuBy0MxeN7PXzOxeMxtb63plmNmwsG7/WWK9RP1HNDM3szfCuu80s38ws0G93Mf5ZrajijqMBa4GTnf3UXHvvxbM7Eoz++9a10PKo4BeGxe5+wnAO4GXgW/VuD7Zfgd4C/iImb2z1pXppTPD8zoT+CTw2X4+/ruAfe7+Sj8fVwRQQK8pd38TuAM4vVC5mZ1kZveY2Z6wNX+PmY3JKv+Jmf21mf3UzH5jZveZWUtW+e+b2Qtmts/MriuzWlcA/wz8AlgYUa85wF8CnwhbxE9k1ecmM3s4XH63mY00s2VmdsDMHjWz8Vn7+aaZbQ/LHjOz384qu9HMlpvZ98P3tsXMOsp5A+7+NPBfwNQCdR9qZt8ws13hzzfCZcOAHwOjw7q/bmajC2w/IqzTnvDcXm9mx4Upnfuztl+as13B/YdXay3hOteb2SEzGx6+vsnMvlHsuBF/nxvN7HYz+0F47p40s0lmdq2ZvRKe84/kvKfvmtlL4dXNTWY2yMx+i+CzcE5Y386sw5wUXl3+xsx+ZmanZu1vevi33h/+Oz2rbIKZrQ23ux/I/rweH9Z5n5l1htu+I+LPLAUooNeQmTUBnwAeiVjlOOBfCVp+44CDwD/lrPNJ4CrgFGAI8Bfhvk8Hvg38PjAaGAmMoQgzGwecDywLfz5VaD13XwX8LfAjdz/B3c/MKr4sPGYbcCqwPnwPJwNPATdkrfso0B6W/Ttwu5kdn1V+MfBDoBm4q8B7j3ofpwO/DWwsUHwdcHZ43DOBDwDXu/sbwEeBXeF7OsHddxXY/lvACGAicB7BObrK3R/I2f7K7I2K7P/RcD8AM4AXgHOzXq8tdtwip+Ei4N+Ak8LzsJrg89QG/BXwL1nrfg84BLwbmAZ8BPgDd38K+CNgfVjf5qxtLge+Eu7/V8DfAJjZycC9wD8SfOb+AbjXzEaG2/078BhBIP9rggZExhXhexwbbvtHBJ95KZe766cff4BfA68DnQT/iXYBZ2SVLwVuiti2HXgt6/VPCIJR5vWfAKvC378M/DCrbBjwNjCrSN2uBzaFv48GDgPTcuo+K/z9RuAHOdv/BLgu6/XfAz/Oen1RZv8Rx3+NIG2S2f8DWWWnAweLbOvAgXAfzwE3AccVqPdzwMeytpsN/Dr8/XxgR5FjDCJIR52etewPgZ+UuX1eOUFQ+0dgMLAb+BxwM3A8QTBrKXXcAse5Ebg/57y/DgwKX58Ynq9m4B3hvhuz1r8ceCj8/Urgv3P2vxT4TtbrjwFPh7//PvDznPXXh/sZR/CZH5ZV9u+ZzxHwaeBh4L21/n9arz9qodfGPA9aO0OB/wmsNbNCnWhNZvYv4SX2AWAd0Gw9O/uyR1N0ASeEv48GtmcKPGgh7sva9+tZP+PCxZ8iaJnjQetxLT1bUOV4Oev3gwVeZ+qHmV1tZk+Fl+adBK2zlqz1c9/b8WY2uMix3+fuJ7n7qe5+vbsfKbDOaIJWcMYL4bJytBBcBeVu31bm9oWsJQj07wOeJEjbnEdwFfErd99b4XFzz/tedz+c9RqCv8W7gAbgpTDN0UnQej+lRL2Lfe5eyFk3U9fRBA2SN3LKMv6N4Erih2E67O/MrKFEPSSLAnoNufthd19B0BL+UIFVrgYmAx909+EEl+AAVsbuXyK4dA02CNI7mcte/Nhl/wnu/mKY53wPcK2Z7bZg2N0HgcsjgmhV03SG+fIvAQuAk8IvuP2U996qsYsgiGWMC5dB6fe0F+gusP3OMo9daP8PE/yNLwXWuvsvw31+nGPplmqPW8x2ghZ6i7s3hz/D3X1KkToXk3t+4VhdXyLIvQ/LKQsO5N7t7l9x99OB6cBcItJ+UpgCeg1Z4BKCPORTBVY5kaA11RnmJm8osE6UO4C5ZvYhMxtCkDct9ve+gqB1eDpBaqedoFOxiSD3m+tlYHxUx1wZTiS4/N4DDDazLwPDK9xXb9wGXG9mrWFn5JeBzPDLl4GRZjai0IZhC3c58DdmdqKZvQv4X1nbl5K3f3fvIsgp/ynHAvjDBCmVtTEdN5K7vwTcB/y9mQ0PO3hPNbNMXv9lYEz4GSrHfwKTzOyTZjbYzD5B8Jm6x91fADYAXzGzIWb2IYJ0EABmdoGZnRFegR4g+BI7XOAYEkEBvTbuNrPXCT60fwNc4e5bCqz3DaCRoIX2CLCq3AOE+/tTghzlSwS55YJjoMOOyAXAt9x9d9bPNoLL4EJpl9vDf/eZ2ePl1ivLaoJRH88QXHa/SVaKqA/dRBBUfkGQ4ng8XIYHo2NuA54P0w+FUjF/BrwBPA/8N8H5vbWcAxfZ/1qCtMfPs16fSJBiq/q4ZfgUQUrnlwSfkzsIhtQCPAhsAXab2d5SO3L3fQQt66sJUnxfBOaGqSMIOvE/CLxK0ED5ftbmo8JjHyBo4Kwlhi+tgcTCzggREalzaqGLiKSEArqISEoooIuIpETJgG5mt4a3C2/OWtZuZo+Y2SYz22BmH+jbaoqISCklO0XNbAbBXWbfd/ep4bL7gK+7+4/N7GPAF939/FIHa2lp8fHjx1ddaRGRgeSxxx7b6+6tpdYrdtcdAO6+zrImVMos5tiY4REcuzGjqPHjx7Nhw4ZyVhURkZCZ5d59W1DJgB7h88BqM/saQdpmetSKZrYIWAQwbty4qNVERKRKlXaK/jHwBXcfC3wB+G7Uiu5+i7t3uHtHa2vJKwYREalQpQH9CmBF+PvtBFOQiohIDVUa0HdxbA7nC4Fn46mOiIhUqmQO3cxuI5jes8WC5yHeQPBor2+Gs/C9SZgjFxGR2ilnlMvlEUXvj7kuIiJ1aeXGnSxZvZVdnQcZ3dzI4tmTmTetmmnyK1PpKBcRESEI5teueJKD3cFMvzs7D3LtiicB+j2o69Z/EZEqLFm99WgwzzjYfZglq7f2e10U0EVEqrCrs/BzrKOW9yUFdBGRKoxubuzV8r6kgC4iUoXFsyfT2DCox7LGhkEsnj253+uiTlERkSpkOj41ykVEJAXmTWurSQDPpZSLiEhKKKCLiKSEArqISEoooIuIpIQCuohISiigi4ikhAK6iEhKKKCLiKSEArqISEoooIuIpIQCuohISpQM6GZ2q5m9Ymabc5b/mZltNbMtZvZ3fVdFEREpRzkt9KXAnOwFZnYBcAnwXnefAnwt/qqJiEhvlAzo7r4OeDVn8R8DN7v7W+E6r/RB3UREpBcqzaFPAn7bzH5mZmvN7KyoFc1skZltMLMNe/bsqfBwIiJSSqUBfTBwEnA2sBhYbmZWaEV3v8XdO9y9o7W1tcLDiYhIKZUG9B3ACg/8HDgCtMRXLRER6a1KA/pK4EIAM5sEDAH2xlUpERHpvZKPoDOz24DzgRYz2wHcANwK3BoOZXwbuMLdvS8rKiIixZUM6O5+eUTR78VcFxERqYLuFBURSQkFdBGRlFBAFxFJCQV0EZGUUEAXEUkJBXQRkZRQQBcRSQkFdBGRlFBAFxFJCQV0EZGUUEAXEUmJknO5iEjtrdy4kyWrt7Kr8yCjmxtZPHsy86a11bpakjAK6CIJt3LjTq5d8SQHuw8DsLPzINeueBJAQV16UEAXSbglq7ceDeYZB7sPs2T11lgDuq4C6p8CukjC7eo82KvlldBVQDqoU1Qk4UY3N/ZqeSWKXQUkycqNOzn35geZcM29nHvzg6zcuLPWVUoUBXSRhFs8ezKNDYN6LGtsGMTi2ZNjO0Z/XAVUK3MVsbPzIM6xqwgF9WNKBnQzu9XMXgkfN5db9hdm5mamB0SL9JF509r46vwzaGtuxIC25ka+Ov+MWFMh/XEVUK16uYqopXJy6EuBfwK+n73QzMYCHwZejL9aIpJt3rS2Ps1lL549uUcOHeK/CqhWPVxF1FrJFrq7rwNeLVD0deCLgB4OLVLn+uMqoFr1cBVRaxWNcjGzi4Gd7v6EmZVadxGwCKCtrY2lS5f2KJ8yZQpnnXUW3d3dLFu2LG/79vZ22tvb6erqYvny5XnlHR0dTJ06lf3793PnnXfmlZ9zzjlMnjyZvXv3cs899+SVz5gxg4kTJ7J7925WrVqVVz5z5kzGjh3L9u3bWbNmTV75nDlzGDVqFM8//zzr1q3LK587dy4tLS1s3bqV9evX55VfeumljBgxgs2bN7Nhw4a88gULFtDU1MSmTZvYtGlTXvnChQtpaGjg0UcfZcuWLXnlV155JQAPP/wwzzzzTI+yhoYGFi5cCMDatWvZtm1bj/KmpiYWLFgAwAMPPMCOHTt6lA8fPpz58+cDsGrVKnbv3t2jfOTIkVx00UUA3H333ezbt69H+ahRo5gzZw4AK1as4MCBAz3Kx4wZw6xZswBYvnw5XV1dPconTJjAeeedB8CyZcvo7u7uUT5p0iSmT58OkPe5A332Cn32PjsKGBWUf2Ty2QCJ+exdesJbbHvzDQ4fcbp8COu6J9LYMIhPtr2a9/dN+2cvSq87Rc2sCbgO+HI567v7Le7e4e4dDQ0NvT2ciAgALScMZULLMIYODjqIM1cRp7aeUOOaJYe5l86YmNl44B53n2pmZwBrgMxX1hhgF/ABd99deA+Bjo4OL9QKFRGRaGb2mLt3lFqv1ykXd38SOCXrQL8GOtx9b2/3JSIi8Sln2OJtwHpgspntMLPP9H21RESkt0q20N398hLl42OrjYiIVEx3ioqIpIQCuohISmi2RRGJjabgrS0FdBGJhabgrT2lXEQkFpo8q/YU0EUkFpo8q/YU0EUkFpo8q/YU0EUkFv3xII4kSPJTk9QpKiKxyHR8pnmUS9I7fhXQRSQ2ff0gjlor1vGbhPetgC4iA1IlY+aT3vGrHLqIDDiVPnA66R2/CugiMuBUOmY+6R2/SrmIpIRuuy9fpamTpHf8KqCLpEDSR1/UUqEvutHNjewsELzLSZ0kueNXKReRFNBt94VF5covOK010amTSimgi6RA0kdf1ErUF91DT+/hq/PPoK25EePYA6eT2vIul1IuIinQ3NTAa13decuTMvqiVop90SU5dVKpkgHdzG4F5gKvuPvUcNkS4CLgbeA54Cp37+zLiopIYSs37uT1Nw/lLW8YZL1KIaSxU7XSXHm9notyUi5LgTk5y+4Hprr7e4FngGtjrpeIlGnJ6q10H/G85cOGDC47CFU6LjvpKhlmWM/nomRAd/d1wKs5y+5z90yT4BFgTB/UTUTKEJVW2H8wPwUTJa2dqvOmtfU6V17P5yKOHPqngR9FFZrZImARwLhx42I4nIjAsbRAfts80Jv8eZo7VXubK6/nc1HVKBczuw44BCyLWsfdb3H3DnfvaG1treZwIhLKTgsU0tsheEm/pb0/1fO5qDigm9kVBJ2lC909qpEgIn2gUFogo5IheEm/pb0/1fO5qCjlYmZzgC8B57l7V7xVEpFSoi7/DfjpNRf2en9Jv6W9P9XzuShn2OJtwPlAi5ntAG4gGNUyFLjfzAAecfc/6sN6ikiWam5dl9LqdYx6yYDu7pcXWPzdPqiLpES9juGtJ4tnT+4xdwtUlxbQXDDpoFv/JVb1PIa3nlQyHK+Yeh6qJ8fo1n+JVdIf0ZUmcaYF6nmonhyjFrrESoGhPtXzUD05RgFdYqXAUJ/qeaieHKOALrFKemBYuXEn5978IBOuuZdzb35Quf1Q3Dl5qQ3l0CVWSR7Dq5EcxdXrUD05RgFdYpfUwBDVYXv18ieA8oN6UoZlJqUekhwK6DJgRHXMHnYvu6WelFZ+UuohyWL9OQ1LR0eHb9iwod+OJ5Lt3JsfjJzMKqOtREs3ah/NjQ0MGzq431rLUfVoa26s6NZ/STYze8zdO0qtpxa6DBiF7q7MVaqlG9XK7zzYTWc4/3hmHxteeJWHnt4Te5BfuXFn5BeThocObBrlIgNGZiTHoGD+oUjF7pAsd/jlwe7DLHvkxdjvmM2kWqJoeOjApoAudaeaoYfzprXx9wvOzBtamSuqpVtoWGaU3GRmHLfSF5s2N0nDQ6U2lHKRxCk2eiOOzsDsoZVRqYuolm6hYZldbx/ita7yHvdWTUqkWKoFODpuXKNfBi4FdEmUUgE7rrliMusuvv2JvAcsNwyyoi3d7GGZKzfu5Ma7tuStY+S30AmXnXvzgwWDbDlfZFHamhuPBvNyvvAU9NNJKRdJlFKz/sU5V8yS1VvzgjnAsCGDywpumeDZmfMw5pOaGlh49rjI1EyhfHqpWSrLTbWUM2uiZsRMLwX0lEjLLe2lAna1c8Vkn6eo9MX+g+WlT6KCbNOQwdw074yjt9IXkhtkK/0iA3rcol/OF56myk0vBfQUSFOLq1TArmaumNzz1Ns65CoVPOdNa+On11xI1JiaXZ0Hj37BlBqGGFWnTKoFgvd3XMQInuztNSNm5ZLecCoZ0M3sVjN7xcw2Zy072czuN7Nnw39P6ttqSjFpanGVCtjVTCJVLG1R6FillHu10NzUUHA9Bxbf8UTRjs5yv8gyX1aHC9womPueNCNmZeqh4VROC30pMCdn2TXAGnd/D7AmfC01kqYWVzkBO9Py3Xbzx/npNReW3ZlX7HxUMsNguVcLxW7G7j4cXdibL7KoL6tBZnnvKekzYiZVPTScynmm6DozG5+z+BKCB0cDfA/4CfClGOslvZC2Bwb31eReUeep0tvlM3W88a4tRztGj2/IbyOVm5PPrVPuyJPcIZOZQDJvWlvkl9UR97xzmeQZMZOsHhpOlQ5bfIe7vwTg7i+Z2SlRK5rZImARwLhx4yo8nBQT9wOD06qvztNbh44c/f21ru68YYLNTQ1lj1OH6C+YYkMSo76sRjQWTvckdUbMJKuHhlOfd4q6+y3u3uHuHa2trX19uAFJDycoT1+cp1KX4Ss37uT1Nw+Vvb9iXzDFjrV49mQajsvvEH3j7UOJyvHWs3pIVVXaQn/ZzN4Zts7fCbwSZ6Wk99TiKk/c56nUZXjUWPdsFt6FVCr1UexY86a18ZW7t+RdCXQfdj2gOyb1kKqqNKDfBVwB3Bz++x+x1UgkQhLvbix1GV5WftVh280fr/pYnRFpnew6JPEc1pOkN5zKGbZ4G7AemGxmO8zsMwSB/MNm9izw4fC1SJ+N0631kLGo91XqMryc/Gq5OdhKj5VZXutzKH2vnFEul0cUzYy5LlLn+vIpOnHN4VKJ61c+ybJHXjx6M1Kh9xXV6i01B3tvcrCVHKvcaQGS1uosdSWhK43CNDmXxKYvA0athoyt3LizRzDPyH5fxS7Dc4PwiMYGzIL0SCWBqDfHyt1/PQy7g9INAz1+L5oCekLVYwskroBR6L2XO2SskvNWbJslq7dGThNQ7vvqz7xrsWPVw7A7KN0wqKcrjf6muVwSqF5znXHcUh713i84rbXkkLFKzlupbYoF7cxUuHH2E/TlPCH1MOwOSjcM6uVKoxb69SHREyZM8BtuuKHHsilTpnDWWWfR3d3NsmXL8rZpb2+nvb2drq4uli9fnlfe0dHB1KlT2b9/P3feeWde+TnnnMPkyZPZu3cv99xzT175jBkzmDhxIrt372bVqlV55TNnzmTs2LFs376dNWvW5JXPmTOHUaNG8fzzz7Nu3bq88rlz59LS0sLWrVtZv359Xvmll17KiBEj2Lx5M5kHaG98sZO3DgUtkIfePpW3aODdg/YyZehrTBvX3GP7hQsX0tDQwKOPPsqWLfnzcl955ZUAPPzwwzzzzDM9yhoaGli4cCEAa9euZdu2bT3Km5qaWLBgAQAPPPAAO3bs6FE+fPhw5s+fD8CqVavY/KsX2Lb3DQ6Hw/QO+PFs5NTgsW87Hmffvn3sff0ttr96kLcPHeatIcO5+OMfY960NlasWMGBAwd6vPc9R07gsUNjAPjYsG1MPKnh6LZDBg/i9Env5g8vmwvAsmXL+Plze45uC7DjSDObD42irbmRz456Me/cTJkyhc+v+Q27O1/nw0Oe7VE2dPAgrrr4fP501avs7TzABUOey9t+6+FT2Hb4ZEY2HOLyd+ym5YShPcp789n7wR3/0ePcAWz2cfzF/Om8v+VIbJ+97PO/tfF0/nzOe3l3w2sUenj7ggULaGpqYtOmTWzatCmvvK8+extf7OTAIeOht98NwPsH76D1uNcZOngQ08Y1s/HFTl7rHsS67okAfKDhRU62rqPlACNHjuSiiy4C4O6772bfvn09jj9q1CjmzAlmNMl89rKNGTOGWbNmAbB8+XK6urp6lE+YMIHzzjsPCD573d09RxhNmjSJ6dOnA7B06dK8c9PbuHfVVVeV9ZBotdAT6O1DhTvQopYnRcsJQ5nQMoyhgwdhBPOKZ9+4s/f1t9i29w3eOnQYB15/61BeCzrqPR7sPkzLCUOZNq6ZD04cybRxzZz2zhN7rBO1bbGWW1RZZl+LZ0/m+BKPnHvz0BG2v1pd63D7qwd7BHOAtw8fiX2ekOxzeOefTE9kimLsyY15z30ddJwx9uTGY+XHRZcPZP3aQu/o6PBCLQHpKWo61UrnHOkLleSqo97XSU0NbPzyR4quU857r2TbcrbJfq9R/1uM8saSR5lwzb0F913tfuuVRrn0ZGZltdDVKZpASZ+bpdJRBlGt4de6ulm5cSfzprVV9d4r2bacbbI7GqO+AKrtWKyXDsv+UqojOek3+NSKUi4JlPS5WSqdRrRYcLp6+RNMuOZelqzeyv94f1tF772S89bbbXrTsdibTs566bCU3unvB2Io5SK9Vml6YOXGnXz+R/mda7kaGwYl6gssVzmX+7lXMVD6fQ20NELaVfIZiFJuykUBXXqtmjx3+1fuy3uociFJ6i+oRD30g0jfivMzUG5AV8ol5frikq+a9MCNF0/J27aQeh9THPdY6aQ/y1Ly1WK8vDpFU6yvbpGudhrRoYOPO1qncObYPPXeGVisk7O3qRXd6l6fatHRrRZ6ivXlMxArea5nJjBlp1wKBfOGQVb3nYFRVzEXnNba67tZ6+FZlpKvFh3dCugpVotLvmKpgagHGecaNmRw3bc8o0bPPPT0nl4HZ93qXp9qMVpNKZcU6+9LvlKpgXIDUCUPVU6iQmOlvxAxyqfYudEY9frV3+Pl1UJPsf6+5CuVGig3AKU5UFUygZnGqEu5FNBTrL8v+UqlBgoFplxpD1SVBOek32gmyaGUS0rljqT4+ifaex0ACo3GgOjRLaVSA4VGx1xwWisPPb1nwNxMU+kIoVrf6q6bnupDVTcWmdkXgD8gGKzwJHCVu78Ztb5uLOofcdyhVmgfDccZWPAk+UL7jfPOOEkO/V1rr89vLDKzNuDPgQ53nwoMAi6rdH8SnziGuRXaR/cR7xHMc/er1EA6adhk/ag25TIYaDSzbqAJ2FV9laRacQxzq3TdWqcGJH4aNlk/Km6hu/tO4GvAi8BLwH53vy93PTNbZGYbzGzDnj17Kq+plC2OR8H11boSj/6cCiCOz5P0j2pSLicBlwATgNHAMDP7vdz13P0Wd+9w947W1tbKaypl6+1IikLBodA+Go4zGgb1fFJM2kelJNH1K5/kCz/a1G/PnNWwyfpRzbDFWcA2d9/j7t3ACmB6PNWSavQmlx31kGQgbx9LfvdMlvzOmcqR19DKjTtZ9siLeVMm9GVOO2l9I5qoLFrFo1zM7IPArcBZwEFgKbDB3b8VtY1GuSSPpnmtL1F/LxgYj6sbqCNu+vwRdO7+MzO7A3gcOARsBG6pdH9SG8U6vNI69rie31epKQLSrtiIm3r5G/alqka5uPsNwA0x1UUKiCP4FNtH1M1AIxobUjlla71PRRv19zIYEDltjbgpTrf+J1hUfrs3OcNS+4jq8DIjlWOP631MdaG/lwELzx5XF19I1dKIm+IU0BOsr24QKudmoM6uwjMe7uw8WNcdUpW28JLSEVfo7/X1T7Rz07wzalKf/qYRN8VpLpcE68sbhErdDLRk9dbIS/vM8npLV0BlU9EmLU0zkG/eqvZpWWmnFnqC9eUNQqX2EXVpH+dwuVq0eitp4dV7miZtKnla1kChgJ5gcVxeVrqPQpf2UQNcK+mQiqN/oBLFxlRHfcGUe6WUlLSMDFxKuSRYHJeX1ewj99I+agx0JR1StRx+VihlUSytUk6aJmlpGRmYFNATLo58aVw518WzJxe8qaOSDqmkDT8r9gVTzvvW+GhJAgV0KVucHVJJe05msS+Yct530r6gZGBSQJdeSWJrPw7lPG2p2PsudoOWSH9Rp6jURNImfKq2A3rx7MnBE51yvPH2IXWOSr+p6hF0vaXJuSTJqp1mYdpf3cdrBW7I0kRnUq0+n5xLJG2qTSdF3V2rPLr0F6VcRGKieUak1hTQJdHq6WadqLtrM/PfJLnukg5KuUhi1dvNOtnDG3d2HuwxVULS6y7poBa6JFY9zqGSmWek0FQJSa+71D8FdEmser5Zp57rLvWrqpSLmTUD3wGmElxdftrd18dRMZGk3U3aG3HVvZ4fl1eugfAe+0u1LfRvAqvc/TTgTOCp6qskEqjnhxnEUfdazUjZnwbCe+xPFQd0MxsOzAC+C+Dub7t7Z1wVE0na3aS9EUfdk9qHEOfIo6S+x3pVTcplIrAH+FczOxN4DPicu7+RvZKZLQIWAYwbN66Kw0naRV1610MAL6TauicxDx/3yKMkvsd6Vk3KZTDwPuDb7j4NeAO4Jncld7/F3TvcvaO1tbWKw0ma6dJo+P1eAAAJaklEQVQ7XxJvVIq7RZ3E91jPqgnoO4Ad7v6z8PUdBAFepNd06Z0viX0Icbeok/ge61nFAd3ddwPbzSxz5mcCv4ylVjLg6NI7XxL7EOJuUSfxPdazau8U/TNgmZkNAZ4Hrqq+SjIQ1fMQxVKqGZaXtD6EvpjHPmnvsZ5VNWzR3TeF+fH3uvs8d38trorJwJLWS++09Q2oRZ1smstFEiHOx9slSRqfNaoWdXIpoEtipDFQqG9A+pPmchHpQxqWJ/1JAV16qKf5x+tBWvsGJJmUcpGj6m3+8XqQ1r4BSSYFdDkqjR14SZDGvgFJJqVc5Ch14InUNwV0OUodeCL1TQFdjlIHnkh9Uw5djlIHnkh9U0CXHtSBJ1K/lHIREUmJfm2h79u3j6VLl/ZYNmXKFM466yy6u7tZtmxZ3jbt7e20t7fT1dXF8uXL88o7OjqYOnUq+/fv584778wrP+ecc5g8eTJ79+7lnnvuySufMWMGEydOZPfu3axatSqvfObMmYwdO5bt27ezZs2avPI5c+YwatQonn/+edatW5dXPnfuXFpaWti6dSvr1+c/P/vSSy9lxIgRbN68mQ0bNuSVL1iwgKamJjZt2sSmTZvyyhcuXEhDQwOPPvooW7ZsySu/8sorAXj44Yd55plnepQ1NDSwcOFCANauXcu2bdt6lDc1NbFgwQIAHnjgAXbs2NGjfPjw4cyfPx+AVatWsXv37h7lI0eO5KKLLgLg7rvvZt++fT3KR40axZw5cwBYsWIFBw4c6FE+ZswYZs2aBcDy5cvp6urqUT5hwgTOO+88AJYtW0Z3d3eP8kmTJjF9+nSAvM8d6LOnz179fvaiqIUuIpIS5u79drCOjg4v1BIQEZFoZvaYu3eUWk8tdBGRlFBAFxFJCQV0EZGUqHqUi5kNAjYAO919bvVV6qma5zEm+VgiInGLY9ji54CngOEx7KuH/pzOVVPHiki9qyrlYmZjgI8D34mnOj0Vm861no8lItIXqs2hfwP4InAkagUzW2RmG8xsw549e3q18/6czlVTx4pIvas4oJvZXOAVd3+s2Hrufou7d7h7R2tra6+O0Z/TuWrq2PTQY/RkoKqmhX4ucLGZ/Rr4IXChmf0gllqF+nM6V00dmw6ZvpCdnQdxjvWFKKjLQFBxQHf3a919jLuPBy4DHnT334utZgSdkV+dfwZtzY0Y0NbcyFfnn9EnnZT9eSzpO+oLkYEs8dPn9ud0rpo6tv6pL0QGslhuLHL3n/TFGHSR3lJfiAxkulNUUkV9ITKQJT7lItIbeoyeDGQK6JI66guRgUopFxGRlFBAFxFJCQV0EZGUUEAXEUkJBXQRkZRQQBcRSQkFdBGRlFBAFxFJCd1YVCU9h1REkkIBvQp6DqmIJIlSLlXQ3NsikiQK6FXQ3NsikiQK6FXQ3NsikiQK6FXQ3NsikiQVB3QzG2tmD5nZU2a2xcw+F2fF6oGeQyoiSVLNKJdDwNXu/riZnQg8Zmb3u/svY6pbXdDc2yKSFBW30N39JXd/PPz9N8BTgCKbiEiNxJJDN7PxwDTgZwXKFpnZBjPbsGfPnjgOJyIiBVQd0M3sBOD/AZ939wO55e5+i7t3uHtHa2trtYcTEZEIVQV0M2sgCObL3H1FPFUSEZFKVDPKxYDvAk+5+z/EVyUREamEuXtlG5p9CPgv4EngSLj4L939P4tsswd4oaID9p0WYG+tK1Em1TV+9VJPUF37Sj3U9V3uXjJnXXFATwsz2+DuHbWuRzlU1/jVSz1Bde0r9VTXUnSnqIhISiigi4ikhAI63FLrCvSC6hq/eqknqK59pZ7qWtSAz6GLiKSFWugiIimhgC4ikhIDKqCbWbOZ3WFmT4fT/p5jZieb2f1m9mz470kJqOdkM9uU9XPAzD6fxLoCmNkXwimUN5vZbWZ2vJlNMLOfhXX9kZkNqXU9Aczsc2E9t5jZ58NliTivZnarmb1iZpuzlhWsmwX+0cx+ZWa/MLP3JaCuvxue1yNm1pGz/rVhXbea2ewa13NJGAN+YWZ3mllzresZlwEV0IFvAqvc/TTgTIIZIq8B1rj7e4A14euacvet7t7u7u3A+4Eu4E4SWFczawP+HOhw96nAIOAy4H8DXw/r+hrwmdrVMmBmU4HPAh8g+PvPNbP3kJzzuhSYk7Msqm4fBd4T/iwCvt1PdcxYSn5dNwPzgXXZC83sdILPxJRwm/9rZoPoH0vJr+f9wFR3fy/wDHBtAuoZiwET0M1sODCDYLoC3P1td+8ELgG+F672PWBebWoYaSbwnLu/QHLrOhhoNLPBQBPwEnAhcEdYnpS6/hbwiLt3ufshYC1wKQk5r+6+Dng1Z3FU3S4Bvu+BR4BmM3tn/9S0cF3d/Sl3L/SE9EuAH7r7W+6+DfgVwZdqn4uo533h3x/gEWBMresZlwET0IGJwB7gX81so5l9x8yGAe9w95cgmOMdOKWWlSzgMuC28PfE1dXddwJfA14kCOT7gceAzqz/NDtIxlz5m4EZZjbSzJqAjwFjSeB5zRJVtzZge9Z6STnHhSS5rp8Gfhz+nuR6lmUgBfTBwPuAb7v7NOANEpCyKCbMO18M3F7rukQJc7qXABOA0cAwgnRArpqPj3X3pwhSQfcDq4AnCJ68VY+swLKan+MIiayrmV1H8PdflllUYLWa17M3BlJA3wHscPfMQzjuIAjwL2cuVcN/X6lR/Qr5KPC4u78cvk5iXWcB29x9j7t3AyuA6QQpgMwjDscAu2pVwWzu/l13f5+7zyC4FH+WZJ7XjKi67SC4ushIzDkuIHF1NbMrgLnAQj92M07i6tlbAyagu/tuYLuZTQ4XzQR+CdwFXBEuuwL4jxpUL8rlHEu3QDLr+iJwtpk1hVMqZ87rQ8DvhOskpa6Y2Snhv+MIOvBuI5nnNSOqbncBnwpHu5wN7M+kZhLoLuAyMxtqZhMIOnJ/XqvKmNkc4EvAxe7elVWUqHpWxN0HzA/QDmwAfgGsBE4CRhKMHng2/PfkWtczrGsTsA8YkbUsqXX9CvA0QY7634ChBH0WPyfoWLodGFrreoZ1/S+CL5wngJlJOq8EXy4vAd0ErcXPRNWNID3wf4DnCKaw7khAXS8Nf38LeBlYnbX+dWFdtwIfrXE9f0WQK98U/vxzresZ149u/RcRSYkBk3IREUk7BXQRkZRQQBcRSQkFdBGRlFBAFxFJCQV0EZGUUEAXEUmJ/w9tWamI/KmFZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c1daa0240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "dataA=pd.DataFrame(columns=['subject','values','id'])\n",
    "dataB=pd.DataFrame(columns=['subject','values','id'])\n",
    "\n",
    "A=np.array([    \n",
    "56.9,63.2,65.5,73.6,74.1,77.1,77.3,77.5,77.8,78.9,79.5,80.8,81.2,81.9,82.2,83.1,84.4,84.9,86,86.3,86.3,86.6,86.6,86.6,\n",
    "87.1,87.5,87.8,88.6,89.3,89.6,90.3,91.1,92.1,93.5,94.5,94.6,95,95.2,95.3,95.6,95.9,96.4,97.2,97.5,97.9,98.2,98.5,98.8,\n",
    "98.9,99,99.3,99.3,99.9,100.1,101,101,101.5,101.5,101.5,101.8,101.8,102.8,102.9,103.2,103.8,104.4,104.8,105.1,105.5,\n",
    "105.7,106.1,106.8,107.2,107.4,107.5,107.5,108,108.2,108.6,109.1,110.1,111.2,111.7,111.7,112,113.1,116,116.7,118.8,\n",
    "119.7,120.7,122.8,124.7,126.4,127.6,128.2,129.6,130.4,133.2])\n",
    "\n",
    "B=np.array([52.9,59.2,63,66.2,64.8,69,67.1,70.1,69.2,73.8,71.8,73.3,73.1,74.7,74.1,74.1,76,75.4,74.6,79.2,77.8,80.8,77.6,77.5,78.6,\n",
    "78.7,81.5,79.3,78.9,85.9,80.7,80.6,82.8,86,84.3,87.6,84,85.9,84.4,85.2,85.2,89.2,87.8,88,88.7,91.2,91.8,92.5,88,93.5,\n",
    "89,89.4,89.2,91.3,90.4,91.2,91.4,93,91.2,92,91.8,96.8,92.8,94,93.5,95.8,97.1,97.3,95.1,95.8,95.5,95.9,95.4,97.3,97.7,\n",
    "93,97.6,96.1,96.2,99.5,99.8,105.3,103.6,100.2,100,98.8,110,103.5,109.4,112.1,111.3,108.6,112.4,113.8,115.6,118.1,116.8,\n",
    "121.6,115.8])\n",
    "\n",
    "dataA['subject']=np.arange(0,99,1)\n",
    "dataA['values']=A\n",
    "dataA['id']=np.ones(99)\n",
    "\n",
    "\n",
    "dataB['subject']=np.arange(0,99,1)\n",
    "dataB['values']=B\n",
    "dataB['id']=np.ones(99)*2\n",
    "\n",
    "data=pd.concat([dataA,dataB])\n",
    "data=data.reset_index(drop=True)\n",
    "\n",
    "md=smf.mixedlm('values~subject', data, groups=data['id'], re_formula='~subject')\n",
    "mdf=md.fit()\n",
    "\n",
    "print(mdf.summary())\n",
    "\n",
    "plt.title('Scatterplot for each method')\n",
    "plt.scatter(A,B)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "bland_altman_plot(A, B)\n",
    "plt.title('Bland-Altman Plot of two methods')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem:\n",
    "Potatoes measured by two sclaes, A and B.<br>\n",
    "(a) Compute the Grubbs estimators for these data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_b: 1211.3333333333333\n",
      "beta_0: -17.0\n",
      "sigma_b: 586217.0238095239\n",
      "sigma_e_1: -25305.468253968284\n",
      "sigma_e_2: -51974.13492063503\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib as plt\n",
    "\n",
    "A=np.array([135, 940, 1075, 925, 2330, 2870, 1490, 2110, 650, 1380, 970, 1000, 1640, 345, 310])\n",
    "B=np.array([165, 910, 1060, 925, 2290, 2850, 1425, 2050, 630, 1370, 1000, 1000, 1575, 345, 320])\n",
    "\n",
    "A_mean = A.mean()\n",
    "B_mean = B.mean()\n",
    "S_A = np.std(A)\n",
    "S_B = np.std(B)\n",
    "S_AB = np.sum((A-A_mean)*(B-B_mean))/14\n",
    "\n",
    "mu_b = A_mean\n",
    "beta_0 = B_mean-A_mean\n",
    "sigma_b = S_AB\n",
    "sigma_e_1 = S_A**2-S_AB\n",
    "sigma_e_2 = S_B**2-S_AB\n",
    "\n",
    "print(\"mu_b:\", mu_b)\n",
    "print(\"beta_0:\", beta_0)\n",
    "print(\"sigma_b:\", sigma_b)\n",
    "print(\"sigma_e_1:\", sigma_e_1)\n",
    "print(\"sigma_e_2:\", sigma_e_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the estimators from Grubbs are collected:\n",
    "$ \\hat{\\mu}_b=1211.3$, $\\hat{\\beta}_0=-17$, ${\\hat{\\sigma}}_b^2=586217$, ${\\hat{\\sigma}}_{e1}^2=-25305$, and ${\\hat{\\sigma}}_{e2}^2=-51974$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) Fit mixed-effects model using ML method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    packet  scale   id\n",
      "0        1    135  1.0\n",
      "1        2    940  1.0\n",
      "2        3   1075  1.0\n",
      "3        4    925  1.0\n",
      "4        5   2330  1.0\n",
      "5        6   2870  1.0\n",
      "6        7   1490  1.0\n",
      "7        8   2110  1.0\n",
      "8        9    650  1.0\n",
      "9       10   1380  1.0\n",
      "10      11    970  1.0\n",
      "11      12   1000  1.0\n",
      "12      13   1640  1.0\n",
      "13      14    345  1.0\n",
      "14      15    310  1.0\n",
      "15       1    165  2.0\n",
      "16       2    910  2.0\n",
      "17       3   1060  2.0\n",
      "18       4    925  2.0\n",
      "19       5   2290  2.0\n",
      "20       6   2850  2.0\n",
      "21       7   1425  2.0\n",
      "22       8   2050  2.0\n",
      "23       9    630  2.0\n",
      "24      10   1370  2.0\n",
      "25      11   1000  2.0\n",
      "26      12   1000  2.0\n",
      "27      13   1575  2.0\n",
      "28      14    345  2.0\n",
      "29      15    320  2.0\n",
      "                   Mixed Linear Model Regression Results\n",
      "============================================================================\n",
      "Model:                   MixedLM       Dependent Variable:       scale      \n",
      "No. Observations:        30            Method:                   REML       \n",
      "No. Groups:              2             Scale:                    573880.0614\n",
      "Min. group size:         15            Likelihood:               -234.5395  \n",
      "Max. group size:         15            Converged:                No         \n",
      "Mean group size:         15.0                                               \n",
      "----------------------------------------------------------------------------\n",
      "                           Coef.    Std.Err.   z    P>|z|   [0.025   0.975] \n",
      "----------------------------------------------------------------------------\n",
      "Intercept                  1405.333  609.636  2.305 0.021   210.469 2600.197\n",
      "packet                      -25.313  536.624 -0.047 0.962 -1077.075 1026.450\n",
      "Intercept RE             573880.061                                         \n",
      "Intercept RE x packet RE      0.000                                         \n",
      "packet RE                573880.061                                         \n",
      "============================================================================\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/regression/mixed_linear_model.py:2001: ConvergenceWarning: Gradient optimization failed.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/regression/mixed_linear_model.py:2039: ConvergenceWarning: The Hessian matrix at the estimated parameter values is not positive definite.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py:1029: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return np.sqrt(np.diag(self.cov_params()))\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "dataA=pd.DataFrame(columns=['packet','scale','id'])\n",
    "dataB=pd.DataFrame(columns=['packet','scale','id'])\n",
    "\n",
    "dataA['packet']=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]\n",
    "dataA['scale']=A\n",
    "dataA['id']=np.ones(15)\n",
    "\n",
    "\n",
    "dataB['packet']=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]\n",
    "dataB['scale']=B\n",
    "dataB['id']=np.ones(15)*2\n",
    "\n",
    "data=pd.concat([dataA,dataB])\n",
    "data=data.reset_index(drop=True)\n",
    "print(data)\n",
    "\n",
    "md=smf.mixedlm('scale ~ packet', data, groups=data['id'], re_formula='~ packet')\n",
    "mdf=md.fit()\n",
    "\n",
    "print(mdf.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the error variance estimates are ${\\hat{\\sigma}}_{e1}^2=609.6^2=379881$, and ${\\hat{\\sigma}}_{e2}^2=536.6^2=287939.6$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) The error variance estimation is ${\\hat{\\sigma}}_{e1}^2=-25305$, and ${\\hat{\\sigma}}_{e2}^2=-51974$ for Grubbs estimation, and ${\\hat{\\sigma}}_{e1}^2=609.6^2=379881$, and ${\\hat{\\sigma}}_{e2}^2=536.6^2=287939.6$ for mixed-effects model using ML. The mixed-effects model can be fit by the ML method, and this estimates are used to perform inference on measures of agreements. The Grubbs estimates, lead to a negative estimates of error variance."
   ]
  },
  {
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}