How to fit a line to refraction data.

How_to_fit_a_line_GEOPH_326.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lab 2 GEOPH 326 - February 14, 2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the refraction seismology method one needs to fit lines to $t-x$ picks. I will explain linear regression via the method of least squares. We start by defining the  model that provides the relationship bewtween variables of interest and observations:\n",
    "\n",
    "\n",
    "$$ t_i = a \\, x_i + b + \\epsilon_i \\quad i=1\\dots N$$\n",
    "\n",
    "$t_i, x_i$ are picks extracted from refraction profiles. The parameters $a$ and $b$ are unknows to be found from the pairs $t_i,x_i$. The parameters $a$ and $b$ will be used in your lab to estimate velocity, depth and dip of a geological interface. \n",
    "\n",
    "\n",
    "<b>CASE I</b>: We first assume that the line is such that $t(x=0)=0$. In other words we assume $b=0$. This is\n",
    "the line that you will use to fit $v_1$ to the direct wave. In this case we have  \n",
    "\n",
    "$$t_i = a\\, x_i + \\epsilon_i\\,.$$\n",
    "\n",
    "We now form an error energy function\n",
    "\n",
    "$$J(a) = \\sum_i \\,e_i^2 = \\sum_i \\, (t_i - a \\, x_i)^2$$\n",
    "\n",
    "and we compute $a$ such that the $J(a)$ is minimized. The condition for the minimum of $J(a)$ is given by \n",
    "\n",
    "$$ \\frac{\\partial J(a)}{\\partial a}= 2 \\sum_i \\,(t_i - a \\, x_i) \\,x_i =0\\,.$$\n",
    "\n",
    "The latter leads to\n",
    "\n",
    "$$ \\sum_i \\, t_i x_i - a \\sum_i \\,x_i^2$$ \n",
    "\n",
    "therefore, \n",
    "\n",
    "$$ {\\hat{a}} \\,\\,= \\frac{\\sum_i \\, t_i x_i }{\\sum_i \\, x_i x_i}\\,.$$\n",
    "\n",
    "The result $\\hat {a}$ is the estimated slope.\n",
    "\n",
    "<b> CASE II</b> Fit  both $a$  and  $b$. This is the procedure for fitting a line to the head wave.\n",
    "In this case we have the following model \n",
    "\n",
    "$$t_i = a\\, x_i + b \\epsilon_i\\,.$$\n",
    "\n",
    "We now form an error energy function than depends on $a$ and $b$\n",
    "\n",
    "$$J(a,b) = \\sum_i \\,e_i^2 = \\sum_i \\, (t_i - a \\, x_i - b)^2$$\n",
    "\n",
    "and we compute $a$ and $b$ such that $J(a,b)$ is minimized. This is accomplished by finding the minimum of $J(a,b)$ with respect to $a$ and $b$ \n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial J(a,b)}{\\partial a}&=& 2 \\sum_i \\,(t_i - a \\, x_i-b) \\,x_i =0\\\\ \n",
    "\\frac{\\partial J(a,b)}{\\partial b}&=& 2 \\sum_i \\,(t_i - a \\, x_i-b)  =0 \n",
    "\\end{align}\n",
    "\n",
    "The last two equations can be written down via a system of two equations with two unknows \n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\sum_i t_i x_i \\\\\n",
    "\\sum_i t_i \n",
    "\\end{pmatrix} \n",
    "=\n",
    "\\begin{pmatrix}\n",
    "\\sum_i x_i^2 &  \\sum_i x_i\\\\\n",
    "\\sum_i x_i &  N\n",
    "\\end{pmatrix} \n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b \n",
    "\\end{pmatrix}\\,.\n",
    "$$\n",
    "\n",
    "We now write the last system in matix-times-vector form as follows\n",
    "\n",
    "$${\\bf v} = {\\bf M} {\\bf u}$$ from where we can estimate the solution \n",
    "\n",
    "$${\\hat {\\bf u}} = {\\bf M}^{-1} {\\bf v}$$\n",
    "\n",
    "Once we have computed the vector $\\hat{\\bf u}$, we can estimate $\\hat{a}$ and $\\hat{b}$\n",
    "from the elements of $ \\hat{\\bf u}$\n",
    "\n",
    "\\begin{align}\n",
    "{\\hat a} = {\\hat u}(1)\\\\\n",
    "{\\hat b} = {\\hat u}(2)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Next, I will write a function to compute $a$ or $a$ and $b$. Then, we will test the function with synthetic data.\n",
    "You need to understand and use the function <code>fit_a_line</code> and use it to fit direct and head waves from the data provided in the lab. Then, estimated slopes and intercepts will be used to estimate velocities, critial angle, depth under the source and dip of the interface. Easy! \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fit_a_line (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fit_a_line(xi,ti;flag=true)\n",
    "#\n",
    "# Function to estimate the parameters of a regression line  \n",
    "# using the least-squares method\n",
    "#\n",
    "# t[i]  = a x[i] + b + e[i]\n",
    "# \n",
    "# flag==true   --- >  Estimate only a (assume t=0 for x=0)\n",
    "# flag==flase  --- >  Estimate a and b\n",
    "#\n",
    "# ----------------------------------\n",
    "# M D Sacchi\n",
    "# GEOPH 326, University of Alberta\n",
    "# 2018\n",
    "# ----------------------------------\n",
    "    \n",
    "    if flag==true\n",
    "        a_estimated = sum(ti.*xi)/sum(xi.*xi);\n",
    "    else\n",
    "        M = zeros(2,2); v = zeros(2,1)\n",
    "        \n",
    "        M[1,1]=sum(xi.^2)\n",
    "        M[1,2] = sum(xi)\n",
    "        M[2,1] = M[1,2]\n",
    "        M[2,2] = length(xi)*1.0\n",
    "        v[1] = sum(ti.*xi)\n",
    "        v[2] = sum(ti)\n",
    "        u = M\\v\n",
    "        a_estimated=u[1]\n",
    "        b_estimated=u[2]\n",
    "    end\n",
    "    if flag==true;\n",
    "        return a_estimated\n",
    "    else\n",
    "        return a_estimated,b_estimated\n",
    "    end\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAFKCAYAAAB/xnmMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl8TNf/x/HXnQhZxRpbQ4JISrUErTWW1FL9oraiTcUuShBUWyrK92unipJSVXxpUK3aWlsQsQWJ7dvaFaUatSQSa5bz+2N+mRqTqJDkTpLP8/GYBzn3zNz3HTL55Nx7ztWUUgohhBBCCJEnGPQOIIQQQgghso4Ud0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0IIIYQQeYgUd0JkwoULFzAYDPTs2VPvKFnmWY4pL74PQgiRV0hxlw8ZDAazR4ECBShevDhNmzZlyZIlesfLFTRN0ztClnv8mJ6mgMuK92H37t3UqVOHkiVL0qFDh+d+vad16tQpOnbsyLBhwxg+fDgBAQFcu3Ytx/YvhBDZpYDeAYQ+NE1j7NixACQlJXHmzBnWrFlDREQEhw4dYs6cOTontE4vvPACJ0+exMXFRe8oWSajY0or3LK7kG3YsCF79uzB0dGRunXrZuu+0ty+fZvXX3+dKVOm8M477wAwadIkXn/9daKjo7G1tc2RHEIIkR2kuMvHQkJCzL7eu3cvvr6+zJs3j+HDh+Pu7q5PMCtWoEABqlSponeMLJXRMSmlzP7MTjExMaSkpNCkSZNs3xfA1KlTSUpKomvXrqa2wMBAxo4dy6JFi+jfv3+O5BBCiOwgp2WFSf369fHy8kIpRUxMjMX2qKgoOnXqROnSpSlUqBDly5cnMDCQq1evpvt6SilmzZpF1apVsbe354UXXiAoKIj4+Hjc3d3x8PAw9X30FODp06fp0qULrq6u2NjYEBER8cw51q1bh5+fH2XKlMHOzo5y5crRpEkTQkNDn6nfk05Vrlq1Cl9fX1xcXHBwcODll19m8uTJPHz40KLvo69z4cIFunbtSokSJbC3t6dOnTps3Lgx3ff0cYmJiRQsWJCGDRuatd+7dw87OzsMBgPLli0z2xYaGorBYGDx4sUZHtOnn35KxYoVAViyZInZafylS5daHMuz5k8TERGBs7MztWvXztTzntWqVat47bXXMBj+/ggsWrQo3t7efPfddzmSQQghsouM3Il0PX5aatGiRfTr1w97e3vatm2Lm5sbp0+fZuHChaxfv579+/fj5uZm9pyBAwfy5ZdfUq5cOfr374+trS3r1q3jwIEDJCcnU7BgQYv9njt3jrp16+Ll5cV7773HvXv3zE4XZibHggULCAwMpEyZMrRr144SJUpw7do1jh49yuLFixkwYECm+j3q8VOVo0aNYvLkyZQsWRJ/f3+cnJz46aefGDVqFJs3b2bLli3pnuq7ePEir732GpUqVSIgIIAbN26wcuVK2rVrx7Zt2/5xJMvJyYnXXnuNAwcOkJiYiJOTEwB79uwxFZXh4eH4+/ubnhMeHo6mafj5+WV4TE2bNiU+Pp5Zs2ZRo0YN3nrrLdO2GjVqmP5+4cKF58qfJiIiggYNGpgVW9klISGBs2fP0rJlS4ttZcqU4cCBA9meQQghspUS+Y6macpgMFi0R0REKIPBoOzs7NSff/5paj916pSytbVVnp6e6o8//jB7Tnh4uLKxsVHt27c3a9+1a5fSNE15e3ur+Ph4U/vDhw+Vr6+v0jRNeXh4mNp/++03pWma0jRNjR49Ot3cmc3h4+Oj7Ozs1F9//WXxWjdu3Mh0v0dz9uzZ09S2d+9epWmaqlChgoqNjTW1JycnqzZt2ihN09TEiRPTfR1N09T48ePNtm3evFlpmqZat26d7vvwuJCQEKVpmtq4caOp7aOPPlIFChRQfn5+ys3NzdSekpKiihUrpipXrvzEY1JKqQsXLqTbntX5k5OTVeHChVVgYKAaMWKEGj58uGrZsqU6fPjwUz0/s44fP640TVOjRo2y2Na+fXtlMBhUUlJStuxbCCFygozc5VNKKcaNG4dSiqSkJM6ePcuaNWvQNI3p06dTqlQpU9/Q0FCSk5OZNWsWZcqUMXudZs2a0aZNG9avX8+dO3dwdHQEMM26HT16NIULFzb1t7W1ZdKkSRanEdOULl3aNNHjcc+Sw8bGhgIFLP+bFytWzOzrp+2XnkWLFgHwySef4OrqavaaM2bM4KeffmLhwoV8/PHHFs91d3fnk08+MWtr0aIFbm5uHDx48B/3DeDn58e///1vwsPDad26NWAcnatduzYdOnRg0KBBnDlzBk9PT44cOcKtW7fo3LnzP76ueopr7bIif0xMDAkJCfz++++sXbsWGxsbFi5cSPv27Tl79iw2NjZm/QcMGMCRI0ee6rXTTJ06lUaNGgHGyRRAuiPHjo6OKKWIi4ujRIkSmdqHEEJYCynu8rFx48aZfa1pGosWLSIgIMCsfd++fQDs3LmTqKgoi9e5du0aKSkpnDp1Ch8fHwAOHz6MpmnpFnGvvfaaxQ/sNK+88kqGMxWfNsfp06epWbMm/v7+DB8+nKpVq9K1a1d8fX1p0KABJUuWNHve0/bLSExMDJqm0axZM4ttnp6elCtXjgsXLpCQkICzs7PZ9ho1aqQ7G9XNzS3dY0xPvXr1sLe3Jzw8HID4+HgOHz7Mhx9+aMq0fft2PD092b59O0C6WZ9FVuSPiIjAxcWF7777zvT/ws3NjYsXLxITE0OdOnXM+j9+HWRmpe0jvdxJSUkAJCcnP9c+hBBCT1Lc5VOappGSkgIYL77fu3cvvXv3pn///pQvX56mTZua+t64cQOAadOmPfH17ty5Y/o6Pj4ewGwEMI2NjQ3FixdP93VKly6d4T6eNkdiYiIAwcHBlChRgnnz5jF79mw+//xzNE2jcePGTJs2jVq1amWqX0bSjvXx0cQ0ZcqU4fLly8TFxVkUd0WKFEn3OQUKFCA1NfWJ+01ja2tLgwYNCA8P5/r16+zZs4eUlBT8/Pzw9vamTJkyhIeH079/f8LDwzEYDFlW3GVF/p07d9KoUSPs7e1Nbb///juA2f+prPKkov3OnTtomma6dlEIIXIjmS0rsLe3x8/Pj/Xr15OSkkJAQAD37t0zbXdxcUHTNG7fvk1qamq6j5SUFNNpL8B0KvbPP/+02F9KSoqpUHvck9ZUe5Yc7733Hvv27ePGjRts3LiR3r17s2vXLlq2bMn169cz3S+jXECGs4bT2rNzbTw/Pz+UUoSHhxMeHo6dnR0NGjQAjKN027dv5+HDh0RGRlKtWjWrOuW4b98+fH19zdr2798PgLe3d5bvr1SpUmiaxq1btyy23blzhyJFikhxJ4TI1aS4EybVq1enb9++XL58mZkzZ5ra69Wrh1KKXbt2PfVr+fj4oJRi9+7dFtv2799vGjXMjGfJkcbFxYU33niDBQsW0KNHD27evElkZOQz93tU2rHu3LnTYtvZs2e5fPkyHh4eZtceZrW0ma/h4eHs2LGDBg0amK4p8/Pz4+bNm8ybN4+7d+9azJLNSNrpy2f5t3pa169f59atW7z66qumttTUVDZt2kSjRo3SHcnt168fderUydTj0f+Hjo6O1KxZk0uXLlm89tmzZ81mAwshRK6k21QOoZuMZssqpdSVK1eUnZ2dKlq0qLp165ZSSqmTJ0+qggULqipVqqjTp09bPOfBgwdq165dZm0RERFK0zTl5eVlNlv2wYMHT5wtm97MzDSZzbF9+/Z0X+df//qX0jRNbdq0KVP9MsqZNlvWw8PDbMZtcnKyateu3RNny2Z0vI0bN87w3yg9ycnJqkiRIsrV1VVpmqYmTZpk2nbx4kWlaZoqVaqU0jRNrV+//qmyJCQkKE3TVOPGjS32l1X579y5owwGgzpz5oyp7ccff1QGg0FFRkb+4/OfVUhIiCpXrpxZ29mzZ5WmaWrevHnZtl8hhMgJcs2dMFO2bFkCAwOZNWsWU6dOZeLEiXh5ebFo0SJ69epFtWrVaNWqFZ6eniQlJXHp0iUiIyMpVaoUv/76q+l1fH196devHwsWLKBatWp06NABW1tb1q9fT9GiRSlbtmym1zTLbI727dvj7OxM3bp1qVChAkopIiMjOXToELVr1+b111/PVL+M1KtXj5EjRzJ16lReeuklOnXqhIODAz///DO//PILjRo14oMPPsjkv0Tm7gxhY2NDkyZNWLt2LYDZ6Fz58uWpVKkS586do0CBAjRu3PipXtPJyYm6desSGRmJv78/np6e2NjY0K5dO4trB581v4ODA6+//jq//vorlStX5urVqwwYMICJEydmOKM6KwwYMIDZs2ezbNky0xqAc+bMoVq1avTt2zfb9iuEEDlC39rSUkBAgGn9rPQej69vJjLvSSN3SikVGxurHB0dlZOTk7p27Zqp/fjx46pHjx6qQoUKqlChQqp48eKqevXqKjAwUO3YscPidVJTU9XMmTOVt7e3KlSokCpXrpwaNGiQio+PV05OTqpmzZqmvk8zcpfZHF9++aVq3769qlixonJwcFDFihVTPj4+atq0aSoxMTHT/f4p54oVK1TDhg2Vs7OzsrOzUy+99JKaOHGievDggUXffzreJk2aZGrkTiml5syZozRNU0WKFFGpqalm2/r37680TVN169bNVJazZ8+qNm3aqOLFiyuDwaAMBoNasmTJE9fAy2z+S5cuqS5duqjg4GDVpk0b9cMPPzzV857XkSNH1JtvvqmCg4NVnz59VMeOHdXly5dzZN9CCJGdNKVy4MaRmbB//37Onz9v1paamkpgYCAeHh4cP35cp2Qiq5w5cwYvLy+6devG8uXL9Y4jhBBC5ClWd1q2bt261K1b16xt9+7d3L17l3fffVenVOJZxMbGUrJkSbPTr3fv3mXo0KGA8XSoEEIIIbKW1RV36fn222/RNI133nlH7ygiE2bOnElYWBhNmzaldOnS/Pnnn4SHh3PlyhVat25Np06d9I4ohBBC5DlWd1r2cUlJSZQpU4aqVas+0xIYQj/bt29n+vTpHDlyhJs3b2Jra0uVKlV45513GDp0aIZ3qRBCCCHEs7P6kbvNmzdz8+ZNOSWbCzVr1izL7oQghBBCiKdj9cXdt99+S8GCBXn77bfT3X79+nU2b96Mu7u72e2LhBBCZK979+5x4cIFWrZs+Vx3Pbl79y4nT57MwmRC5G3e3t44ODhk3EHfybpPlpCQoBwcHFTbtm0z7LNs2TIFyEMe8pCHPHR6LFu27Lk+66Ojo3U/BnnIIzc9oqOjn/g9ZdUjdz/++CP37t174ilZd3d3AJYtW8aLL76YQ8n0ExwcbHZrsLwsvxxrfjlOkGPNa06cOIG/v7/pc/h55ZfPcSGeVdr33D+x6uJu+fLlODs707Zt2wz7pJ2KffHFF/Hx8cmpaLpxcXHJF8cJ+edY88txghxrXpVVl8Tkl89xIbJb5u7/lIP++usvtm3bRvv27bGzs9M7jhBCCCFErmC1xd3KlStJSUmRWbJCCCGEEJlgtcXdt99+S6lSpf7xpu1CCCGEEOJvVlvc7d27l6tXr6Jpmt5RrEq3bt30jpBj8sux5pfjBDnWrKCUYujQUSjrXn9eCKEjqy3uRPrkh2Pek1+OE+RYs0J0dDRz584hJiYmW15fCJH7SXEnhBC5SGjodyQnzyA09Du9owghrJQUd0IIYeVCQibh6uqFp+cbbNx4GujLhg2nqFy5Fa6uXoSETNI7otDR4sWLMRgMLFmyRO8ouUJ+eL+kuBNCCCs3ZswIxowZSUJCAWJj1wAasbFrSEy0ZcyYkYwZM0LviCILHTp0iJ49e1KxYkUcHBxwcXHh5ZdfZuTIkfzxxx8ZPk+uUTfauXMnBoOBcePGpbtd0zTTI6+S4k4IIaycra0tQUG9cXEx/2Hk4qIRFNQbW1tbnZKJrPbhhx/y6quv8u2331K1alWGDBlCnz59cHBwYPr06VSpUoXvv/9e75i5QkbFW/v27Tlx4gRvvfVWDifKOVLcCSFELpGamoy9/ULKlGmOvf1CUlOT9Y6kq5yeOZzd+xs/fjzTpk3Dw8ODI0eOsGHDBiZNmsSMGTPYv38/q1evJjU1la5du7Jz585syZCXZPTvVLhwYapUqULhwoVzOFHOkeJOCCFyCR8fD6ZM0Th7di1Tpmj4+HjoHUlXOT1zODv3d+HCBf79739TsGBB1q1bl+49djt06MDMmTNJSUlhwIABFsWLUoqNGzdSv359nJycKFasGJ07d+bs2bMWrxUbG8uIESPw8vLCycmJokWL4u3tTc+ePfntt98s+m/evJnWrVtTokQJ7OzsqFy5MiNHjiQ+Pt6ir7u7Ox4eHiQkJDBs2DDc3d0pWLAg48aNIzAwEIPBwLp169J9H6KiojAYDLz99tumttOnT/PRRx9Ru3ZtSpYsiZ2dHe7u7vTv358rV66YPb9Hjx40a9YMgHHjxmEwGEyPXbt2AU++5i46OpqOHTvi6upq2s/AgQP5888/Lfr26NEDg8HAxYsXmT9/PtWrV8fe3p7SpUvTv39/bt++bfGcY8eO0a1bN9zd3bGzs8PV1ZVatWoRHBxMcnLW/bJm1feWFUII8beVK+ea/h4U1JugoN46ptHfozOHFy6slav3980335CSkkKnTp2oVq1ahv369OnDuHHjOHXqFBERETRp0sS07YcffuDnn3+mQ4cONGvWjMOHD/P999+zY8cO9u7dS5UqVQC4e/cuDRo04Pz587Ro0YJ27dqhlOLChQusW7eOzp074+Hx9y8O48aNY9y4cRQvXpw2bdrg6urK0aNHmT59Oj/99BP79u3D2dnZ1F/TNB4+fEjTpk2Ji4ujVatWFC5cmIoVK9KyZUsWLFjA0qVL071vfFrB1aNHD7Pjmj9/Ps2aNaNhw4YULFiQ//3vfyxcuJD169dz6NAhypYtCxhPuWqaxpIlS2jSpInZ++Pu7m62r8dP227YsIGOHTuiaRqdOnWiQoUKHDp0iNDQUNauXcvu3bstXgPggw8+YMuWLbRt25ZWrVqxfft2vvrqK86ePUt4eLip37Fjx3jttdewsbGhbdu2eHh4cPv2bc6cOUNoaCgTJkygQIEsKstULhcdHa0AFR0drXcUIYTIV7Lq8zczrzNmzERVsmQVVblyK1Wq1FsKUlWpUm+pSpVaqpIlq6gxYyY+Vxa99tesWTOlaZpauHDhP/Z99913laZpasKECUoppb755hulaZrSNE1t3LjRrO+sWbOUpmnKz8/P1LZu3TqlaZoaNmyYxWsnJSWphIQE09fbt29XmqapBg0aqPj4eLO+ixcvVpqmqeDgYLP2ChUqKE3TVPPmzdXdu3ct9uHl5aUKFSqkbt68adZ+//59VbRoUVW6dGmVkpJiar9y5Yp6+PChxets2bJF2djYqAEDBpi179ixQ2mapsaNG2fxHKX+fr+WLFliaktISFDFihVTBQoUULt37zbrP2XKFKVpmmrRooVZe0BAgNI0TVWoUEH9/vvvpvbk5GTl6+urNE1TBw4cMLUPGzZMaZqm1q1bZ5EpLi5Opaamppv3UU/7vSKnZYUQIjdZtQr8/SGf3qEip2cO59T+rl69CoCbm9s/9n3hhRcALGbO+vn50bp1a7O2QYMGUbFiRbZv386lS5fMttnZ2Vm8doECBXBycjJ9PXv2bAC++uori2vUAgICeOWVV1i+fLnF62iaxowZM7C3t7fYFhAQwMOHDwkLCzNrX79+PXFxcbz77rsYDH+XJ2XLlk130lDz5s2pWrUqmzdvttiWWWvXruXWrVt06dKFBg0amG0bPnw4FSpUYOvWrfz+++8Wzw0JCTH9mwDY2NjQs2dPAA4ePGjRP7333cXFJUtn70pxJ4QQuUFcHLz3HnTpAg8fwv37eifSRU7PHM5NM5UbN25s0WYwGGjYsCEAR44cAaBJkyaUK1eOyZMn88YbbzB79mxiYmJITU21eP6+ffuwtbVl1apVfPrppxaPhw8f8tdff3Hr1i2z59nZ2VG9evV0c3bv3j3da97SOyWbZtmyZbz++uuULFkSW1tb03V0//vf/564PMzTSruOMu16vUfZ2Njg6+sLwOHDhy22165d26Itrdh79H3p2rUrNjY2vPXWWwQEBLB06VLOnTv33NnTI9fcCSGEtdu5E7p3h/h4WLrUOHKXh9foehppM4eLFFlJXFyXbJ85nN37K126NCdPnrQYXUtP2uhR2nVmaUqVKpXhawOmyQ/Ozs7s37+fsWPHsm7dOtPIV4kSJXj//ff55JNPTNd+3bhxg5SUlAzXjAPjKF1iYiJFixY1tbm6umbYv1y5cvj5+bF161ZOnjyJt7c3165dY9OmTdSsWZOXXnrJrH9wcDCzZs2ibNmyvPHGG5QrV840IvjNN9881Xv2T9LemzJlyqS7Pa09vQkkRYoUsWhLe/9SUlJMbXXq1CEyMpIJEyawevVq/vvf/wLg5eXF2LFj6dq16/MdxCNk5E4IIazVgwfwwQfQrBlUrAjHjhlH7/J5YQc5P3M4u/fXqFEjALZt2/bEfikpKaZlUB4/fRgbG5vuc9Jmerq4uJjaypUrx8KFC7l27Rr/+9//mD17NsWLF2f8+PGMHz/e1M/FxYVixYqRmpqa4SMlJcXidPI/nWIMCAgA/h6tW758OSkpKab2NNeuXWP27NlUr16dU6dOsXTpUiZNmkRISAghISEULFjwift5WmnvTXqzYuHv0+aPvofPom7duqbTz3v27GHMmDHExsbyzjvvmE2+eF5S3AkhhDU6fhxefRVmzYIpUyA8HCpU0DuV1Vi5ci5BQb1xcHAgKKi32Uzi3Li/Hj16YGNjw5o1a/j1118z7Ldo0SKuXr2Kt7e3xWnY9Na+S0lJYffu3WiaRs2aNdN9zapVqzJo0CC2bt0KGK8/S1OvXj1u3rz5xEzPokOHDjg7O7N8+XKUUixZsgRbW1veeecds37nz59HKUWLFi1wdHQ023b58mXOnz9v8do2NjaA+ajZP/Hx8QFgx44dFtuSk5OJjIxE0zRTv+dla2tLvXr1GDdunOm6xoyWh3kWUtwJIYQ1SU2Fzz6D2rUhJQUOHDCO3v3/DyyRN3l4eDBq1CiSkpJo27YtJ06csOjz448/MmTIEAoUKEBoaKjF9u3bt7Nx40azti+++ILz58/TtGlT0+jar7/+mu4oX9qolYODg6ktODgYgL59+5pGrx51584doqKiMnGkRnZ2dnTp0oXLly/z2WefcezYMdM6eo9KW5IlMjLS7JrAxMRE+vbtm24BV7x4cQAuXrz41HneeustihUrRlhYmMXxfP7551y4cIHXX3/dbOJEZu3du5f76Vwrm977/rzkmjshhLAWv/8OPXrA9u0QHAwTJ0I6M+tE3vTpp59y584dPvvsM1555RVatmxJ1apVSUpKYu/evRw4cAAHBwfCwsLSnTzRpk0b2rdvT/v27alUqRJHjhxh06ZNFC9enHnz5pn6bdmyhQ8++ID69evj6emJq6srly9fZu3atdjY2PDBBx+Y+jZr1ozJkyfz8ccf4+npSevWrXF3dycxMZGLFy+ya9cuGjVqxE8//ZTp4w0ICGDhwoWMGjXK9PXjSpUqRdeuXVmxYgU1atSgefPmxMfHs3XrVhwcHKhRo4Zpokgab29vypUrx4oVK7C1taV8+fJomkb37t0pX758ulkcHR1ZtGgRnTt3pnHjxnTu3Bk3Nzeio6PZunUrZcqUYf78+Zk+xkdNnTqVHTt20KhRI9zd3XFycuKXX35h06ZNFCtWjH79+j3X65v5x0VVrJyscyeEyBO+/VapIkWUKldOqW3b9E7zVPRY5y4/OHDggAoICFAeHh7K3t5eOTs7q+rVq6sPPvhAXblyxaL/4sWLlcFgUEuWLFEbNmxQ9erVU46Ojqpo0aKqU6dO6syZM2b9T5w4oYYNG6Zq166tSpYsqQoVKqQ8PDxU586d1b59+9LNtHv3bvX222+rsmXLqoIFCypXV1dVs2ZNNXz4cIt/N3d3d+Xh4fFUx+rp6akMBoMqUaKESkpKSrfP3bt31ejRo1XlypWVnZ2dKl++vBo0aJC6ceOGatKkiTIYDBbPOXjwoPLz81MuLi7KYDAog8GgIiIilFLGde7S3q/0nte+fXtVsmRJVbBgQVWhQgX1/vvvq6tXr1r07dGjhzIYDOrixYsW29Jba2/Lli2qZ8+eqmrVqsrFxUU5Ojoqb29vNWTIEHXp0qWner+e9ntFUyp3L5YUExNDrVq1iI6OzrJz4UIIkWNu3YKBAyEsDLp2hXnz4JFZh9Ysqz5/5XNciKfztN8rclpWCCH0sn07BARAQgIsXw6PXUwuhBDPQiZUCCFETrt/H4YNAz8/8PQ0LnEihZ0QIotYZXEXExND27ZtKV68OI6OjlSvXp05c+boHUsIIZ7fsWNQpw7MnQvTp8O2bZDBRd5CCPEsrO607JYtW2jTpg21atUiJCQEJycnzp49y5UrV/SOJoQQzy4lBWbOhNGjwcsLDh2CDG7PJIQQz8Oqirvbt2/TvXt32rRpw+rVq/WOI4QQWePiReO1dbt2GU/HTpgAhQrpnUoIkUdZVXH37bffcu3aNSZMmAAYF0e0t7fHYLDKs8dCCPFkShknSgwcCC4uxgkUTZronUoIkcdZVdW0bds2ChcuzO+//46XlxfOzs64uLjw/vvv8+DBA73jCSHE07t507i0yXvvQZs2xmvtpLATQuQAqyruzpw5Q3JyMm+99RZvvPEGP/zwA7169eLLL7+kZ8+eescTQoins22b8Xq6LVtgxQpYtgyKFNE7lRAin7Cq07KJiYncvXuXAQMG8PnnnwPG+709fPiQ+fPnM378eCpXrqxzSiGEyMC9e/DxxzBrlnGZk8WL4TnuRSmEEM/Cqoo7e3t7ALp162bW3q1bN+bPn8/+/fszLO6Cg4NxcXGxeN7jryWEENni8GHw94dz54yzYgcPhjx0vXBYWBhhYWFmbfHx8TqlEUKD1JD1AAAgAElEQVQ8iVUVd2XLluXXX3+lVKlSZu2urq4A3Lp1K8Pnzpw5U25bI4TIeSkpMG0ahIRA1arGJU5eeknvVFkuvV+W026FJISwLlb1a2Xt2rUBuHz5sln7H3/8AUDJkiVzPJMQQmTowgVo2hRGjYLgYIiKypOFnRAid7Gq4u7tt98G4OuvvzZrX7hwIba2tjSRmWZCCGugFCxZAi+/DJcuwc6dMGWKrF0nhLAKVnVatkaNGvTq1YtFixaRnJyMr68vO3fuZPXq1YwaNYrSpUvrHVEIkd/duAH9+8P33xuXOZkzx7iGnRBCFwaDgcaNG7Njxw69o1gNqyruAL788kvKly/PN998w5o1a3B3d+fzzz9n8ODBekcTQuR3mzdDz57w4AGsWgWdO+udSOQRjy/WbzAYcHFx4eWXX6ZHjx4EBATolCx30DRN7whWxeqKuwIFChASEkJISIjeUYQQwujuXfjwQ/jiC2jRAr75BsqW1TuVyGM0TWPs2LEAJCUlcebMGdasWUNERASHDh1izpw5Oie0TidPnsTBwUHvGFbF6oo7IYSwKjEx8O67xskTs2cbbyWWh5Y4Edbl8YGNvXv34uvry7x58xg+fDju7u76BLNiVapU0TuC1ZFPKCGESE9KCkycCK+9Bvb2EB0NQUFS2IkcVb9+fby8vFBKERMTY7E9KiqKTp06Ubp0aQoVKkT58uUJDAzk6tWr6b7ewYMHadGihen2ns2bN2f//v18+umnGAwGdu3aZdbfYDDQtGlTYmNj6dOnD+XKlaNAgQIsWbLkmTKcP3+efv36UblyZRwcHChevDgvv/wyAwYM4ObNm6Z+Dx8+ZPbs2fj4+FCsWDEcHR3x8PDgrbfeIjw8PN2Mj4uPj+fjjz/Gy8sLe3t7ihUrRqtWrSyeD7Bz504MBgPjxo3jyJEjvPnmmxQpUgRHR0eaNGnCvn370n0/rZWM3AkhxON++804WWLfPuPp2E8/hYIF9U4l8jlbW1uzrxctWkS/fv2wt7enbdu2uLm5cfr0aRYuXMj69evZv38/bm5upv67du2iRYsWKKXo0KEDlSpV4tixYzRt2pRmzZpluN+bN29St25dnJ2d6dSpEwaDwTTBMTMZrl69Sp06dUhISODNN9+kc+fO3L9/n/Pnz7Ns2TKCgoIoVqwYAD169GDFihVUr16dgIAA7O3tuXLlCnv27GHz5s34+fmZZXz8mru4uDgaNGjAiRMnePXVV+nYsSN//fUXq1atokWLFoSGhtKvXz+LYz106BBTp06lfv369OvXj4sXL/L999/j5+fHkSNHcs0ooRR3QgiRRinjLcMGD4YSJSAiAho21DuVyMjdu3DyZM7v19sbcugar127dnHy5EkKFSrEq6++amo/ffo0gYGBVKxYkYiICMqUKWPatn37dlq0aMGQIUP44YcfAEhNTaV3794kJSXx008/0bJlS1P/+fPnM2DAgAwnJRw/fpzu3buzaNEis4kfmc2wevVqbt26xaxZswgKCjLbx71790z7j4+PZ8WKFdSuXZuoqCiLXI+O8GXkww8/5MSJE/Tv35/Q0FCz9tq1azN48GBatmxJhQoVzJ63ceNGFi9eTPfu3U1tCxYsIDAwkFmzZjF37tx/3Lc1kOJOCCEArl+Hfv1gzRrjjNjPP4fChfVOJZ7k5EnQ4w4Z0dGQDXdEUkoxbtw4lFIkJSVx9uxZ1qxZg6ZpTJ8+3ezuTaGhoSQnJzNr1iyzogqgWbNmtGnThvXr13Pnzh0cHR3Zu3cv586do1mzZmaFHUC/fv2YOXMmp0+fTjdXoUKFmD59usWM3sxmSGNnZ2exj7Tbj8Lfo3CFChVKt+BMG93LyMOHD1m2bBnOzs5MmjTJbFvlypUZPHgw//nPf1i6dCljxowx296wYUOzwg6gV69eDBw4kIMHDz5xv9ZEijshhPj5Z+jVC5KSjOvXdeigdyLxNLy9jYWWHvvNJuPGjTP7WtM0Fi1aZLEUSto1YDt37iQqKsrida5du0ZKSgqnT5+mZs2aHD58GDAWL4/TNI169eplWNy5u7tTokQJi/anzXDq1Cl8fHxo164do0ePZuDAgWzevJkWLVrQsGFDqlatava8woULmwrDGjVq0LFjRxo1asSrr776VLNiT506xb1792jYsCFFihSx2N6sWTP+85//cOTIEYttaXfKelSBAgUoVarUE2+Bam2kuBNC5F9378IHH8C8edCqFSxaBI+NQAgr5uCQLSNoetE0jZSUFMB4mnLv3r307t2b/v37U758ebNJAzdu3ABg2rRpT3y9xMREwHiqE7C4d3uajNqBDG8g8LQZ7ty5A0D58uU5cOAAn376KZs2bTKdrnVzc2PEiBFmp2pXrlzJlClT+Pbbb03Lw9jZ2dGpUyemT59uuud8etKO9fHRxMePJy4uzmJbesUgGAu8tH+b3ECmfQkh8qeDB6FmTeOadXPnwk8/SWEnrIa9vT1+fn6sX7+elJQUAgICuHfvnmm7i4sLmqZx+/ZtUlNT032kpKTQqFEjwDgaBhAbG5vu/jJqh4wXCM5sBgBvb29WrFjBjRs3OHToEJMnTyY1NZUhQ4awaNEiUz87OzvGjh3LqVOnuHTpEsuWLaNhw4YsW7aMTp06PfG9c/n/O8b8+eef6W5Pm8XrkofvLCPFnRAif0lOhn//G+rXN15Td/gwvP8+yAr3wgpVr16dvn37cvnyZWbOnGlqr1evHkopi6VLMuLz/yOckZGRFttSU1PZu3dvprNlNsOjbGxs8PHxYeTIkYSFhQGwdu3adPu+8MILvPPOO2zevJlKlSqxe/fuJ54i9fb2xt7enqNHj5pG8R6Vdpsynzw06vs4Ke6EEPnHuXPg62tc2uSjj2DvXvDy0juVEE/0ySefmCY1pJ1KHDRoELa2tgQHB3PmzBmL5zx8+NCskGvQoAGVKlVix44dbNq0yazvggULOHPmTKZv4ZXZDDExMekWW2kjbGnX012/fp3jx49b9EtMTCQxMRFbW1sKPmFpIltbW/z9/bl9+7bFhIlz584xe/ZsChYsyHvvvfd0B5oLyTV3Qoi8Tynj9XRDhkCpUhAZaRy5EyIXKFu2rGkpjqlTpzJx4kS8vLxYtGgRvXr1olq1arRq1QpPT0+SkpK4dOkSkZGRlCpVil9//RUwnlpduHAhrVq1om3btnTs2JGKFSty7Ngxtm3bxhtvvMHPP/9sMSP2STKbYenSpSxYsICGDRtSsWJFihYtyrlz51i/fj12dnYMHToUgMuXL+Pj40P16tWpXr06bm5u3L59mw0bNhAbG8uQIUPMZt+mZ/LkyURGRvLFF19w8OBBmjRpwvXr11m1ahV37tzhiy++sFgGJU9RuVx0dLQCVHR0tN5RhBDWKDZWqXbtlAKl+vRR6vZtvRPlGVn1+Suf40ppmqYMBkOG22NjY5Wjo6NycnJS165dM7UfP35c9ejRQ1WoUEEVKlRIFS9eXFWvXl0FBgaqHTt2WLxOVFSUat68uXJ2dlbOzs6qefPmav/+/WrgwIFK0zR19OhRi1xNmzZ9YvanzRAVFaUGDBigXnnlFVWsWDFlb2+vPD09Va9evdQvv/xi6hcXF6fGjx+vmjVrpsqVK6cKFSqkypYtq5o2bapWrFiR7nuXXsa4uDj14YcfKk9PT1WoUCFVtGhR1aJFC7V161aLvjt27FCapqlx48ale4zu7u7Kw8Pjie9DTnja7xVNKaX0LS+fT0xMDLVq1SI6OjpPnz8XQjyDjRuNS5ykpsLChdCund6J8pSs+vyVz3H9NWjQgIMHDxIfH2+25pywLk/7vSLX3Akh8p47dyAwEP71L6hTB/73PynsRL537969dJf/WLx4Mfv27aNFixZS2OURcs2dECJviYoCf3/44w/48kvjXSdkJqwQXLx4kZo1a9KiRQsqVapEcnIyhw8fZs+ePRQtWpQZM2boHVFkESnuhBB5Q1ISTJgA//mPcWHbjRshl9zkW4icULp0afz9/YmIiGDHjh08ePCAMmXK0KtXL0aPHo2Hh4feEUUWkeJOCJH7nTljHK2LjoZPPoHRo8HWVu9UQliVIkWK8NVXX+kdQ+QAKe6EELmXUrBgAQwbBmXLwp498NpreqcSQghdyYQKIYSulFIMHTqKTE/cj42Ftm2NEyf8/Y13mpDCTgghpLgTQugrOjqauXPnEBMT8/RPWrsWXnoJDhyAdetg/nxwcsq+kEIIkYtYXXG3c+dODAZDuo8DBw7oHU8IkcVCQ78jOXkGoaHf/XPnxETo2xfeegvq1YPjx6FNm+wPKYQQuYjVXnM3ZMgQ6tSpY9ZWqVIlndIIIbJSSMgkvvxyMS4uFUlIsAMms2FDBypXbsXt278RGNiD8eM/Nn/Svn3G06+xscbr7Pr0kSVOhBAiHVZb3DVq1IgOHTroHUMIkQ3GjBlByZKuTJjwI7GxawD+/882jBkzksDA7n93TkqC8eNh4kR49VXYvBkqV9YnuMhWJ06c0DuCEFbtab9HrLa4U0qRkJCAvb09BQpYbUwhxDOwtbUlKKg3X3yxltjYv9tdXDSCgnr/3XDq1N+TJT79FD7+GOTzIM/y9/fXO4IQeYLVfkr27NmTxMREbGxsaNSoEdOmTaNWrVp6xxJCZKHU1GTs7RdSpMhK4uK6kJqabNygFISGwogR4OZmPCX72GUaIu/w9vYmOjpa7xhC5Bre3t5P3G51xV2hQoXo1KkTrVu3pkSJEvzyyy9Mnz6dRo0asXfvXmrUqKF3RCFEFvHx8WDwYI3evdfy9ddh7N7tAVevQq9esGkTDBgA06aBo6PeUUU2cnBweOJN0IUQmaOpTC8ulfPOnTvHyy+/jK+vLz///LPZtpiYGGrVqkV0dLR8OAiR261ZY5wNW6AALFoErVvrnUg8gXz+CmGdrG7kLj2VKlWiXbt2/PDDDyil0NKZIRccHIyLi4tZW7du3ejWrVtOxRRCPKvbt2HoUPjmG+MyJwsWQMmSeqcSjwgLCyMsLMysLT4+Xqc0QognyRXFHcALL7zAw4cPuXPnDk7pLFY6c+ZM+c1RiNxozx7jpInr1+Hrr6FnT1nixAql98ty2sidEMK6WN0ixhk5f/489vb26RZ2Qohc6OFDGDUKfH2N94U9etR4rZ0UdkII8Vysrrj766+/LNqOHj3KunXraNGihQ6JhBBZ7sQJ4x0mpk0zrmEXEQEVK+qdSggh8gSrOy3bpUsXHBwcqFevHq6urvz6668sWLAAJycnJk+erHc8IcTzSE2FuXNh5Ehwd4f9+0FO6wkhRJayuuKuffv2LF++nJkzZ3L79m1cXV3p1KkTY8eOpaL8Zi9E7vXHH8br6bZsgUGDYMoUcHDQO5UQQuQ5VlfcBQUFERQUpHcMIURW+u476N8f7OyM69e1bKl3IiGEyLOs7po7IUQeEh8P3bvD229Ds2Zw/LgUdkIIkc2sbuROCJFH7NplLOxu3oTFi41/l5mwQgiR7WTkTgiRtR48gA8/hCZNjPeFPXoUAgKksBNCiBwiI3dCiKzzyy/w7rvw668wcSJ88AHY2OidSggh8hUZuRNCPL/UVJg1y7isycOHEBUFH30khZ0QQuhAijshxPO5csU4SWLoUAgMhOhoqFlT71RCCJFvyWlZIcSzW7XKWNDZ2xvXr2veXO9EQgiR78nInRAi8+Lj4b33oEsXY0F3/LgUdkIIYSVk5E4IkTkREcZlTeLi4L//NU6gkJmwQghhNWTkTgjxdB48MM5+bdoUPDzg2DHw97co7JRSDB06CqWUTkGFECJ/k+JOCPHPjh+HV181zoidMgXCw6FChXS7RkdHM3fuHGJiYnI4pBBCCJDiTgjxJKmp8NlnULs2pKTAwYP/uHZdaOh3JCfPIDT0uxwMKoQQIo0Ud0KI9P3+u3GSxPDhMHAgHDoEr7ySbteQkEm4unrh6fkGGzeeBvqyYcMpKlduhaurFyEhk3I2uxBC5GMyoUIIYSksDN5/H5ycjKdgmzV7YvcxY0ZQsqQrEyb8SGzsGoD//7MNY8aMJDCwew6EFkIIATJyJ4R41K1b8M47xscbbxgnTfxDYQdga2tLUFBvXFzMJ1e4uGgEBfXG1tY2uxILIYR4jIzcCSGMtm+HgABISIBvv4Vu3TL9EqmpydjbL6RIkZXExXUhNTU5G4IKIYR4Ehm5EyK/u38fhg0DPz/w9DTOjH2Gwg7Ax8eDKVM0zp5dy5QpGj4+HlkcVgghxD+RkTsh8rOjR41r1Z0+DTNmGO8Pa3j23/lWrpxr+ntQUG+CgnpnRUohhBCZICN3QuRHKSkwbZpx7TqDwTgTdtiw5yrshBBCWAf5JBciv7l40ThJ4sMPYcgQOHAAqlfXO5UQQogsYvXF3YQJEzAYDFSXHz5CPB+ljPeCffll+O034wSKqVOhUCG9kwkhhMhCVl3cXb58mYkTJ+Lo6IgmNyYX4tndvAldukD37tC2rXGJkyZN9E4lhBAiG1j1hIoRI0ZQv359kpOTuX79ut5xhMidtm0zLnFy7x6sXAlvv613IiGEENnIakfudu3axffff8/nn3+OUkpG7oTIrHv3jLNfmzeHqlWNS5xIYSeEEHmeVY7cpaSkEBQURN++falWrZrecYTIfY4cgXffhXPnYOZMGDxYZsIKIUQ+YZXF3ZdffsmlS5fYvn273lGEyF1SUmD6dBgzBqpVg+ho459CCCHyDav7Vf7GjRuEhIQQEhJC8eLF9Y4jRO5x4QI0bQoff2xcsy4qSgo7IYTIh6xu5O6TTz6hRIkSBAUFZep5wcHBuLi4mLV169aNbs94GyUhco20JU4GDYJixSAiAho10juVyGPCwsIICwsza4uPj9cpjRDiSayquDtz5gxfffUVn3/+OZcvXza1379/n4cPH3Lx4kUKFy5M0aJFLZ47c+ZMfHx8cjKuEPq7cQP694fvvzcuczJnDhQurHcqkQel98tyTEwMtWrV0imRECIjVlXcXblyhdTUVAYPHszgwYMttnt4eDB06FA+++wzHdIJYWU2b4aePeHBA1i9Gjp21DuREEIIK2BVxV316tVZs2aN2bInSik++eQTEhMTmTVrFpUqVdIxoRBW4O5d463DvvgCWraERYugbFm9UwkhhLASVlXcFS9enHbt2lm0z5w5E4C2bdvmdCQhdKGUIjh4NDNnTjBf4zE6Gvz9jZMn5syBgQNB1oAUQgjxCKubLZseTdNkEWORr0RHRzN37hxiYmKMDSkpMGEC1K0LDg4QE2OcQCHfF0IIIR6TK4q7HTt2cOzYMb1jCJFjQkO/Izl5BqGh38H58+DrCyEhMHIk7NsHL76od0QhhBBWKlcUd0LkByEhk3B19cLT8w02bjwN9MHpuy0kVq7Cxf0HWdi9r3H0rmBBvaMKIYSwYlZ1zZ0Q+dmYMSMoWdKVCRN+JDn2K76nEx1uHybMzo34cSPpHdxf74hCCCFyARm5E8JK2NraEhTUmzY21/gfL+HLLjrwPZ+Wr0HgyEHY2trqHVEIIUQuIMWdENbi7l14/32++uMARw0led21Kpvsb5Kamqx3MiGEELmIFHdCWIODB6FmTVi8mK99fDk9cwh7f9vElCkaPj4eeqcTQgiRi8g1d0LoKTkZJk2C8eOhRg04fJjeXl6mzUFBvQkK6q1jQCGEELlNpoq7jh074uHhgZeXFy+++CINGzbMrlxC5H3nzsF770FUFIweDWPGgFxXJ4QQ4jllqrgLDw9n3bp1lC5dmtjY2OzKJETephR8/TUMHQqlS8Pu3VCvnt6phBBC5BGZKu6aNm2Kr68vAFWqVMmWQELkadeuQd++sG4d9OkDM2eCk5PeqYQQQuQhmSrunOSHkBDPbsMG6N3bOHK3di3IvZKFEEJkg0zNlv3tt9+4fv36U/VNSUl5pkBC5Dl37kD//tCmDdSpA8ePS2EnhBAi22SquNu7dy+urq5Uq1aNwMBAli9fzqVLl9Lt26tXrywJKESuFhVlnAW7bBl8+SWsXw+lSumdSgghRB6WqdOyvr6++Pv7ExkZyZYtW1iwYAEAbm5uNGrUyPSoWrUqCQkJ2RJYiFwhKcl4H9j//Adq1YKNG0GuUxVCCJEDMlXcubm50adPH/r06QPAlStXiIyMZPfu3URGRrJixQpSU1MpWrQo9+7dy5bAQli9M2fA3x+io43Lm4weDQVkSUkhhBA5I1M/cW7evGn2dbly5ejatStdu3YFIC4ujj179hAREcHcuXOzLqUQuYFSsGABDBsGZcvCnj3w2mt6pxJCCJHPZKq4i4iIIDY2llIZXDNUpEgR3nzzTd58801Onz6dJQGFyBViY41Lm2zYAP36wYwZssSJEEIIXWRqQsXAgQPx9/dn2bJl3L9//4l9HR0dnyuYELnGunVQvTocOGCcMDF/vhR2QgghdJOp4m7KlCls3bqVF198kXnz5j2x7+TJk58rmBBWLzHROFrXrp3xDhPHj8O//qV3KiGEEPncM13lXatWLWrVqvXEPm5ubs8USIhcYd8+431h//wTvvrKuDixpumdSgghhMjcyJ0Q+V5SknEGbMOGULIkHDliHL2Twk4IIYSVsLri7pdffqFz585UqlQJR0dHSpYsSePGjdmwYYPe0UR+d+oU1K8PkybBuHEQGQmVK+udSgghhDBjdYtvXbp0icTERHr06EHZsmW5e/cuq1evpm3btsyfP5++ffvqHVHkN0pBaCiMGAFubsZTsnXq6J1KCCGESJemlFJ6h/gnqamp1KpVi/v373PixAmzbTExMdSqVYvo6Gh8fHx0SijyrKtXjdfT/fwzvP8+TJsGDg56pxLCKsjnrxDWyepOy6bHYDDwwgsvEB8fr3cUkZ+sWWNc4uTwYePtw+bOlcJOCCGE1bO607Jp7t69y927d4mPj2fdunVs2rTJdCcMIbJVQgIMGQLffAPt2xvvOlGihN6phBBCiKditcXdsGHDWLBgAWAcuevYsSNffPGFzqlEnrdnj3GJk7/+gkWLoEcPmQkrhBAiV7Ha4i44OJi3336bK1eusGrVKpKTk3nw4IHesURe9fChcQbs5MnGBYm3bYOKFfVOJYQQQmRarphQAdCyZUvi4uKIiooya0+7oNfX1xcXFxezbd26daNbt245GVPkRidOgL8/HDtmLPA+/BBsbPROJYRVCQsLIywszKwtPj6eXbt2yYQKIayM1Y7cPa5jx44EBgZy5swZPD09LbbPnDlTPlxE5qSmGidJjBwJHh4QFQXyf0iIdKX3y3LaL9dCCOuSK2bLAty7dw9AZsyKrPHHH/DGGzB4MPTtC9HRUtgJIYTIE6yuuPvrr78s2pKSkli6dCkODg5UrVpVh1QiT1m92rjEyfHjsGkTzJ4N9vZ6pxJCCCGyhNWdlu3Xrx8JCQn4+vpStmxZ/vzzT5YvX87p06eZMWMGDrLOmHhW8fHGkbqlS6FjR5g/H4oX1zuVEEIIkaWsrrjr2rUrX3/9NaGhody4cQNnZ2dq167NtGnT+Ne//qV3PJFb7doF3bvDzZuwZIlxuRNZ4kQIIUQeZHXFXZcuXejSpYveMURe8eABhIQYbxvWsCHs3Anu7nqnEkIIIbKN1RV3QmSZX34xLnHyyy8waRKMGCFLnAghhMjzpLgTeU9qKsyZg/rwQ/50dKb0/v1oMhNWCCFEPmF1s2WFeC6XL0PLljB0KNc6dKBK/D1icsc63UIIIUSWkOJO5B0rVxqXODlxArZuZZS9G4kpnxEa+p3eyYQQQogcI8WdyP3i4ozX1nXtyvGyblS5b4fngBls3Hga6MuGDaeoXLkVrq5ehIRM0jutEEIIka2kuBO5286d8PLLsH49LFuG9+FDBI39mISEAsTGrgE0YmPXkJhoy5gxIxkzZoTeiYUQQohsJcWdyJ0ePIAPPoBmzaBSJePdJt59F9uCBQkK6o2Li/kadi4uGkFBvbG1tdUpsBBCCJEzZLasyH2OHzeehj15EqZOhWHDwGD+e0pqajL29gspUmQlcXFdSE1N1imsEEIIkbNk5E7kHqmp8NlnULu28e8HDxrXrjNY/jf28fFgyhSNs2fXMmWKho+Phw6BhRBCiJwnI3cid/j9dwgIgB07jCN1EyaAnV2G3VeunGv6e1BQb4KCeudESiGEEEJ3UtwJ6/ftt/D+++DsDOHhxuvshBBCCJEuOS0rrNetW9CtG7z7LrRuDceOSWEnhBBC/AMZuRPWaft242nYhATjyF23bnonEkIIIXIFGbkT1uX+feM1dX5+UKWKcWasFHZCCCHEU5ORO2E9jh41noI9e9Y4K3bIkHRnwgohhBAiY/KTU+gvJcW4Xl2dOmBjA4cOQXCwFHZCCCHEM5CfnkJfFy8aJ0l89BEMHQoHDsBLL+mdSgghhMi15LSs0IdSsGwZDBoERYoY169r3FjvVEIIIUSuJyN3IufdvAldukD37tCunfFaOynshBBCiCwhI3ciZ23dCj16wL17sGoVdO6sdyIhhBAiT7G6kbuDBw8yaNAgqlWrhpOTExUqVKBLly6cOXNG72jiedy7Z5z92qIFVKtmXOKkc2eUUgwdOgqllN4JhRBCiDzB6oq7KVOmsGbNGpo3b87s2bPp168fu3btwsfHh19++UXveOJZHD4MtWvD/PkwaxZs2gTlygEQHR3N3LlziImJ0TmkEEIIkTdY3WnZ4cOHU6dOHQoU+Dtaly5dqF69OpMnT+a///2vjulEpqQtcTJ2rHG0Ljra+OcjQkO/Izl5BqGh37FwYS2dggohhBB5h9WN3NWrV8+ssAOoXLkyVatW5eTJkzqlEpn222/QpAmMHg3Dh0NUlKmwCwmZhKurF56eb7Bx42mgLxs2nKJy5Va4unoREjJJ1+hCCCFEbmZ1xV16lFLExsZSokQJvaOIf6IULFkCr7wCly9DRARMmgQFC5q6jBkzgjFjRkSLWPcAABZsSURBVJKQUIDY2DWARmzsGhITbRkzZiRjxozQL78QQgiRy+WK4m758uX88ccfdOnSRe8o4kmuXzfOfu3RAzp2NC5x0qiRRTdbW1uCgnrj4qKZtbu4aAQF9cbW1jaHAgshhBB5j9Vdc/e4kydPMnDgQOrXr09AQIDecURGNm2Cnj3h4UNYvdpY3P2D1NRk7O0XUqTISuLiupCampwDQYUQQoi8zaqLuz///JM333yTokWLsnr1ajRNy7BvcHAwLi4uZm3dunWjW7du2R0zf7t7F0aOhLlzoWVLWLQIypZ9qqf6+HgweLBG795r+frrMHbv9sjmsEKIZxUWFkZYWJhZW3x8vE5phBBPoikrXWAsPj6eJk2acPnyZSIjI/H29k63X0xMDLVq1SI6OhofH58cTpnPRUeDvz9cuADTp8P778MTCnAhRN4in79CWCervObu/v37tGnThrNnz7Jhw4YMCzuhk+RkmDAB6tYFR0fjOnYDB0phJ4QQQlgBqzstm5KSQpcuXYiKimLt2rW89tprekcSjzp/Ht57D/bvh48+Mq5h98hMWCGEEELoy+qKu+HDh7N+/XratGnD9evXWbZsmdl2f39/nZLlc0rBN98YbyFWsiTs2gUNGuidSgghhBCPsbri7ujRo2iaxvr161m/fr3ZNk3TpLjTw19/Qb9+8OOP0KsXfP45ODvrnUoIIYQQ6bC64m7Hjh16RxCP+uknY0GXkgI//ADt2+udSAghhBBPYJUTKoQVuHMHBgyAN98EHx84flwKOyGEECIXsLqRO2EFDh40LnHy++8wbx4EBspMWCGEECKXkJE78bfkZBg/HurVAxcX4xInAwZIYSeEEELkIjJyJ4zOnjUucXLgAHzyifEh93gVQgghch0p7vI7pWDhQggOhtKlYc8e4+LEQgghhMiV5LRsfnbtGrz1lnGZk3fegSNHpLATQgghcjkZucuv1q+H3r2Nf1+7lv9r796DqqwTP46/QRBBlFlREmxWCm+bYXiXrcyx9RKt5ChqWpKG2i8TSy1NN3BacpXSpbxkum5KmuaK2aptjmNKN68Lamaa2oqsitcRBBW5Pb8/npXNpNb0nPOc85zPa4YZ/Z4483mGRj7n+zzf75f4eGvziIiIiENo5s7blJTAM8+YZa5zZ3OLExU7ERER29DMnTfZvt1cNHHyJCxYACNHaiWsiIiIzWjmzhuUl8PUqfDAAxAaaj5bN2qUip2IiIgNaebO7g4dMjckzs2FlBT4wx/ATz92ERERu9JvebsyDPPW64QJ0KQJbN0KnTpZnUpEREScTLdl7ejUKejTxzxdIjHRPGlCxU5ERMQraObObv7+dxgxAnx9Yf16ePRRqxOJiIiIC2nmzi6Ki81S17cv3H8/fPONip2IiIgX0sydHWzdam5xcvq0eZTY009rJayIiIiX0sydJysvh1degQcfhDvugL17zVMnVOxERES8lmbuPNXBg+YWJ3v3wquvwssva4sTERER0cydxzEMmDcP2rUzjxLbts2cvVOxExEREVTuPEtBAcTFwZgxMHy4uTFxhw4/+Z8bhsELL0zBMAwXhhQREREruWW5u3TpElOnTqV37940aNAAX19fMjMzrY5lrQ8/hOhoc8+6jz82Z++Cgn72W3Jycpg3bw65ubkuCikiIiJWc8tyd/bsWdLS0vjuu++IiYkBwMdbFwlcvGjO0vXvby6c2LfPnL27CfPnr6KiYhbz569yckgRERFxF25Z7iIiIjh16hRHjx7ljTfesDqOdb78Eu67D7Ky4N13zdm7Ro1+9ltSU6cTFtaS5s0f4eOPDwEjWb/+O5o1601YWEtSU6e7JruIiIhYwi3LXe3atQkLCwPwzufFyspg8mTo2tU8F3bvXnP27iZmL1NSXiQlZSLFxX6cPr0G8OH06TWUlPiTkjKRlJQXnZ9fRERELOOW5c6rHTgAsbEwcya89hp89hncffdNf7u/vz/JyUmEhFxfBENCfEhOTsLf39/RiUVERMSNqNy5i6oqmD3b3OLk8mXYvh2mTIFatW7x7SoIDFxEeHgPAgMXUVVV4eDAIiIi4o5U7tzBiRPQuzc8/zyMHAk5OdC+/W29Zbt2d5Ge7sORI38nPd2Hdu3uclBYERERcWe22fl23LhxhISEXDc2ePBgBg8ebFGim7RqFTzzDNSpAxs2QK9eDnnblSvnVf85OTmJ5OQkh7yviHinFStWsGLFiuvGioqKLEojIj/HNuUuIyODdu3aWR3j5hUVQXIyLF0KCQnwzjsQGmp1KhGRGtX0YTk3N5f2t3mXQUQczzblzqN8/jkMHQoXLsCSJZCYeFMrYUVERET+Fz1z50pXr8KkSdCtG/z61/D11/DUUyp2IiIi4jBuO3M3d+5cCgsLOXnyJABr164lPz8fgLFjx1K/fn0r4/1y+/fDE0/At9/C9Onw4ou3vBJWRERE5Ke4bbmbNWsWx44dA8yjx9asWcOHH36Ij48PiYmJnlPurm1x8vLLEBUFO3ZA27ZWpxIRERGbcttyd/ToUasj3L7jx2HYMPj0U3Obk+nTITDQ6lQiIiJiY25b7jzeBx/As89CUBBs3Ag9elidSERERLyAFlQ4WmGh+Wzd4MHQsyfs26diJyIiIi6jmTtH2rLFXP168SIsWwZDhmglrIiIiLiUZu4cobTUXP368MNw993mFidPPKFiJyIiIi6nmbvb9fXX8OST8N138PrrMH48+Kozi4iIiDXUQm5VVRXMnAkdO4JhwK5d5uydip2IiIhYSE3kVuTnw+9+By+9BGPGmMWuTRurU4mIiIjotuwvtnw5jB4N9eqZ+9d17251IhEREZFqmrm7WRcumNubPPEExMWZz9qp2ImIiIib0czdzfj0U3OLk0uXzJm7wYOtTiQiIiJSI83c/ZzSUhg3zny+rmVLc7ZOxU5ERETcmGbufsqePeYWJ0eOwJ//bJ4Nq5WwIiIi4ubUVn6sstLcr65TJ6hVC/75T3P2TsVOREREPIBXNhbDMHjhhSkYhnH9C8eOmYskXn4ZXngBdu6Ee++1JqSIiIjILfDKcpeTk8O8eXPIzc01BwwDli4196rLyzPPiH39dQgIsDSniIiIyC/lleVu/vxVVFTMYv78VXD+PAwcCImJEB9vLpp46CGrI4qIiIjcEq8pd6mp0wkLa0nz5o/w8ceHgJEUr/6C03c04cLqD1nZ73Fz9i4kxOqoIiIiIrfMa8pdSsqLpKRMpLjYj6LTy3mL51lZuJUDteqx5tV0+n3wntURRURERG6b15Q7f39/kpOTiA0oIof2jGIhY3mL/2vahadTXsTf39/qiCIiIiK3zWvKHZWVMGMGK/O/pNynhN4N27IoMIhKo9LqZCIiIiIO43bl7urVq0yaNImIiAiCgoLo0qULmzZtur03PXoUunWDKVP4R6sYvpz1B/5x7FPS031o1+4uh+QWERERcQduV+6GDRtGRkYGQ4cOZfbs2dSqVYu4uDi++uqrX/5mhgFLlsB998Hx4/DZZ/Q9kMtz454hKCiI5OQkVq6c5/BrEBEREbGKW5W7nTt3snLlSmbMmEF6ejojRoxg8+bNNG3alIkTJ/6yNzt3DhISYPhw6N8f9u6FBx90TnAXWrFihdURXMZbrtVbrhN0rSIiruBW5S4rKws/Pz9GjRpVPRYQEEBSUhLbtm3jxIkTN/dGGzaYGxJ/9hlkZcHixVC/vpNSu5Y3/cLwlmv1lusEXauIiCu4VbnbvXs3LVq0IDg4+Lrxjh07ArBnz56ff4PLl2HMGHjkEfNW7L595qydiIiIiJfwszrADxUUFBAeHn7D+LWxkydP/vQ3HzgAQ4aY58POnQujR4OPj7OiioiIiLgltyp3V65cIaCG81zr1KlT/fpPeuopiImB3buhVStnRRQRERFxa25V7gIDA7l69eoN46WlpdWv/9i1wnfg0UdhyhTz1mxurnODWqioqIhcG1/fD3nLtXrLdYKu1W4OHDgA/I8P3iLicm5V7sLDw2u89VpQUABARETEDa/l5eUB8OTatbB2rVPzuYv27dtbHcFlvOVaveU6QddqR3l5edx///1WxxCR/3Crcte2bVuys7MpLi6mXr161eM7duwAICYm5obv6dWrF8uWLSMyMrLGmT0REXGOK1eukJeXR69evayOIiI/4GMYhmF1iGt27txJly5deOONN5gwYQJgnlhx77330qhRI7Zu3WpxQhERERH35lYzd506dWLAgAFMnjyZM2fOEBUVRWZmJvn5+SxevNjqeCIiIiJuz61m7sCcqUtJSWHZsmVcuHCB++67j7S0NHr06GF1NBERERG353blTkRERERunVudUCEiIiIit8djy93Vq1eZNGkSERERBAUF0aVLFzZt2mR1LIe7dOkSU6dOpXfv3jRo0ABfX18yMzOtjuUUu3btYsyYMbRu3Zrg4GCaNm3KoEGDOHz4sNXRHGr//v0MGDCAqKgo6tatS6NGjXjooYdYv3691dFcYtq0afj6+hIdHW11FIfKzs7G19e3xq+dO3daHc/hcnNziY+PJzQ0lLp16xIdHc2cOXOsjiUiuNmCil9i2LBhrF69mnHjxtG8eXMWL15MXFwcW7ZssdV+S2fPniUtLY2mTZsSExNDdnY2PjY9Vi09PZ1t27YxYMAA2rRpQ0FBAXPnzqVdu3Zs376d1q1bWx3RIfLz8ykpKWHYsGFERERw+fJlsrKyiI+PZ8GCBYwcOdLqiE5z/Phx/vSnP1G3bl3b/n/8/PPPV5+HfU1UVJRFaZxj48aN9OnTh/bt25OamkpwcDBHjhzhxIkTVkcTETz0mbtrW6bMnDmT8ePHA//dMiUsLIyvvvrK4oSOU1ZWRmFhIWFhYeTk5NCxY0eWLFlCYmKi1dEcbtu2bXTs2BE/v/9+5jhy5AjR0dEkJCSwdOlSC9M5V1VVFe3bt6e0tLR61387evzxxzl//jwVFRWcO3eOffv2WR3JYbKzs+nevTtZWVn069fP6jhOc/HiRVq0aMEDDzxAVlaW1XFEpAYeeVs2KysLPz8/Ro0aVT0WEBBAUlIS27Zts9Wnx9q1axMWFgaAB/bwXyQ2Nva6YgfQrFkz7rnnHg4ePGhRKtfw9fXlzjvvpKioyOooTvP555+zevVq3nzzTQzDsO3MnWEYFBcXU1FRYXUUp1i+fDlnzpxh2rRpgPnoSFVVlcWpROSHPLLc7d69mxYtWhAcHHzd+LVbIXv27LEiljiBYRicPn2ahg0bWh3F4S5fvsy5c+f4/vvvycjIYMOGDTz88MNWx3KKyspKkpOTGTlypG1ur/+U4cOHExISQmBgIN27dycnJ8fqSA61adMm6tevz7///W9atmxJvXr1CAkJYfTo0TWeDS4irueRz9wVFBQQHh5+w/i1sZrOpxXP9P7773Py5Elee+01q6M43Pjx41m4cCFgztz179+fuXPnWpzKOd555x3y8/PZvHmz1VGcJiAggISEBOLi4mjYsCH79+9n5syZPPjgg2zdurXG4xM90eHDh6moqKBv376MGDGC9PR0tmzZwpw5cygsLGT58uVWRxTxeh5Z7q5cuUJAQMAN43Xq1Kl+XTzfwYMHee655/jtb3/LU089ZXUchxs3bhwDBw7kxIkT/O1vf6OiosKWMx/nz58nNTWV1NRUQkNDrY7jNLGxscTGxlb//fe//z0JCQm0adOGyZMn88knn1iYznFKSkq4fPkyzz77LG+++SYAffv2paysjAULFvDHP/6RZs2aWZxSxLt55G3ZwMDAGn8JlpaWVr8unu3UqVM8+uij/OpXvyIrK8uWz2e1bNmS7t27M3ToUNatW0dJSQl9+vSxOpbDvfLKKzRs2JDk5GSro7hcVFQUjz32GFu2bLHNM7PX/n0dPHjwdePX/r59+3aXZxKR63lkuQsPD6/x1mtBQQEAERERro4kDlRUVMQjjzzCxYsX2bBhA40bN7Y6kkv079+fXbt22Wpfv8OHD/OXv/yF5ORkjh8/Tl5eHnl5eZSWllJWVsaxY8e4cOGC1TGd6s4776SsrIxLly5ZHcUhrv37escdd1w3fm3hl91/niKewCPLXdu2bTl06BDFxcXXje/YsQPANs+2eKPS0lL69OnDkSNHWL9+Pa1atbI6kstce5zATitmT5w4QVVVFWPHjuXuu++u/tq5cyeHDh3irrvuIi0tzeqYTvWvf/2LwMDAGxaAeaoOHToA5p6FP3TtA3ejRo1cnklErueR5S4hIYHKysrqh9HB3Odu8eLFdOnShSZNmliYTm5VZWUlgwYNYseOHaxatYrOnTtbHckpzp49e8NYeXk57733HkFBQdxzzz0WpHKO6Oho1qxZw0cffVT9tWbNGlq3bk3Tpk356KOPSEpKsjqmQ9T0c927dy9r166lZ8+eFiRyjoEDBwLw17/+9brxRYsW4e/vT7du3SxIJSI/5JELKjp16sSAAQOYPHkyZ86cISoqiszMTPLz81m8eLHV8Rxu7ty5FBYWVn8yXrt2Lfn5+QCMHTuW+vXrWxnPYSZMmMC6devo06cP586dY9myZde9/uSTT1qUzLFGjRpFcXExXbt2JSIiglOnTvH+++9z6NAhZs2aRVBQkNURHSY0NJTHHnvshvGMjAwA4uPjXR3JaQYNGkRQUBCxsbGEhYXx7bffsnDhQoKDg5kxY4bV8RwmJiaGp59+mnfffZeKigq6du1KdnY2WVlZTJkyxWseoxBxa4aHKi0tNV566SUjPDzcqFOnjtG5c2dj48aNVsdyisjISMPHx8fw8fExfH19DV9f3+o/Hzt2zOp4DtOtW7fqa/vxl6+vr9XxHOaDDz4wevToYTRu3Njw9/c3GjRoYPTs2dNYt26d1dFcplu3bkZ0dLTVMRxq9uzZRufOnY3Q0FDD39/faNKkiZGYmGh8//33VkdzuPLycuPVV181IiMjjdq1axstWrQw3nrrLatjich/eOTxYyIiIiJSM4985k5EREREaqZyJyIiImIjKnciIiIiNqJyJyIiImIjKnciIiIiNqJyJyIiImIjKnciIiIiNqJyJyIiImIjKnciIiIiNqJyJyIiImIjKnciIiIiNqJyJyIiImIjKnciIiIiNqJyJyIiImIjKnciIiIiNqJyJyIiImIjflYHEPFG5eXlZGRk4OPjQ05ODtOmTSMzM5OSkhI6dOjAkCFDrI4oIiIeSuVOxAJvv/02AwcOJDIykokTJ9KrVy/27dtHjx49OHnypMqdiIjcMpU7EQv4+fkRGRkJwLlz54iPjycwMJC3336bJk2aWBtOREQ8mo9hGIbVIUS82W9+8xumTZtGv379rI4iIiI2oHInYqEzZ87QuHFjTp06RVhYmNVxRETEBrRaVsTFysvL2bx5MwBffPEFkZGR1cVuw4YNHDx40Mp4IiLi4VTuRFxs0aJFxMXFceXKFT755JPqYldWVsbmzZtp1aqVxQlFRMST6basiIt98803pKWl0bJlS/r378/ChQtp0KABhmEwduxY3Z4VEZHbonInIiIiYiO6LSsiIiJiIyp3IiIiIjaiciciIiJiIyp3IiIiIjaiciciIiJiIyp3IiIiIjaiciciIiJiIyp3IiIiIjaiciciIiJiIyp3IiIiIjaiciciIiJiIyp3IiIiIjby/2euKut5kqahAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x32324b910>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated slope = 1.204958783179411\n"
     ]
    }
   ],
   "source": [
    "using PyPlot\n",
    "\n",
    "# Compute synthetic coordinates (xi,ti) and then use them to estimate the slope of \n",
    "# of the regression line\n",
    "\n",
    "xi = [0.22, 1.2, 1.9, 2.3, 3.1, 3.9, 4.1, 5.4]\n",
    "ti = 1.2*xi + 0.3*randn(size(xi))\n",
    "\n",
    "# Estimate the slope \n",
    "\n",
    "a_estimated = fit_a_line(xi,ti)\n",
    "\n",
    "# Use estimated slope to compute a regression line \n",
    "\n",
    "x = collect(0:0.1:maximum(xi))\n",
    "\n",
    "t = a_estimated*x;\n",
    "\n",
    "# Plot observations and regression line \n",
    "\n",
    "subplot(221);plot(xi,ti,\"*b\",label=\"Observations\"); \n",
    "subplot(221);plot(x ,t ,\"-r\",label=\"Regression\"); \n",
    "\n",
    "title(L\"Regression with $b=0$\"); \n",
    "xlabel(L\"$x$\"); \n",
    "ylabel(L\"$T$\")\n",
    "legend(bbox_to_anchor=[1.05,1],loc=2,borderaxespad=0)\n",
    "println(\"Estimated slope = \",a_estimated)\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAFKCAYAAAB/xnmMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlYVGX7wPHvGUB2UVxwSQUVIY1K1HLPLTX7uZtKWbiLCypqVia4lGsWiilqvKZmoWaZW7nhhruCZuWeaWKGuSG4hMDz+2Ne5nWcQUWBGeD+XNdcMs95zjn3GQVvnlVTSimEEEIIIUSBoLN0AEIIIYQQIudIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcieEEEIIUYBIcidENp07dw6dTkevXr0sHUqOeZJnKoifgxBCFASS3BVSOp3O6GVra0uJEiVo2rQpixcvtnR4+YKmaZYOIcc9+EyPk8DlxOewa9cu6tSpQ6lSpejUqdNTX+/ChQvMmjXroXVOnjxJ586dGTFiBCNHjiQwMJDLly8/9b2FEMLSbC0dgLAcTdMYN24cAPfu3eP06dOsWrWKHTt2cOjQIWbPnm3hCK3TM888w4kTJ3Bzc7N0KDkmq2fKTNxyO5Ft2LAhu3fvxtnZmbp16z719cLDw3F2ds7y+M2bN2nRogXTpk3jzTffBGDKlCm0aNGCuLg47OzsnjoGIYSwFEnuCrmwsDCj93v27KFx48bMnTuXkSNH4unpaZnArJitrS3VqlWzdBg5KqtnUkoZ/Zmb4uPjSU9Pp0mTJk91nZSUFJYuXcrRo0ezrDN9+nTu3btH9+7dDWVBQUGMGzeOhQsXMmDAgKeKQQghLEm6ZYWR+vXr4+Pjg1KK+Ph4k+P79++nS5culClTBnt7eypWrEhQUBCXLl0yez2lFLNmzaJ69eo4OjryzDPPEBwcTFJSEp6ennh5eRnq3t8FeOrUKbp160bp0qWxsbFhx44dTxzHmjVraN68OWXLlsXBwYHy5cvTpEkTIiMjn6jew7oqV6xYQePGjXFzc8PJyYnnn3+eqVOnkpqaalL3/uucO3eO7t27U7JkSRwdHalTpw7r1683+5k+KCUlhSJFitCwYUOj8jt37uDg4IBOp2Pp0qVGxyIjI9HpdCxatCjLZxo/fjyVK1cGYPHixUbd+EuWLDF5lieNP9OOHTtwdXWldu3a2TrvQV9++SWtW7emTJkyWdZZsWIFL7/8Mjrd/34EFi9eHF9fX7799tunur8QQliatNyJLD3YNbVw4UL69++Po6Mj7dq1o0KFCpw6dYqoqCjWrl3Lvn37qFChgtE5gwcPZt68eZQvX54BAwZgZ2fHmjVrOHDgAGlpaRQpUsTkvr///jt169bFx8eHt99+mzt37hh1F2YnjgULFhAUFETZsmVp3749JUuW5PLly/z8888sWrSIgQMHZqve/R7sqhwzZgxTp06lVKlS9OjRAxcXF3788UfGjBnDxo0b2bRpk9nuvvPnz/Pyyy9TpUoVAgMDuXr1KsuXL6d9+/Zs2bLlkS1ZLi4uvPzyyxw4cICUlBRcXFwA2L17tyGpjImJoUePHoZzYmJi0DSN5s2bZ/lMTZs2JSkpiVmzZvHiiy/SoUMHw7EXX3zR8PW5c+eeKv5MO3bsoEGDBkYJV3Yppfj888+Jjo7Osk5ycjJnzpyhVatWJsfKli3LgQMHnvj+QghhFZQolDRNUzqdzqR8x44dSqfTKQcHB/X3338byk+ePKns7OyUt7e3+uuvv4zOiYmJUTY2Nqpjx45G5Tt37lSapilfX1+VlJRkKE9NTVWNGzdWmqYpLy8vQ/kff/yhNE1TmqapDz/80Gzc2Y3D399fOTg4qH/++cfkWlevXs12vfvj7NWrl6Fsz549StM0ValSJZWYmGgoT0tLU23btlWapqnJkyebvY6maWrixIlGxzZu3Kg0TVNt2rQx+zk8KCwsTGmaptavX28oe//995Wtra1q3ry5qlChgqE8PT1dubu7q6pVqz70mZRS6ty5c2bLczr+tLQ0VbRoURUUFKRGjRqlRo4cqVq1aqUOHz78WOdnWrt2rWrYsOFD6/zyyy9K0zQ1ZswYk2MdO3ZUOp1O3bt3L1v3FUIIayLJXSGV+Z/y+PHj1bhx49SYMWNU165dlZ2dnbKxsVGff/65Uf3hw4crTdPUjz/+aPZ6HTp0ULa2tiolJcVQ1qdPH6Vpmvrqq69M6u/evTvL5K5s2bIqNTXV7H2yG4e/v79ydnZW169ff+jn8bj17o/z/oSnb9++StM09cUXX5jUP3XqlLKxsVGVK1c2ex0vLy+VkZFhcl7FihVVqVKlHhmPUvqkXNM0NWLECENZnTp1VN26ddWcOXOUpmnq1KlTSiml4uLilKZpasCAAQ99poeV53T8Bw4cUJqmqddff12lpaUppZT64osvlKenp+F9pjt37hj+/T7Oq2nTpkbnZ/7bGz9+vEkcPXr0UJqmmU3yhRAiv5Bu2UJuwoQJRu81TWPhwoUEBgYale/duxeA7du3s3//fpPrXL58mfT0dE6ePIm/vz8Ahw8fRtM0k7FgAC+//DI2NjZmY3rhhReynK34uHGcOnWKmjVr0qNHD0aOHEn16tXp3r07jRs3pkGDBpQqVcrovMetl5X4+Hg0TaNZs2Ymx7y9vSlfvjznzp0jOTkZV1dXo+Mvvvii2dmoFSpUMPuM5tSrVw9HR0diYmIASEpK4vDhw7z33nuGmLZu3Yq3tzdbt24FMBvrk8iJ+Hfs2IGbmxvffvut4d9FhQoVOH/+PPHx8dSpU8dQ18HBgRMnTphcY/fu3Xz66ad8//33RuVOTk5G7zOvby7me/fuAZCWlvZYcQshhDWS5K4Q0zSN9PR0QD/4fs+ePfTp04cBAwZQsWJFmjZtaqh79epVAD755JOHXu/WrVuG90lJSQB4eHiY1LWxsaFEiRJmr/OwgfCPG0dKSgoAISEhlCxZkrlz5xIREcHMmTPRNI1XXnmFTz75hFq1amWrXlYyn7Vs2bJmj5ctW5aEhARu3LhhktwVK1bM7Dm2trZkZGQ89L6Z7OzsaNCgATExMVy5coXdu3eTnp5O8+bN8fX1pWzZssTExDBgwABiYmLQ6XQ5ltzlRPzbt2+nUaNGODo6GsouXLgAYPRvKpO5mb3jx49n4MCBj5zJ/LCE/datW2iaZhi3KIQQ+ZHMlhUAODo60rx5c9auXUt6ejqBgYHcuXPHcNzNzQ1N07h58yYZGRlmX+np6TRq1MhwTtGiRQH4+++/Te6Xnp5uSNQe9LA11Z4kjrfffpu9e/dy9epV1q9fT58+fdi5cyetWrXiypUr2a6XVVxAlrOGM8tzc2285s2bo5QiJiaGmJgYHBwcaNCgAaBvpdu6dSupqanExsZSo0YNSpYsmWuxZNfevXtp3LixUdm+ffsA8PX1feT5//zzDxs3buTtt99+ZF0PDw80TeP69esmx27dukWxYsUkuRNC5GuS3Akjfn5+9OvXj4SEBMLDww3l9erVQynFzp07H/ta/v7+KKXYtWuXybF9+/YZWg2z40niyOTm5sZrr73GggUL6NmzJ9euXSM2NvaJ690v81m3b99ucuzMmTMkJCTg5eVlSHhzQ+bM15iYGLZt20aDBg0Ms5GbN2/OtWvXmDt3Lrdv3zaZJZuVzC7MJ/m7elxXrlzh+vXrvPTSS4ayjIwMNmzYQKNGjUxacu/evWuyw4qHhwfXr1+nWLFiJscebKF0dnamZs2a/PnnnyaxnDlzxmgmsBBC5EeS3AkTY8eOxd7enhkzZnDjxg0AhgwZgp2dHSEhIZw+fdrknMwWofu98847AEyaNImbN28a1R0zZswTxZbdOLZt22b2OomJicD/xmM9br2s9O7dG4CPP/7YqJUvPT2dUaNGoZSiT58+D73G06pZsyZubm6sXr2a3377zSiBy0xwpk6davT+UYoXLw7ol2vJLU5OTmiaRrly5Qxla9eu5dKlS0yaNMmkfuaYu8zXgQMHcHFxYd++fUblma8H1+QDeP31103GA/7+++8kJCTQpUuXnH9IIYTIQzLmTpgoV64cQUFBzJo1i+nTpzN58mR8fHxYuHAhvXv3pkaNGrRu3Rpvb2/u3bvHn3/+SWxsLB4eHhw7dsxwncaNG9O/f38WLFhAjRo16NSpE3Z2dqxdu5bixYtTrly5bK9plt04OnbsiKurK3Xr1qVSpUoopYiNjeXQoUPUrl2bFi1aZKteVurVq8fo0aOZPn06zz33HF26dMHJyYmffvqJ3377jUaNGvHuu+9m828ieztD2NjY0KRJE1avXg1glNxVrFiRKlWq8Pvvv2Nra8srr7zyWNd0cXGhbt26xMbG0qNHD7y9vbGxsaF9+/YmYwefNH4nJydatGjBsWPHqFq1KpcuXWLgwIFMnjzZ7GQcMB5zFxERQbt27Yxa/h5l4MCBREREsHTpUsP6f7Nnz6ZGjRr069fvsa8jhBBWyVLTdLNy6NAh1apVK1W0aFHl6uqqWrZsqY4cOWLpsAqcrNa5y5SYmKicnZ2Vi4uLunz5sqH8l19+UT179lSVKlVS9vb2qkSJEsrPz08FBQWpbdu2mVwnIyNDhYeHK19fX2Vvb6/Kly+vhgwZopKSkpSLi4uqWbOmoe7Dlt140OPGMW/ePNWxY0dVuXJl5eTkpNzd3ZW/v7/65JNPjJZtedx6j4pz2bJlqmHDhsrV1VU5ODio5557Tk2ePFn9+++/JnUf9bxNmjR56N+RObNnz1aapqlixYqZLE8yYMAApWmaqlu3brZiOXPmjGrbtq0qUaKE0ul0SqfTqcWLFz90Dbzsxv/nn3+qbt26qZCQENW2bVv1/fffP9Z5GRkZytvbWx08ePCx6t/vyJEj6vXXX1chISGqb9++qnPnziohISHb1xFCCGujKZUHm0Y+pvj4eBo0aEClSpUYMGAA6enpzJ07l2vXrnHgwIECt59nYXb69Gl8fHwICAjg66+/tnQ4Ip9au3YtU6dOZffu3ZYORQghrIZVdcuGhobi7OzM3r17DWN9evToQbVq1RgzZgwrV660cIQiuxITEylVqpRR9+vt27cZPnw4oO8OFeJJzZo1y/BvSQghhJ5VJXexsbG0adPGkNiBfs2zxo0bs27dOm7fvv3Ige3CuoSHhxMdHU3Tpk0pU6YMf//9NzExMVy8eJE2bdrI4HXxxK5du0aZMmXo3LmzpUMRQgirYlXJXWpqqtEippmcnJxITU3ll19+4eWXX7ZAZOJJtWzZkqNHj7Jp0yauXbuGnZ0d1apVY/jw4dLiIp6Ku7s7S5cutXQYQghhdawqufPx8WHv3r1kZGQYuvFSU1MNSxb89ddflgxPPIFmzZrl2E4IQgghhHg0q0ruBg0axMCBA+nTpw+jR48mPT2djz/+2LDDwf07JmS6cuUKGzduxNPT02yrnxBCiNxx584dzp07R6tWrZ5qx5Pbt2+b3S9YCGGer6/vw4epWXq67oM+/PBDVaRIEaVpmtI0Tb300ktq7NixStM0tXr1apP6S5cuVYC85CUvecnLQq+lS5c+1c/9uLg4iz+DvOSVn15xcXEP/Z6yqpY70K/wP2rUKI4dO4abmxs1atQw7GZgbikUT09PAJYuXcqzzz6bl6FaREhIiNG2YAVZYXnWwvKcIM9a0Bw/fpwePXoYfg4/rcLyc1yIJ5X5PfcoVpfcARQrVoz69esb3m/ZsoUKFSqY3UA8syv22Wefxd/fP89itBQ3N7dC8ZxQeJ61sDwnyLMWVDk1JKaw/BwXIrdZ/d6yy5cv59ChQzKzUgghhBDiMVhVy93OnTuZOHEirVq1wt3dnX379rFo0SJee+01hg0bZunwhBBCCCGsnlUld8888wy2trZ88sknJCcnU7lyZSZNmsSIESOyvcG8EEIIIURhZFXJXeXKldmwYYOlw7BqAQEBlg4hzxSWZy0szwnyrEIIkRekOSyfKUz/YRSWZy0szwnyrEIIkRckuRNCiHxEKcXw4WNQSlk6FCGElZLkTggh8pG4uDjmzJlNfHy8pUMRVmLRokXodDoWL15s6VDyhcLweUlyJ4QQ+Uhk5LekpX1KZOS3lg5F5JJDhw7Rq1cvKleujJOTE25ubjz//POMHj36oXusa5qWh1Far+3bt6PT6ZgwYYLZ45qmGV4FlSR3Qghh5cLCplC6tA/e3q+xfv0poB/r1p2katXWlC7tQ1jYFEuHKHLIe++9x0svvcQ333xD9erVGTZsGH379sXJyYkZM2ZQrVo1vvvuO0uHmS9klbx17NiR48eP06FDhzyOKO9IcieEEFYuNHQUoaGjSU62JTFxFaCRmLiKlBQ7QkNHExo6ytIhWkRejz/M7ftNnDiRTz75BC8vL44cOcK6deuYMmUKn376Kfv27WPlypVkZGTQvXt3tm/fnisxFCRZ/T0VLVqUatWqUbRo0TyOKO9IcieEEFbOzs6O4OA+uLkZt0S4uWkEB/fBzs7OQpFZVl6PP8zN+507d46PPvqIIkWKsGbNGrN77Hbq1Inw8HDS09MZOHCgSfKilGL9+vXUr18fFxcX3N3deeONNzhz5ozJtRITExk1ahQ+Pj64uLhQvHhxfH196dWrF3/88YdJ/Y0bN9KmTRtKliyJg4MDVatWZfTo0SQlJZnU9fT0xMvLi+TkZEaMGIGnpydFihRhwoQJBAUFodPpWLNmjdnPYf/+/eh0Orp27WooO3XqFO+//z61a9emVKlSODg44OnpyYABA7h48aLR+T179qRZs2YATJgwAZ1OZ3jt3LkTePiYu7i4ODp37kzp0qUN9xk8eDB///23Sd2ePXui0+k4f/488+fPx8/PD0dHR8qUKcOAAQO4efOmyTlHjx4lICAAT09PHBwcKF26NLVq1SIkJIS0tDSzn8mTsKp17oQQQmQtIyMNR8coihVbzo0b3cjIyLn/DPKj+8cfRkXVytf3+/LLL0lPT6dLly7UqFEjy3p9+/ZlwoQJnDx5kh07dtCkSRPDse+//56ffvqJTp060axZMw4fPsx3333Htm3b2LNnD9WqVQPg9u3bNGjQgLNnz9KyZUvat2+PUopz586xZs0a3njjDby8vAzXnTBhAhMmTKBEiRK0bduW0qVL8/PPPzNjxgx+/PFH9u7di6urq6G+pmmkpqbStGlTbty4QevWrSlatCiVK1emVatWLFiwgCVLltCuXTuT58tMuHr27Gn0XPPnz6dZs2Y0bNiQIkWK8OuvvxIVFcXatWs5dOgQ5cqVA/RdrpqmsXjxYpo0aWL0+Xh6ehrd68Fu23Xr1tG5c2c0TaNLly5UqlSJQ4cOERkZyerVq9m1a5fJNQDeffddNm3aRLt27WjdujVbt27liy++4MyZM8TExBjqHT16lJdffhkbGxvatWuHl5cXN2/e5PTp00RGRjJp0iRsbXMoLVP5XFxcnAJUXFycpUMRQohc1bXrIBUREaVu3bqlIiKiVNeugywaT079/M3OdUJDJ6tSpaqpqlVbKw+PDgoylIdHB1WlSitVqlQ1FRo6+alisdT9mjVrpjRNU1FRUY+s+9ZbbylN09SkSZOUUkp9+eWXStM0pWmaWr9+vVHdWbNmKU3TVPPmzQ1la9asUZqmqREjRphc+969eyo5OdnwfuvWrUrTNNWgQQOVlJRkVHfRokVK0zQVEhJiVF6pUiWlaZp69dVX1e3bt03u4ePjo+zt7dW1a9eMyu/evauKFy+uypQpo9LT0w3lFy9eVKmpqSbX2bRpk7KxsVEDBw40Kt+2bZvSNE1NmDDB5Byl/vd5LV682FCWnJys3N3dla2trdq1a5dR/WnTpilN01TLli2NygMDA5WmaapSpUrqwoULhvK0tDTVuHFjpWmaOnDggKF8xIgRStM0tWbNGpOYbty4oTIyMszGe7/H/V6RblkhhMgnli+fQ3C/t3D69VeCg/uwfPkcS4eU5/J6/GFe3e/SpUsAVKhQ4ZF1n3nmGQCTmbPNmzenTZs2RmVDhgyhcuXKbN26lT///NPomIODg8m1bW1tcXFxMbyPiIgA4IsvvjAZoxYYGMgLL7zA119/bXIdTdP49NNPcXR0NDkWGBhIamoq0dHRRuVr167lxo0bvPXWW0ZbjpYrV87s0INXX32V6tWrs3HjRpNj2bV69WquX79Ot27daNCggdGxkSNHUqlSJTZv3syFCxdMzg0LCzP8nQDY2NjQq1cvAA4ePGhS39zn7ubmlqOzdyW5E0KI/GLTJvDzgzZt4M4dS0djEXk9/jA/jXd85ZVXTMp0Oh0NGzYE4MiRIwA0adKE8uXLM3XqVF577TUiIiKIj48nIyPD5Py9e/diZ2fHihUrGD9+vMkrNTWVf/75h+vXrxud5+DggJ+fn9k433nnHbNj3sx1yWZaunQpLVq0oFSpUtjZ2RnG0f36668PXR7mcWWOo8wcr3c/GxsbGjduDMDhw4dNjteuXdukLDPZu/9z6d69OzY2NnTo0IHAwECWLFnC77///tSxmyNj7oQQwtolJEBICKxcCU2awOrVYKZFpDDJ6/GHuX2/MmXKcOLECZPWNXMyW48yx5ll8vDwyPLagGHyg6urK/v27WPcuHGsWbPG0PJVsmRJBg0axNixYw1jv65evUp6enqWa8aBvpUuJSWF4sWLG8pKly6dZf3y5cvTvHlzNm/ezIkTJ/D19eXy5cts2LCBmjVr8txzzxnVDwkJYdasWZQrV47XXnuN8uXLG1oEv/zyy8f6zB4l87MpW7as2eOZ5eYmkBQrVsykLPPzS09PN5TVqVOH2NhYJk2axMqVK/nqq68A8PHxYdy4cXTv3v3pHuI+0nInhBDW6t49mDEDfH0hNhaWLoWtW6F6dUtHZnH+/l5Mm6Zx5sxqpk3T8Pf3evRJVny/Ro0aAbBly5aH1ktPTzcsg/Jg92FiYqLZczJnerq5uRnKypcvT1RUFJcvX+bXX38lIiKCEiVKMHHiRCZOnGio5+bmhru7OxkZGVm+0tPTTbqTH9XFGBgYCPyvte7rr78mPT3dUJ7p8uXLRERE4Ofnx8mTJ1myZAlTpkwhLCyMsLAwihQp8tD7PK7Mz8bcrFj4X7f5/Z/hk6hbt66h+3n37t2EhoaSmJjIm2++aTT54mlJcieEENYoNhZq1oT33oPeveHkSXjrLSjAq+pnx/LlcwgO7oOTk1OejD/M7fv17NkTGxsbVq1axbFjx7Kst3DhQi5duoSvr69JN6y5te/S09PZtWsXmqZRs2ZNs9esXr06Q4YMYfPmzYB+/FmmevXqce3atYfG9CQ6deqEq6srX3/9NUopFi9ejJ2dHW+++aZRvbNnz6KUomXLljg7OxsdS0hI4OzZsybXtrGxAYxbzR7F398fgG3btpkcS0tLIzY2Fk3TDPWelp2dHfXq1WPChAmGcY1ZLQ/zJCS5E0IIa3L5MgQGQuPG4OIChw5BRAQ8ZYuBsG5eXl6MGTOGe/fu0a5dO44fP25S54cffmDYsGHY2toSGRlpcnzr1q2sX7/eqOzzzz/n7NmzNG3a1NC6duzYMbOtfJmtVk5OToaykJAQAPr162dovbrfrVu32L9/fzaeVM/BwYFu3bqRkJDAZ599xtGjRw3r6N0vc0mW2NhYozGBKSkp9OvXz2wCV6JECQDOnz//2PF06NABd3d3oqOjTZ5n5syZnDt3jhYtWhhNnMiuPXv2cPfuXZNyc5/705Ixd0IIYQ3S02H+fPjwQ9DpYMEC6NNH/7UoFMaPH8+tW7f47LPPeOGFF2jVqhXVq1fn3r177NmzhwMHDuDk5ER0dLTZyRNt27alY8eOdOzYkSpVqnDkyBE2bNhAiRIlmDt3rqHepk2bePfdd6lfvz7e3t6ULl2ahIQEVq9ejY2NDe+++66hbrNmzZg6dSoffPAB3t7etGnTBk9PT1JSUjh//jw7d+6kUaNG/Pjjj9l+3sDAQKKiohgzZozh/YM8PDzo3r07y5Yt48UXX+TVV18lKSmJzZs34+TkxIsvvmiYKJLJ19eX8uXLs2zZMuzs7KhYsSKapvHOO+9QsWJFs7E4OzuzcOFC3njjDV555RXeeOMNKlSoQFxcHJs3b6Zs2bLMnz8/2894v+nTp7Nt2zYaNWqEp6cnLi4u/Pbbb2zYsAF3d3f69+//VNc38shFVSzg1KlTqlu3buqZZ55RTk5OytfXV02cONHsejmyzp0QIt87cECpWrWUAqX69lXqn38sHdFjscQ6d4XBgQMHVGBgoPLy8lKOjo7K1dVV+fn5qXfffVddvHjRpP6iRYuUTqdTixcvVuvWrVP16tVTzs7Oqnjx4qpLly7q9OnTRvWPHz+uRowYoWrXrq1KlSql7O3tlZeXl3rjjTfU3r17zca0a9cu1bVrV1WuXDlVpEgRVbp0aVWzZk01cuRIk783T09P5eXl9VjP6u3trXQ6nSpZsqS6d++e2Tq3b99WH374oapatapycHBQFStWVEOGDFFXr15VTZo0UTqdzuScgwcPqubNmys3Nzel0+mUTqdTO3bsUErp17nL/LzMndexY0dVqlQpVaRIEVWpUiU1aNAgdenSJZO6PXv2VDqdTp0/f97kmLm19jZt2qR69eqlqlevrtzc3JSzs7Py9fVVw4YNU3/++edjfV6P+72iKZVHm/I9pgsXLvD8889TvHhxgoKCcHd3Z8+ePSxatIh27drxww8/GNWPj4+nVq1axMXF5VhfuBBC5Inr12HMGH2L3QsvwNy5UK+epaN6bDn181d+jgvxeB73e8XqumW/+uorkpKS2LNnj2Fvvb59+5KRkcGSJUtISkp66tkqQghhUUrB4sUwejTcvQvh4TB4MOTU1kNCiELN6gZzZG60++AaOWXKlMHGxibHpj0LIYRF/PKLfrJEr17w6qv6WbDDhkliJ4TIMVaX3DVt2hSAPn368PPPP3PhwgWWL1/OvHnzGDp0qNmtTIQQwuolJ8PIkfrlTa5c0a9X9/XXkMWiqUII8aSs7lfFVq1m0yy0AAAgAElEQVRa8dFHHzF58mSjNV/Gjh1rtLCiEELkC0rBt9/qd5i4cQMmTdJ/Lb0QQohcYnXJHUClSpV45ZVX6Ny5MyVKlGDdunVMmjQJDw8PBg8ebOnwhBDi8Zw6BUOGwObN0LEjzJwJWSzFIIQQOcXqkrtly5YxYMAATp8+bdg3r0OHDmRkZPDee+8REBCAu7u7yXkhISEmEy0CAgIICAjIk7iFEMLg9m2YMgWmT4fy5WHdOnj9dUtH9VSio6OJjo42KjO3z6YQwvKsLrmbO3cu/v7+Jhsit23blkWLFnHkyBGaNWtmcl54eLhMoRdCWN66dTB0KFy8CO+/r38VgLHC5n5ZzlyWQQhhXaxuQkViYqLZ7UTu3bsH6Pd4E0IIq3PuHLRvD23bgrc3/PorTJhQIBI7IUT+YnXJnY+PD/Hx8Zw+fdqoPDo6GhsbG55//nkLRSaEEGb8+y9MngzVq0NcnH7yxIYN+gRPCCEswOq6Zd99911++uknGjVqxJAhQ3B3d2fdunVs2LCBfv36UaZMGUuHKIQQejEx+sWHz5yB4cNh3DhwdbV0VEKIQs7qWu4aNWrEnj17qFWrFnPnziUkJIQ//viDyZMnExkZaenwhBAC/voLAgKgRQsoVQoOH4YZMySxE0JYBatruQOoU6cO69evt3QYQghhLC0NPv8cwsLAwUG/hdjbb4OmWToyIYQwsLqWOyGEsEp79kCtWjBihD6hO3kS3nlHEjshhNWR5E4IIR7myhXo0wcaNAB7ezhwAObMgeLFLR2ZEALQ6XSGrUuFniR3QghhTkYGLFgAPj6wahVERsLevVC7tqUjEwWQTqczetna2lKiRAmaNm3K4sWLLR2e1dOkBd2IVY65E0IIi4qPh0GDYP9+6NULpk3TT5wQIhdpmsa4ceMA/dqup0+fZtWqVezYsYNDhw4xe/ZsC0donU6cOIGTk5Olw7AqktwJIUSmpCQYOxbmzoUaNSA2Fho2tHRUohAJCwszer9nzx4aN27M3LlzGTlyJJ6enpYJzIpVq1bN0iFYHemWFUIIpWDpUn0X7KJF+mVN4uMlsRMWV79+fXx8fFBKER8fb3J8//79dOnShTJlymBvb0/FihUJCgri0qVLZq938OBBWrZsiaurK25ubrz66qvs27eP8ePHo9Pp2Llzp1H9zPFsiYmJ9O3bl/Lly2Nra2vUVZydGM6ePUv//v2pWrUqTk5OlChRgueff56BAwdy7do1Q73U1FQiIiLw9/fH3d0dZ2dnvLy86NChAzExMWZjfFBSUhIffPABPj4+ODo64u7uTuvWrU3OB9i+fTs6nY4JEyZw5MgRXn/9dYoVK4azszNNmjRh7969Zj9PayUtd0KIwu3YMX0X7I4d0K0bfPoplC9v6aiEMGFnZ2f0fuHChfTv3x9HR0fatWtHhQoVOHXqFFFRUaxdu5Z9+/ZRoUIFQ/2dO3fSsmVLlFJ06tSJKlWqcPToUZo2bWp2z/ZM165do27duri6utKlSxd0Op1hQ4HsxHDp0iXq1KlDcnIyr7/+Om+88QZ3797l7NmzLF26lODgYNzd3QHo2bMny5Ytw8/Pj8DAQBwdHbl48SK7d+9m48aNNG/e3CjGB8fc3bhxgwYNGnD8+HFeeuklOnfuzD///MOKFSto2bIlkZGR9O/f3+RZDx06xPTp06lfvz79+/fn/PnzfPfddzRv3pwjR47kn1ZClc/FxcUpQMXFxVk6FCFEfpKcrNTo0UrZ2ipVrZpSmzdbOqJ8J6d+/j7xdW7dUiouLu9ft2491fOao2ma0ul0JuU7duxQOp1OOTg4qL///ttQfvLkSWVnZ6e8vb3VX3/9ZXROTEyMsrGxUR07djSUpaenq6pVqyqdTqc2bNhgVH/evHmG++/YscMkLk3TVGBgoEpPTzc6lt0YIiIilKZpKiIiwuQ5b9++re7cuaOUUurGjRtK0zRVp04dlZGRYVL36tWrJjE2bdrUqKx///5K0zQVFBRkVH769Gnl5uam7O3t1blz5wzl27ZtMzzr4sWLjc6ZP3++0jRNDRo0yCSWvPa43yvScieEKFyU0s9+HTZMv8zJ+PEwapR+mRORv5w4oV97MK/FxYG/f45fVinFhAkTUEpx7949zpw5w6pVq9A0jRkzZuDh4WGoGxkZSVpaGrNmzaJs2bJG12nWrBlt27Zl7dq13Lp1C2dnZ/bs2cPvv/9Os2bNaNWqlVH9/v37Ex4ezqlTp8zGZW9vz4wZM9DpjEdyZTeGTA4ODib3cHR0NHyd2Qpnb29vdhZsZuteVlJTU1m6dCmurq5MmTLF6FjVqlUZOnQoH3/8MUuWLCE0NNToeMOGDXnnnXeMynr37s3gwYM5ePDgQ+9rTSS5E0IUHmfOQHAwbNgA//d/EBEBXl6Wjko8KV9ffaJlifvmkgkTJhi91zSNhQsXEhgYaFSeOQZs+/bt7N+/3+Q6ly9fJj09nVOnTlGzZk0OHz4M6JOXB2maRr169bJM7jw9PSlZsqRJ+ePGcPLkSfz9/Wnfvj0ffvghgwcPZuPGjbRs2ZKGDRtSvXp1o/OKFi1qSAxffPFFOnfuTKNGjXjppZcea1bsyZMnuXPnDg0bNqRYsWImx5s1a8bHH3/MkSNHTI7VNrPUka2tLR4eHly/fv2R97YWktwJIQq+O3f0y5lMnQplysDq1dCunaWjEk/LySlXWtAsRdM00tPTAbhz5w579uyhT58+DBgwgIoVKxpNGrh69SoAn3zyyUOvl5KSAugnFwBGrX/3y6ocMIyve9DjxnDr1i0AKlasyIEDBxg/fjwbNmzg+++/B6BChQqMGjWK4OBgw3nLly9n2rRpfPPNN4blYRwcHOjSpQszZsygdOnSWd4z81kfbE188Hlu3LhhcsxcMgj6BC/z7yY/kNmyQoiC7aef4LnnYPJkGDlSP4FCEjth5RwdHWnevDlr164lPT2dwMBA7ty5Yzju5uaGpmncvHmTjIwMs6/09HQaNWoE6FvDABITE83eL6tyyHqB4OzGAODr68uyZcu4evUqhw4dYurUqWRkZDBs2DAWLlxoqOfg4MC4ceM4efIkf/75J0uXLqVhw4YsXbqULl26PPSzc3NzA+Dvv/82ezxzFm9mvYJIkjshRI5SSjF8+BiUUpYN5MIF6NwZ2rTRd70ePQqTJulbe4TIJ/z8/OjXrx8JCQmEh4cbyuvVq4dSymTpkqz4/7eFMzY21uRYRkYGe/bsyXZs2Y3hfjY2Nvj7+zN69Giio6MBWL16tdm6zzzzDG+++SYbN26kSpUq7Nq166FdpL6+vjg6OvLzzz8bWvHut23bNuB/n0lBJMmdECJHxcXFMWfObLNrcuWJ1FSYPl0/LmrvXoiOhs2bc3WclBC5aezYsYZJDZldiUOGDMHOzo6QkBBOnz5tck5qaqpRItegQQOqVKnCtm3b2LBhg1HdBQsWcPr06Wxv4ZXdGOLj480mW5ktbJnj6a5cucIvv/xiUi8lJYWUlBTs7OwoUqRIlnHZ2dnRo0cPbt68aTJh4vfffyciIoIiRYrw9ttvP96D5kMy5k4IkaMiI78lLe1TIiO/JSoqj2cybt+uX7Pu1Cn9xIkJE+C/3VFC5FflypUjKCiIWbNmMX36dCZPnoyPjw8LFy6kd+/e1KhRg9atW+Pt7c29e/f4888/iY2NxcPDg2PHjgH6rtWoqChat25Nu3bt6Ny5M5UrV+bo0aNs2bKF1157jZ9++slkRuzDZDeGJUuWsGDBAho2bEjlypUpXrw4v//+O2vXrsXBwYHhw4cDkJCQgL+/P35+fvj5+VGhQgVu3rzJunXrSExMZNiwYUazb82ZOnUqsbGxfP755xw8eJAmTZpw5coVVqxYwa1bt/j888+pVKnSE/6N5AO5vihLLpN17oSwvNDQyapUqWqqatXWysOjg4IM5eHRQVWp0kqVKlVNhYZOzt0ALl1SqkcPpUCp+vWVOnIkd+8nlFJWsM5dAZLVOneZEhMTlbOzs3JxcVGXL182lP/yyy+qZ8+eqlKlSsre3l6VKFFC+fn5qaCgILVt2zaT6+zfv1+9+uqrytXVVbm6uqpXX31V7du3Tw0ePFhpmqZ+/vlnk7geXEPuQY8bw/79+9XAgQPVCy+8oNzd3ZWjo6Py9vZWvXv3Vr/99puh3o0bN9TEiRNVs2bNVPny5ZW9vb0qV66catq0qVq2bJnZz85cjDdu3FDvvfee8vb2Vvb29qp48eKqZcuWarOZNS0z17mbMGGC2Wf09PRUXl5eD/0c8sLjfq9YXXIXGBhoWEjQ3OvBhRLlh4IQlpeamqoiIqKUh8f/Kf1CcvqXh8f/qYiIKJWampo7N05LU2r2bKWKFlWqRAml/vMfpR5YaFXkHknuCo769esrOzs7dfv2bUuHIh4i3y5iHBQURMuWLY3KMjIyCAoKwsvLK8upzUIIy7GzsyM4uA+ff76a+yfdublpBAf3yZ2b7tun74I9cgT699fPhn3E4qZCFGZ37tzh33//NVnuY9GiRezdu5c2bdoYLSYs8i+rS+7q1q1L3bp1jcp27drF7du3eeuttywUlRDicWRkpOHoGEWxYsu5caMbGRlpOX+Tq1fhgw/giy+gZk39pImXX875+whRwJw/f56aNWvSsmVLqlSpQlpaGocPH2b37t0UL16cTz/91NIhihySL2bLfvPNN2iaxptvvmnpUIQQD+Hv78W0aRpnzqxm2jQNf/8c3P0hIwP+8x/w8YEVK+Dzz+HgQUnshHhMZcqUoUePHhw/fpyoqCjmz5/PhQsX6N27N4cOHcLHx8fSIYocYnUtdw+6d+8eK1asoEGDBlSsWNHS4QghHmL58jmGr4OD++Rcl+yRI/ou2L17oUcP+OQT/U4TQojHVqxYMb744gtLhyHygNW33G3cuJFr165Jl6wQhdHNmzB8uH5z+KQk/VInX30liZ0QQjyE1bfcffPNNxQpUoSuXbtaOhQhRF5RCpYt028XlpSk3xN2+HCws7N0ZEIIYfWsOrlLSUlh9erVtGrViuLFiz+0bkhIiMk+cQEBAQQEBORmiEKInHbiBAweDFu3QqdOMHMmVKhg6agKvejoaMM2UZnM7TYghLA8q07ufvjhB+7cufNYXbLh4eEFep84IQq8W7f0e7/OmAEVK8JPP0Hr1paOSvyXuV+W4+PjqVUrj3chEUI8klUnd19//TWurq60a9fO0qEIIXKLUrBmDQwdComJ8OGH8N574OBg6ciEECJfstrk7p9//mHLli289dZbOMgPeSEKpj/+QA0dirZuHap1a7StW6FKFUtHJSzk+PHjlg5BCKv2uN8jVpvcLV++nPT0dJklK0RB9O+/MH06TJ7MvaJFeUvnwPsffUQtSewKtR49elg6BCEKBKtN7r755hs8PDxo0aKFpUMRQuSkTZtgyBD44w8YMYLhl1JZ+dWzuM1bSVRUbUtHJyzA19eXuLg4S4chRL7h6+v70ONWm9zt2bPH0iEIIXLSxYswYgSsWMEflbx427U8id8fJTnZAfiMdes6UbVqa27e/IOgoJ5MnPiBpSMWecTJyUkmxAmRg6x+EWMhRD537x589hn4+sKOHbB0Kc+cOkG3CaEkJ9uSmLgK0EhMXEVKih2hoaMJDR1l6aiFECLfkuROCJF7du0Cf394913o2VO/ht1bb2FXpAjBwX1wc9OMqru5aQQH98FOFisWQognJsmdECLnXb4MvXpBo0bg5AQHD8Ls2VCsmFG1jIw0HB2jKFv2VRwdo8jISLNQwEIIUXBIcieEyDnp6RAZCT4+sHo1LFgAe/fqW+/M8Pf3Yto0jTNnVjNtmoa/v1ceByyEEAWP1U6oEELkM4cOwcCB+j9799bvB1uq1ENPWb58juHr4OA+BAf3ye0ohRCiwJOWOyHE07l+HQYNgpde0k+e2L0b/vOfRyZ2Qgghcoe03AkhnoxSsGSJfrLE3bsQHg6DB4Ot/FgRQghLkpY7IUT2/forvPKKfgZsixb6WbDDhkliJ4QQVkCSOyHE40tOhlGj4MUX9TNit2yBb76BcuUsHZkQQoj/kl+zhRCPphSsXAkhIXDtGnz0kX63CXt7S0cmhBDiAdJyJ4R4uNOnoXVr6NoVatWCY8fggw8ksRNCCCslyZ0Qwrw7dyAsDJ57Dk6dgjVr9GvXeXpaOjIhhBAPId2yQghT69bB0KFw8SKMHq1vqXNysnRUQgghHoMkd0KI/zl/Xj/rdfVqePVV2LABqlWzdFRCCCGyQbplhRCQmqrfUeLZZ/X7wK5YARs3SmInhBD5kFUmd/Hx8bRr144SJUrg7OyMn58fs2fPtnRYQhRMW7fCCy/A2LH67cNOnIA33gBNs3RkQgghnoDVdctu2rSJtm3bUqtWLcLCwnBxceHMmTNcvHjR0qEJUbBcuqRfs+6bb6BBA31rnZ+fpaMSQgjxlKwqubt58ybvvPMObdu2ZeXKlZYOR4iCKS0N5s6F0FD9ciZffgnvvAM6q2zIF0IIkU1W9dP8m2++4fLly0yaNAmAW7dukZGRYeGohChA9u6F2rVh+HAICNB3wfbsKYmdEEIUIFb1E33Lli0ULVqUCxcu4OPjg6urK25ubgwaNIh///3X0uEJkX9duQJ9+0L9+vr9X/fvh3nzwN3d0pEJIYTIYVaV3J0+fZq0tDQ6dOjAa6+9xvfff0/v3r2ZN28evXr1snR4QuQ/GRkQFQU+PvDdd/ru2P37oU4dS0cmhBAil1jVmLuUlBRu377NwIEDmTlzJgAdOnQgNTWV+fPnM3HiRKpWrWrhKIXIJw4fhkGDYN8+CAyE6dOhdGlLRyWEECKXWVXLnaOjIwABAQFG5Znv9+3bl+cxCZHvJCXpd5eoXRuSk2HHDli0SBI7IYQoJKyq5a5cuXIcO3YMDw8Po/LS//1P6fr161meGxISgpubm1FZQECASaIoRIGllH5Zk5EjISVF31I3dCjY2Vk6MlEAREdHEx0dbVSWlJRkoWiEEA9jVcld7dq12bJlCwkJCXh7exvK//rrLwBKlSqV5bnh4eH4+/vneoxCWKXjx/VdsNu36xcg/uwzeOYZS0clChBzvyzHx8dTq1YtC0UkhMiKVXXLdu3aFYD//Oc/RuVRUVHY2dnRpEkTC0QlhBW7dQvefx+efx4SEvR7wa5YIYmdEEIUYlbVcvfiiy/Su3dvFi5cSFpaGo0bN2b79u2sXLmSMWPGUKZMGUuHKIR1UAp++AGGDYN//oGwMHj3XXBwsHRkQgghLMyqkjuAefPmUbFiRb788ktWrVqFp6cnM2fOZOjQoZYOTQjrcPYsBAfDjz9CmzYwezZUrmzpqIQQQlgJq0vubG1tCQsLIywszNKhCGFd7t7VT5KYMgVKlYJVq6B9e9A0S0cmhBDCilhdcieEMGPjRhgyBM6d08+GDQ0FZ2dLRyWEEMIKWdWECiHEAxIS9LNfW7eGChXg6FGYOlUSOyGEEFmS5E4Ia3TvHsyYAb6+sGsXfP01xMTAs89aOjIhhBBWTrplhbA2O3fq16w7flzfFTtxIjywQLcQQgiRFWm5E8JaJCbCO+/AK6+AqyscOgSzZkliJ4QQIlskuRPC0tLTYe5c8PHRL28SFQW7d0PNmpaOTAghRD4kyZ0QlnTgALz8MgweDF27wsmT0KcP6ORbUwghxJOR/0GEsIRr1yAoCOrW1bfc7d0LCxZAiRKWjkwIIUQ+JxMqhMhLGRmweDGMHg2pqRARAQMHgo2NpSMTQghRQEjLnRB55ehRaNwYevfWr1t38qR+NqwkdkIIIXKQJHdC5LbkZBgxAvz99d2x27bBV19BmTKWjkwIIUQBJN2yQuQWpWDFCn1id+MGTJ4Mw4dDkSKWjkwIIUQBJi13QuSGkyehZUvo3l0/aeL4cf04O0nshBBC5DJJ7oTISbdvw9ix4OcHZ8/q16377juoWNHSkQkhhCgkpFtWiJyydi0EB8OlSzBmDLz3Hjg6WjoqIYQQhYwkd0I8rT/+gGHD9Mldq1awZQtUrWrpqIQQQhRSVtctu337dnQ6ndnXgQMHLB2eEP/z778waRJUrw7x8bByJfz0kyR2QgghLMpqW+6GDRtGnTp1jMqqVKlioWiEeMCWLfotw86ehZAQCAsDFxdLRyWEEEJYb3LXqFEjOnXqZOkwhDB28SKMHAnLl+sXJP7uO3juOUtHJYQQQhhkK7nr3LkzXl5e+Pj48Oyzz9KwYcPcigulFMnJyTg6OmJra7U5qCgs0tJg9mx9C52TEyxZAj16gKZZOjIhhBDCSLbG3MXExNCuXTteeeUVlFK5FRMAvXr1ws3NDUdHR5o1a0ZcXFyu3k+ILO3eDbVqwahREBioX8Pu7bdzNbFTSjF8+Jhc/z4TQghR8GQruWvatCmNGzemWrVqNGrUKFcCsre3p0uXLkRERLBmzRo+/vhjfvnlFxo1asSRI0dy5Z5CmPXPP9CrFzRsCPb2cOAAfP45FCuW67eOi4tjzpzZxMfH5/q9hBBCFCzZ6u90yYMB4/Xq1aNevXqG9//3f/9Hly5deP755/nggw/46aefcj0GUcilp0NUFHzwgf79/PnQty/o8m5yeWTkt6SlfUpk5LdERdXKs/sKIYTI/7KV3P3xxx9cuXKFkiVLPrJueno6NjY2TxzY/apUqUL79u35/vvvUUqhmekOCwkJwc3NzagsICCAgICAHIlBFBJxcTBwIBw8qG+1mzYNSpXKk1uHhU1h3rxFuLlVJjnZAZjKunWdqFq1NTdv/kFQUE8mTvwgT2IR4kHR0dFER0cblSUlJVkoGiHEw2QruduzZw+lS5fm2WefpVGjRoZXRTNbK/Xu3ZvFixfnWKDPPPMMqamp3Lp1y2wLYnh4OP7+/jl2P1HIXL+u3zYsMlK/ddiuXdCgQZ6GEBo6ilKlSjNp0g8kJq4C+O+fbQkNHU1Q0Dt5Go8Q9zP3y3J8fDy1aknLshDWJlv9TI0bN2bBggXUrl2bTZs28fbbb+Pp6UmlSpXo0aMH8+fP59ixYwAkJyfnaKBnz57F0dExT7qGRSGilH7mq6+v/s9PP9W33uVxYgdgZ2dHcHAf3NyMW6bd3DSCg/tgZ2eX5zEJIYTIf7LVclehQgX69u1L3759Abh48SKxsbHs2rWL2NhYli1bRkZGBsWLF+fOnTtPFNA///xDqQe6wX7++WfWrFnD66+//kTXFMKsX3/VL0S8cyd0765P7MqVs3RUZGSk4egYRbFiy7lxoxsZGWmWDkkIIUQ+kq3k7tq1a0bvy5cvT/fu3enevTsAN27cYPfu3ezYsYM5c+Y8UUDdunXDycmJevXqUbp0aY4dO8aCBQtwcXFh6tSpT3RNIYykpMCECTBzJlSuDJs3Q4sWlo7KwN/fi6FDNfr0Wc1//hPNrl1elg5JCCFEPqKpbCyk5eLiwu+//46Hh8cj63bo0IEffvgh2wHNnj2br7/+mjNnznDz5k1Kly5N8+bNGTduHJUrVzapnznmIy4uTsbcCSNKKUJCPiQ8fJJ+Eo5S+h0lQkLgyhUIDdXvNmFvb+lQhciX5OevENYpW2PuBg8eTI8ePVi6dCl37959aF1nZ+cnCig4OJh9+/Zx5coVUlNTSUhIYPHixWYTOyEexmituNOn4bXX4I03oGZNOHYMxoyRxE4IIUSBk63kbtq0aWzevJlnn32WuXPnPrSudKEKS4uM/BbbtClc7DdUv//riROwejWsWQNe0tUphBCiYHqiVVlr1arFiBEjHlqnQoUKTxSQEE8jLGwKpUv74O39Gne/38mvhNP68F7mOD1DpVt2hB36zdIhCiGEELkqWxMqhLB2oaGj8LKxwWPKZ7T5N5EtNKcNP5JkP4oPP+wga8UJIYQo8CS5EwVHaip24eH0mj6Ry2npdGMZK+gKaFT771pxQgghREEnyZ0oGLZt069Zd+oUDB1Kq9W/cPJSMmWLtZS14oQQQhQqebcTuhC54e+/oUcPaNYMiheH+Hj47DOq1a7GtGkaZ86sZto0DX9/mUAhhBCicJCWO5E/paXp94EdOxaKFIGFCyEwEHT631eWL//fItrBwX2kS1YIIUShIS13Iv/Ztw/q1IFhwyAgAE6ehF69DImdEEIIUZjJ/4Yi/7h6Ffr3h3r19Incvn0wbx64u1s6MiGEEMJqSLessH4ZGfDll/Dee/ru2DlzYMAAsLGxdGRCCCGE1ZGWO2HdjhyBhg2hb194/XV9F+ygQZLYCSGEEFmQ5E5Yp6Qk/Zi6WrX0X2/fDosXg4eHpSMTQgghrJp0ywrrohQsWwYjRkByMkybpk/y7OwsHZkQQgiRL0jLnbAex49D8+bw5pvQoIH+/ahRktgJIYQQ2SDJnbC8W7fggw/ghRfgwgXYsAFWroQKFSwdmRBCCJHvSLessBylYPVqfbdrYqJ+QeLRo8HBwdKRCSGEEPmW1bfcTZo0CZ1Oh5+fn6VDETnp7Flo2xY6doQaNeC33yAsTBI7IYQQ4ilZdXKXkJDA5MmTcXZ2RtM0S4cjcsK//8JHH+kTuqNH4fvvYf16qFLF0pEJIYQQBYJVd8uOGjWK+vXrk5aWxpUrVywdjnhamzbBkCHwxx8wciSEhoKzs6WjEkIIIQoUq22527lzJ9999x0zZ85EKSUtd/lZQgJ07QqtWkH58vDzzzB1qiR2QgghRC6wyuQuPT2d4OBg+vXrR40aNSwdjnhS9+7Bp5+Cry/s3AlLl8LWrVC9uqUjE0IIIQosq+yWnTdvHn/++Sdbt261dCjiScXG6rcJO3YMBg/Wj7Nzc7N0VEIIIUSBZ3Utd1evXiUsLIywsDBKlChh6XBEdl2+DD17QuPG+m7XQ4cgIkISOyGEECKPWF3L3dixYylZsiTBwcHZOi8kJAS3BxKIgIAAAgICcjI8kZX0dFiwAMaMAZ1O/3WfPvqvhRD5XnR0NNv7IfsAABgVSURBVNHR0UZlSUlJFopGCPEwVpXcnT59mi+++IKZM2eSkJBgKL979y6pqamcP3+eokWLUrx4cZNzw8PD8ff3z8twRaaDB/VdsIcO6RO6qVOhZElLRyWEyEHmflmOj4+nVq1aFopICJEVq2pWuXjxIhkZGQwdOpTKlSsbXgcOHODUqVN4eXnx0UcfWTpMken6dX1S9/LL+skTu3dDVJQkdkIIIYQFWVXLnZ+fH6tWrTJa9kQpxdixY0lJSWHWrFlUkcVuLU8pWLIE3n0X7t6FmTP1SZ5t9v85KaUICfmQ8PBJstyNEEIIkQOsKrkrUaIE7du3NykPDw8HoF27dnkdknjQL7/oE7lduyAgQL/USdmyT3y5uLg45syZzdtvd5buHSGEECIHWFW3bFY0TZNWHUtLTtbvKlGzJly5AjEx8M03T5XYAURGfkta2qdERn6bQ4EKIYQQhVu+SO62bdvG0aNHLR1G4aQUrFihX4g4MhI+/li/w0SzZk98ybCwKZQu7YO392usX38K6Me6dSepWrU1pUv7EBY2JefiF0IIIQqZfJHcCQs5dUq/ZVi3bvDSS3D8OLz/PhQp8lSXDQ0dRWjoaJKTbUlMXAVoJCauIiXFjtDQ0YSGjsqZ+IUQQohCSJI7YerOHQgNBT8/OH0a1q6FVaugUqUcubydnR3BwX1wczPuandz0wgO7oOdnV2O3EcIIYQojKxqQoWwAuvWwdChcPEivPcefPABODrmyq0yMtJwdIyiWLHl3LjRjYyMtFy5jxBCCFGYSMud0Dt/Hjp0gLZtwdtbPyt24sRcS+wA/P29mDZN48yZ1UybpuHv75Vr9xJCCCEKC0nuCjClFMOHj0EplXWl1FSYMgWefVa/08SKFbBhA1SrluvxLV8+h+DgPjg5OREc3Ifly+fk+j2FEEKIgk6SuwIscw25+Ph48xW2boUXXtCPrxs0CE78f3v3HlVVmbhx/DkEKogy3kjRlZgXDMLliIzilLp0vKCBLRWRMrW8TKVYaJelK3HKsTQtHWU0b2PmhTIcEyidxqVWU94WlD/TEDWVQlFxEkFAbvv3xxmZmKy0Duxz9vl+1mIlL+ccn52u43P2ft93Z0kxMRLbzgAA4LIodxb2o3vInT8vPfSQ1L+//VZhn38uLVokNWpkTlAAAOAwLKiwmMTEV/TGG2/Kz+9uFRY2kDRf6enD1aHDYF0r+Fpru3XSkH0fSw0aSG++KY0dy5k6AAAshDN3FvNje8h1/neh/q9+qSL/+YE0Zox0/Lg0bhzFDgAAi6HcWcz/7iHXTPlarYlK/+4ztWh1p2wHD0rLl0tNmpicFAAA1AbKnUUZleV60mussm0BGqFkJbYIkfbvl7p3NzsaAACoRcy5s6LMTL3/3efqWH5BFQ8/rDeDw3T88EnpjjvMTgYAAGoZ5c5Krlyxb2uyfLk6hoRI21Pked99mihpotnZAABAnaDcWYFhSJs3SzNmSNeu2bc1mTpV4h6tAAC4Hebcubpjx6R+/ewrYPv0sW9EnJBAsQMAwE1R7lxVUZH0/PP2O0zk5kr/+If0zjtS69ZmJwMAACZyunJ39OhRxcTEqH379mrYsKFatGihPn36KD093exozsEwpL//XQoOlpYulebMkY4ckQYONDsZAABwAk5X7nJyclRUVKTx48dr6dKlSkxMlCRFR0dr9erVJqcz2alT0tCh0ogR9jN2x45JL7wg1a9vdjIAAOAknG5BRWRkpCIjI2uMTZkyRWFhYXr99dc1adIkk5KZqLRUmj/f/tWypbR9uxQdbXYqAADghJzuzN3NeHh4qE2bNiooKDA7St3bsUO6917p5Zftq2GPHaPYAQCAH+V0Z+5uKC4uVnFxsQoKCpSamqqdO3dq9OjRZseqO998Iz39tH1+Xf/+Unq61Lmz2akAAICTc9pyN336dK1atUqS/czdiBEjlJSUZHKqOlBWJi1ZIr34ouTnJyUnS7Gxks1mdjIAAOACnLbcJSQkaNSoUcrNzdWWLVtUUVGh69evmx2rdn30kfTkk/a96uLjpZdekho3NjsVAABwIU5b7oKCghQUFCRJeuSRRzRo0CBFRUXpwIEDN318QkKC/Pz8aozFxcUpLi6u1rP+anl50rPPShs3Sr16SZmZ9tWwAOAkkpOTlZycXGPMLedBAy7AZhiGYXaIW7Fq1So9/vjjOn78uDp27Fg9npmZqbCwMGVkZKhbt24mJvwFKiulFSvs25l4ekqvviqNHy95uMQ6FwBuzqXffwELc5kWUVJSIslCnxQPHJDCw6Vp0+xz6o4flx57jGIHAAB+FadrEpcuXfrBWHl5ud566y35+PgoODjYhFQOdPmy9Mc/ShER9u/37ZNWrpSaNTM3FwAAsASnm3M3efJkFRYWqnfv3goICFBeXp42bdqk7Oxsvfbaa/Lx8TE74i9TVSWtW2e/H2xFhbRsmfT449Idd5idDAAAWIjTlbvRo0dr7dq1WrFihS5fvqxGjRqpe/fuWrhwoR544AGz4/0yhw/bV8F+9pk0Zoy0cKH9ThMAAAAO5nTlLjY2VrGxsWbHcIyrV6XERPtZuqAgac8eqW9fs1MBAAALc7pyZwmGIb3zjjR9ulRQIL3yiv1uE/XqmZ0MAABYnNMtqKgLhmHo6adnqVZ2gcnKkv7wBykuzr5o4quvpOeeo9gBAIA64ZblLiMjQ3/96zJlZmY67kWLi6VZs6QuXaSzZ6UPPpC2bpXuustxvwcAAMDPcMtyt2LFu6qoeE0rVrzrmBdMTZWCg6XXX7cXvCNHpMhIx7w2AADAbXCbcpeY+Ir8/YPUsWOk3n8/W9IkpacfV4cOg+XvH6TExFdu/0VPn5aioqRhw6R77pG+/FL6058kb29HxwcAALglblPuZs9+RrNnP6fCQk9duLBNkk0XLmxTUZGXZs9+TrNnP3PrL3b9ujRvnv1s3Rdf2C+/fvCB1KFDreUHAAC4FW5T7ry8vBQfP0F+frYa435+NsXHT5CXl9etvdA//2mfV/enP0nx8fYFE8OHSzbbzz4VAACgtrlNubuhqqpC3t5r1KrVAHl7r1FVVcWtPTE3134P2IEDpVat7GfsXn1V8vWt3cAAAAC3we3KXbdu7bRggU0nT27XggU2devW7qefUF5uXyjRubP00UfShg32zYhDQuomMAAAwG1wu02M33nnr9W/jo+foPj4CT/+4H/9y37bsKNH7f+dO1f6zW/qICUAAMAv43Zn7m7JxYvSo49K999vX/l68KD9FmIUOwAA4OTc7szdT6qslFavlmbOtC+QWLlSmjhR8qADAwAA10BruSEjw367sCeesK9+PX5cmjyZYgcAAFwKzeW776QpU6TwcPv+dZ9+Kq1dK7VoYXYyAACA2+a+l2UNw77y9ZlnpNJS+4rYqVMlT/f9XwIAAFyf0525O3TokKZOnaqQkBD5+vqqbdu2io2N1YkTJxz3m3z5pdSnjzRunNS/v5SVJT39NMUOAAC4PKdrMwsWLNC+ffsUExOjLl266Pz580pKSlK3bt20f/9+hfya/eUKC6WXXpIWL7bfKmzXLnu5AwAAsAinK3czZsxQeHi4PL93Fi02NlahoaGaP3++NmzYcPsvahhSSoqUkCD9+9/2gjdjhlS/vgOTAwAAmM/pLstGRETUKHaS1KFDBwUHBysrK+v2X/DECWnwYGnUKCksTDp2TJo1i2IHAAAsyenK3c0YhqELFy6oefPmt/6kkhIpMVG6914pO1tKTZW2b5cCA2stJwAAgNlcotxt2rRJ586dU2xs7K094f337fd+XbBAevZZ++3DoqJqNyQAAIATcLo5d/8rKytLU6ZMUa9evTRu3LiffvDZs9JTT9nP0A0YIO3cKXXqVDdBAQAAnIBTl7u8vDwNHTpUTZo0UUpKimw2248/eN06++bDTZpIW7ZII0fabyEGAADgRpy23BUUFCgyMlJXr17VJ598opYtW/7k4xOSkuR3991SUJB9c+INGxQXF6e4uLg6SgwA1pWcnKzk5OQaYwUFBSalAfBTnLLclZaWKioqSidPntSuXbvUuXPnn33O4rffVrdbnZMHALgtN/uwnJmZqbCwMJMSAfgxTlfuKisrFRsbqwMHDmj79u3q0aPHrT2xY8faDQYAAOACnK7czZgxQ2lpaYqKilJ+fr42btxY4+djxowxKRkAAIDzc7pyd/jwYdlsNqWlpSktLa3Gz2w2G+UOAADgJzhduduzZ4/ZEQAAAFyWS2xiDAAAgFtDuQMAALAQyh0AAICFUO4AAAAshHIHAABgIZQ7AAAAC6HcAQAAWAjlDgAAwEIodwAAABZCuQMAALAQyh0AAICFUO4AAAAshHIHAABgIZQ7AAAAC6HcAQAAWAjlDgAAwEKcstxdu3ZNc+bM0eDBg9W0aVN5eHho/fr1ZscCAABwek5Z7i5duqS5c+fq+PHj6tq1qyTJZrOZnAoAAMD5OWW5CwgIUF5enk6fPq2FCxeaHcepJCcnmx2hzrjLsbrLcUocKwDUBacsd/Xq1ZO/v78kyTAMk9M4F3f6B8NdjtVdjlPiWAGgLjhluQMAAMAvQ7kDAACwEModAACAhXiaHeDXKikpkSR99dVXJiepGwUFBcrMzDQ7Rp1wl2N1l+OUOFarufG+e+N9GIBzcPlyd+bMGUnSmDFjzA1Sh8LCwsyOUGfc5Vjd5TgljtWKzpw5o9///vdmxwDwHy5f7gYNGqSNGzcqMDBQ3t7eZscBALdRUlKiM2fOaNCgQWZHAfA9Ll/umjdvrocfftjsGADgljhjBzgfpy13SUlJunLlis6dOydJSk1NVU5OjiRp2rRpaty4sZnxAAAAnJLNcNJdgtu1a6ezZ89K+u+txwzDkM1m0+nTp3XXXXeZGQ8AAMApOW25AwAAwO1jnzsAAAALcdlyd/36dT3//PMKCAiQj4+PevbsqV27dpkdy+GuXbumOXPmaPDgwWratKk8PDy0fv16s2PVikOHDmnq1KkKCQmRr6+v2rZtq9jYWJ04ccLsaA519OhRxcTEqH379mrYsKFatGihPn36KD093exodWLevHny8PBQaGio2VEcau/evfLw8Ljp18GDB82O53CZmZmKjo5Ws2bN1LBhQ4WGhmrZsmVmxwIgJ15Q8XPGjx+vrVu3KiEhQR07dtS6des0ZMgQ7dmzx1Krty5duqS5c+eqbdu26tq1q/bu3Vs9B9FqFixYoH379ikmJkZdunTR+fPnlZSUpG7dumn//v0KCQkxO6JD5OTkqKioSOPHj1dAQICKi4uVkpKi6OhorVy5UpMmTTI7Yq359ttv9fLLL6thw4aW/Xv81FNPKTw8vMZY+/btTUpTOz788ENFRUUpLCxMiYmJ8vX11cmTJ5Wbm2t2NABy0Tl3Bw8eVM+ePbVo0SJNnz5dkv1M3r333it/f399+umnJid0nLKyMl25ckX+/v7KyMhQeHi43nzzTY0dO9bsaA63b98+hYeHy9Pzv585Tp48qdDQUI0cOVIbNmwwMV3tqqqqUlhYmEpLSy19t5XRo0fr8uXLqqioUH5+vo4cOWJ2JIfZu3ev+vXrp5SUFA0fPtzsOLXm6tWr6tSpk+677z6lpKSYHQfATbjkZdmUlBR5enpq8uTJ1WP169fXhAkTtG/fPkt9eqxXr578/f0l2VcLW1lERESNYidJHTp0UHBwsLKyskxKVTc8PDzUpk0bFRQUmB2l1nz88cfaunWrlixZUr3y3YoMw1BhYaEqKirMjlIrNm/erIsXL2revHmS7FNHqqqqTE4F4Ptcstx9/vnn6tSpk3x9fWuM37gU8sUXX5gRC7XAMAxduHBBzZs3NzuKwxUXFys/P1+nTp3S4sWLtXPnTvXv39/sWLWisrJS8fHxmjRpkmUur/+YRx99VH5+fvL29la/fv2UkZFhdiSH2rVrlxo3bqxvvvlGQUFBatSokfz8/PTkk0/q+vXrZscDIBedc3f+/Hm1atXqB+M3xm5sfAzXt2nTJp07d05//vOfzY7icNOnT9eqVask2c/cjRgxQklJSSanqh1vvPGGcnJytHv3brOj1Jr69etr5MiRGjJkiJo3b66jR49q0aJFuv/++/XZZ5+pa9euZkd0iBMnTqiiokIPPvigJk6cqAULFmjPnj1atmyZrly5os2bN5sdEXB7LlnuSkpKVL9+/R+MN2jQoPrncH1ZWVmaMmWKevXqpXHjxpkdx+ESEhI0atQo5ebmasuWLaqoqLDkmY/Lly8rMTFRiYmJatasmdlxak1ERIQiIiKqv3/ggQc0cuRIdenSRTNnztSOHTtMTOc4RUVFKi4u1hNPPKElS5ZIkh588EGVlZVp5cqVeumll9ShQweTUwLuzSUvy3p7e9/0H8HS0tLqn8O15eXlaejQoWrSpIlSUlIsOT8rKChI/fr10yOPPKK0tDQVFRUpKirK7FgO98ILL6h58+aKj483O0qda9++vYYNG6Y9e/ZYZs7sjffXuLi4GuM3vt+/f3+dZwJQk0uWu1atWt300uv58+clSQEBAXUdCQ5UUFCgyMhIXb16VTt37lTLli3NjlQnRowYoUOHDllqX78TJ05o9erVio+P17fffqszZ87ozJkzKi0tVVlZmc6ePavvvvvO7Ji1qk2bNiorK9O1a9fMjuIQN95f77zzzhrjNxZ+Wf3PE3AFLlnufvvb3yo7O1uFhYU1xg8cOCBJlpnb4o5KS0sVFRWlkydPKj09XZ07dzY7Up25MZ3ASitmc3NzVVVVpWnTpunuu++u/jp48KCys7PVrl07zZ071+yYterrr7+Wt7f3DxaAuaru3btLsu9Z+H03PnC3aNGizjMBqMkly93IkSNVWVlZPRldsu9zt27dOvXs2VOtW7c2MR1+qcrKSsXGxurAgQN699131aNHD7Mj1YpLly79YKy8vFxvvfWWfHx8FBwcbEKq2hEaGqpt27bpvffeq/7atm2bQkJC1LZtW7333nuaMGGC2TEd4mZ/rocPH1ZqaqoGDhxoQqLaMWrUKEnS2rVra4yvWbNGXl5e6tu3rwmpAHyfSy6o+N3vfqeYmBjNnDlTFy9eVPv27bV+/Xrl5ORo3bp1ZsdzuKSkJF25cqX6k3FqaqpycnIkSdOmTVPjxo3NjOcwM2bMUFpamqKiopSfn6+NGzfW+PmYMWNMSuZYkydPVmFhoXr37q2AgADl5eVp06ZNys7O1muvvSYfHx+zIzpMs2bNNGzYsB+ML168WJIUHR1d15FqTWxsrHx8fBQRESF/f38dO3ZMq1atkq+vr+bPn292PIfp2rWrHnvsMf3tb39TRUWFevfurb179yolJUWzZs1ym2kUgFMzXFRpaanx7LPPGq1atTIaNGhg9OjRw/jwww/NjlUrAgMDDZvNZthsNsPDw8Pw8PCo/vXZs2fNjucwffv2rT62//3y8PAwO57DvP3228aAAQOMli1bGl5eXkbTpk2NgQMHGmlpaWZHqzN9+/Y1QkNDzY7hUEuXLjV69OhhNGvWzPDy8jJat25tjB071jh16pTZ0RyuvLzcePHFF43AwECjXr16RqdOnYy//OUvZscC8B8uefsxAAAA3JxLzrkDAADAzVHuAAAALIRyBwAAYCGUOwAAAAuh3AEAAFgI5Q4AAMBCKHcAAAAWQrkDAACwEModAACAhVDuAAAALIRyBwAAYCGUOwAAAAuh3AEAAFgI5Q4AAMBCKHcAAAAWQrkDAACwEE+zAwDuqLy8XIsXL5bNZlNGRobmzZun9evXq6ioSN27d9dDDz1kdkQAgIui3AEmWL58uUaNGqXAwEA999xzGjRokI4cOaIBAwbo3LlzlDsAwC9GuQNM4OnpqcDAQElSfn6+oqOj5e3treXLl6t169bmhgMAuDSbYRiG2SEAd3bPPfdo3rx5Gj58uNlRAAAWQLkDTHTx4kW1bNlSeXl58vf3NzsOAMACWC0L1LHy8nLt3r1bkvTJJ58oMDCwutjt3LlTWVlZZsYDALg4yh1Qx9asWaMhQ4aopKREO3bsqC52ZWVl2r17tzp37mxyQgCAK+OyLFDHvvzyS82dO1dBQUEaMWKEVq1apaZNm8owDE2bNo3LswCAX4VyBwAAYCFclgUAALAQyh0AAICFUO4AAAAshHIHAABgIZQ7AAAAC6HcAQAAWAjlDgAAwEIodwAAABZCuQMAALAQyh0AAICFUO4AAAAshHIHAABgIf8Ps1JveHgkQwoAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x323354950>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated slope     = 1.1973201705016694\n",
      "Estimated intercept = 1.9545255532740005\n"
     ]
    }
   ],
   "source": [
    "# Compute synthetic coordinates (xi,ti) and then use them to estimate the slope and intercept of \n",
    "# of the regression line\n",
    "\n",
    "xi = [0.22, 1.2, 1.9, 2.3, 3.1, 3.9, 4.1, 5.4]\n",
    "ti = 1.2*xi + 2.+ .4*randn(size(xi))\n",
    "\n",
    "# Estimate the slope and intercept \n",
    "\n",
    "(a_estimated,b_estimated) = fit_a_line(xi,ti;flag=false)\n",
    "\n",
    "# Use estimated slope and intercept to compute a regression line \n",
    "\n",
    "x = collect(0:0.1:maximum(xi))\n",
    "\n",
    "t = a_estimated*x + b_estimated\n",
    "\n",
    "# Plot observations and regression line \n",
    "\n",
    "subplot(221);plot(xi,ti,\"*b\",label=\"Observations\"); \n",
    "subplot(221);plot(x ,t ,\"-r\",label=\"Regression\"); \n",
    "\n",
    "title(L\"Regression with $b \\neq 0$\"); \n",
    "xlabel(L\"$x$\"); \n",
    "ylabel(L\"$T$\")\n",
    "legend(bbox_to_anchor=[1.05,1],loc=2,borderaxespad=0)\n",
    "println(\"Estimated slope     = \",a_estimated)\n",
    "println(\"Estimated intercept = \",b_estimated)\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}