Falanga-HW3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions used in this work\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newtonsMethod(interval, guess, f, h, maxIterations = 10, samples = None, debug = False):\n",
    "    '''\n",
    "    newtonsMethod - finds the root of f through estimation via evaluation of\n",
    "    the tangent lines of that function.\n",
    "    f - The function to evaluate, g'(x)  (must be callable)\n",
    "    h - The callable version of g'/g''  (must be callable)\n",
    "    interval - The interval upon which f() is continuously differentiable.\n",
    "    guess - An initial starting point\n",
    "    maxIterations - The number of times to iterate, or loop, looking for the root\n",
    "    '''\n",
    "    if not callable(f) or not callable(h):\n",
    "        raise ValueError(\"f and/or h must be callable objects\")\n",
    "    if not interval[0] <= guess <= interval[1]:\n",
    "       \n",
    "        raise ValueError(\"Initial guess lies outside of the provided interval\")\n",
    "        \n",
    "    iteration = 0\n",
    "\n",
    "    # Contains a list of tuples which may be treated as x,y pairs\n",
    "    pointsList = []\n",
    "    x = guess\n",
    "    while iteration < maxIterations:\n",
    "        y = f(x) if samples is None else f(samples, x)\n",
    "        if debug == True: print(\"DEBUG: x={}, y={}\".format(x, y))\n",
    "        #if np.isclose(y, 0.0, atol = 1e-2):\n",
    "        if np.isclose(y, 0.0):\n",
    "            break\n",
    "        elif iteration >= maxIterations:\n",
    "            raise ArithmeticError(\"No root found for f within max iterations\")\n",
    "\n",
    "        pointsList.append((x, y))\n",
    "        x = x + h(x) if samples is None else h(samples, x)\n",
    "        if x < interval[0] or x > interval[1]:\n",
    "            raise ValueError(\"There is no solution in this interval\")\n",
    "        iteration += 1\n",
    "\n",
    "    # If we got here, a root was found\n",
    "    root = x\n",
    "    return [(root, iteration), pointsList]\n",
    "\n",
    "def bisectionMethod(a, b, g, maxIterations = 10, guess = None, samples = None, debug = False):\n",
    "    # The first check, if g(a) * g(b) > 0, there is no root\n",
    "    atA = g(a) if samples is None else g(samples, a)\n",
    "    atB = g(b) if samples is None else g(samples, b)\n",
    "    if debug == True: print(\"atA:{}; atB:{}\".format(atA, atB))\n",
    "    if atA * atB > 0.0:\n",
    "        raise ValueError(\"No 0 in [{}, {}]\".format(a, b))\n",
    "        \n",
    "    # Generate a guess by means of\n",
    "    #  x(t) = (1/2)*(g(a) * g(b))\n",
    "    nextX = lambda: (1/2)*(a + b)\n",
    "\n",
    "    found = False\n",
    "    iteration = 0\n",
    "    x_t = nextX() if guess is None else guess\n",
    "    while iteration < maxIterations:\n",
    "    \n",
    "        # Ok, the checks are now done, perform the algorithm\n",
    "        candidateVal = g(x_t) if samples is None else g(samples, x_t)\n",
    "\n",
    "        # using default values for isclose() as this defaults to a relative\n",
    "        # tolerance of 1*10^-9 and an absolute tolerance of 0.0\n",
    "        if np.isclose(candidateVal, 0.0):\n",
    "            found = True\n",
    "            break\n",
    "        \n",
    "        # We haven't found a root, so shrink the interval\n",
    "        atA = g(a) if samples is None else g(samples, a)\n",
    "        if atA * candidateVal <= 0:\n",
    "            b = x_t\n",
    "        else:\n",
    "            a = x_t\n",
    "\n",
    "        x_t = nextX()\n",
    "        iteration += 1\n",
    "\n",
    "    if found:\n",
    "        # This little bit of apparent craziness makes quirks with floating point\n",
    "        # representation look correct when we look at it.\n",
    "        root = abs(x_t) if x_t == -0.0 else x_t\n",
    "        return (root, iteration)\n",
    "    else:\n",
    "        raise RuntimeError(\"Unable to find a root for g in [{}, {}]\".format(a, b))\n",
    "        \n",
    "def fixedPointMethod(a, b, gp, L, samples = None, guess = None, debug = False):\n",
    "    G = lambda samples, x: gp(samples, x) + x\n",
    "    \n",
    "    x = (1/2)*(a + b) if guess is None else guess\n",
    "    iteration = 0\n",
    "    contracting = True\n",
    "    while contracting:\n",
    "        # Generate a candidate root value\n",
    "        candidate = G(samples, x)\n",
    "        \n",
    "        # First, check that we've fully converged\n",
    "        if np.isclose(candidate, x):\n",
    "            break\n",
    "            \n",
    "        # Check that the function is contractive\n",
    "        if debug == True:\n",
    "            print(\"a:{}; b:{}; candidate:{}; x:{}\".format(a, b, candidate, x))\n",
    "        if not a <= candidate <= b or not a <= x <= b:\n",
    "            contracting = False\n",
    "            raise RuntimeError(\"G(x) is no longer contractive on interval\")\n",
    "            \n",
    "        # condition 2 for contractiveness\n",
    "        temp = G(samples, candidate)\n",
    "        np.abs(candidate - temp) <= L*np.abs(x - candidate)\n",
    "        \n",
    "        # keep trying\n",
    "        x = candidate\n",
    "        iteration += 1\n",
    "        \n",
    "    return (candidate, iteration)\n",
    "\n",
    "def secantUpdate(x0, x1, f, samples):\n",
    "    if samples is None:\n",
    "        t = x1 - f(x1)*(x1 - x0)/(f(x1) - f(x0))\n",
    "    else:\n",
    "        t = x1 - f(samples, x1)*(x1 - x0)/(f(samples, x1) - f(samples, x0))\n",
    "\n",
    "    return x1, t\n",
    "\n",
    "def secantMethod(x0, x1, f, maxIterations = 10, samples = None, debug = False):\n",
    "    update = lambda x0, x1: x1 - f(x1)*(x1 - x0)/(f(x1) - f(x0))\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < maxIterations:\n",
    "        test = f(x1) if samples is None else f(samples, x1)\n",
    "\n",
    "        if np.isclose(test, 0):\n",
    "            break\n",
    "\n",
    "        x0, x1 = secantUpdate(x0, x1, f, samples)\n",
    "\n",
    "        if debug:\n",
    "            print(\"x0 {}; x1 {}; test {}\".format(x0, x1, test))\n",
    "\n",
    "    else:\n",
    "        raise RuntimeError(\"Unable to converge to solution\")\n",
    "\n",
    "    return (x1, iteration + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework Problem 1: Iteratively Reweighted Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, read the dataset\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipynb.fs.defs.Falanga_HW2 import newtonsMethod\n",
    "df = pd.read_table('datasets/facerecognition.dat')\n",
    "\n",
    "# Now, some nicer lambdas for working with this illdefined dataset\n",
    "weirdKey = 'match eyediff nosecheekdiff variabilityratio'\n",
    "\n",
    "getMatch = lambda i: int(df[weirdKey][i].split()[0])\n",
    "getEyeDiff = lambda i: np.float(df[weirdKey][i].split()[1])\n",
    "getNoseCheekDiff = lambda i: np.float(df[weirdKey][i].split()[2])\n",
    "getVariabilityRatio = lambda i: np.float(df[weirdKey][i].split()[3])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework Problem 2: Problem 2.1 Pg. 54\n",
    "**Cauchy Distribution**: C(\\theta, 1) with X_i = 1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71, -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEgCAYAAACegPWEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xV9f3H8dcnGwJhJQHC3ks2IiAoBQS1iooLZ504W21rtXap7c/a2larVavWvffAgQqI4GCFvfdIIIRNwkjI+Pz+OCf1mia5N8m9Obn3fp6Px30k95x773mfKPnkfL/f8/2KqmKMMcbUVIzXAYwxxoQ3KyTGGGNqxQqJMcaYWrFCYowxplaskBhjjKkVKyTGGGNqxQqJMcaYWokL5EUiEgP0BzKAY8AqVc0NZTBjjDHhQaq6IVFEugB3AeOADcAeIAnoDhwFngJeVNXS0Ec1xhhTH/krJK8D/wa+1nIvFJF04FLggKq+GNKUxhhj6q0qC4kxxhjjj9/OdhFJcZu4ym/vF5pIxhhjwkmVhURELgLWAu+KyCoROdFn9wuhDGaMMSY8+Lsi+Q0wWFUHAFcDL4vIJHefhDSZMcaYsOBv+G+squYAqOoCEfkR8LGItAWsc8UYY4zfK5J83/4Rt6iMBs4B+oQwlzHGmDDhr5DcRLkmLFXNB04HrglVKGNqS0TuFZFX3O/bi8hhEYl1n38lItfV4DOvEpFvfJ4fFpHO7vcviMj/BSt/FRlGi0h2FftPFpENbrZzQ53H57ijRGRdXR3P1C9VFhJVXaaqGyvYXqSqr4YulolkIrJVRMbV1fFUdbuqNlLVkiB/biNV3RzMzwyCPwKPudk+CNVBRERFpGvZc1X9WlV7hOp4pn4LaK4tEbkqxDmMMcHRAVjldQgTXQK5j+QPOFOkGBNyInK9iGwUkf0iMlVEMnz2jReRdSJySESeEJHZgTRRiUhH9y/o/xlcIiKtRWS5iNzhPm8iIs+KSI6I7BCR/ytrEqvgvT/4qxxoJiKfiEi+iMz37V8UkREistDNvlBERvjsy3DPdb977tf77GvgNpsdEJHVgO8Q/PJ5NgGdgY/cpq3E8ld/5Zr8yn4uPxGR7SKyV0R+6/PaWBH5jYhscs9pkYi0E5E57kuWuce5uHyTm4j0cpsQD7q3Dkz02feCiDxe2c/KhB9/95E8DfQCrqybOCaaicgY4AHgIqA1sA14w92XCrwD3A20ANYBIyr+pICP1xGYjdMU9Hd384tAMdAVGAiMBwLtT7kEuA9oBmwE7neP0xz4BHjUzf4Q8ImItHDf9zqQjTMp6gXAn0VkrLvvHqCL+5gA/KSyg6tqF2A7cLbbtFUYYO6RQA9gLPAHEenlbv+Fe05nAik4/aJHVfUUd39/9zhv+n6YiMQDHwFfAOnAT4FXRcS36avCn5UJT/6uSC4F/mSTMpo6chnwnKoudn8J3g0Md3/hn4kz6/R7qlqM80t5Vy2O1Rv4CrhHVZ8GEJGWwBnA7ap6RFV3Aw8DkwP8zPdUdYGb71VggLv9x8AGVX1ZVYtV9XWcG33PFpF2OL/I71LVAlVdCjwDXOG+9yLgflXdr6pZ7nkH232qekxVlwHLcGb6BqeA/k5V16ljmaruC+DzhgGNgL+o6nFV/RL4GKd4lKnsZ2XCkL9Cchbwll12mjqSgXMVAoCqHgb2AW3cfVk++xTnr/iaugzYgXOVU6YDEA/kuE0yB3FmuE4P8DN9C9tRnF+mUO68XNv4/rz2u6Mhy+8re29WuX3BVlnudsCmGnxeBpBV7g9Q33Oq6pgmDPkbtfUVzl9jr9RJGhPtduL8MgdARJJxmoJ2ADlAW5994vu8Bu4F9gKv+fSBZAGFQKqqNnUfKapa23umfnBervY457UTaC4ijSvYB855tyu3rzqOAA19nreqxnuzcJrUqmsn0E6cdYzK+J6TiTB+O9tVdSVwfh1kMdElXkSSfB5xwGvA1SIyQEQSgT8D81V1K04fQ18ROdd97S1U75dieUXAhUAyztQ/Me4Nt18A/xBnstIYEekiIqfW4jgAnwLdReRSEYkTkYtxmtY+dpurvgMecH8O/YBrcZp7AN4C7haRZuLMKPHTah57KTBZROJFZAhOH0ygngH+JCLdxNHPp18nF6djvyLzcQrYne5xRwNn4/Z3mcgT0PBfVd0Z6iAm6nyKs9pm2eNeVZ0J/B54F+cv8S64/ROquhfnF/+DOM1dvYFMnCuIGlHV48AknKar59y/oK8EEoDVwAGcpq/WNT2Ge5x9OM3Ev3Sz3wmc5Z4TOH0HHXH+kn8fp99murvvPpxmoS04Re7lah7+9zg/xwPuZ71Wjfc+hFPIvgDygGeBBu6+e4EX3SbAi3zf5P5cJ+L0N+0FngCuVNW11cxuwoTf9UhEJAnnH8Eovl9qdyXwiaraeHXjCfeXfjZwmarO8jqPMdGsykkbReRenEvSr3AuV3fz/VK7f3GLzC9VdXloYxoDIjIB5//DY8CvcKbvmedpKGOM39l/F6rqvZXse0ic5Xar2/lnTE0Nx2maKWt6OldVj3kbyRhjS+0aY4ypFX9XJACISBpwF04HZ1LZdlUdE6JcxhhjwkRAhQRnKOKbOHfo3ogzTcOeUIUKttTUVO3YsaPXMYwxJmwsWrRor6qmBfLaQAtJC1V9VkRuU9XZwGwRmV3ziHWrY8eOZGZmeh3DGGPChogEPItCoIWkyP2aIyI/xhnvXpu7io0xxkSIQAvJ/4lIE5wbqv6FMxPo7SFLZYwxJmwEWkgOqOoh4BDwI3CW9AxZKmOMMWEjoClScK5CAtlmjDEmyvi7s304zuJBaSLyC59dKUCFq8YZY4yJLv6athJw1gmIA3ynuc6jerOIGmOMiVBVFhKfob4vqOo2EUlW1SN1lM0YY0wYCLSzPUNEpuFcnbQXkf7ADap6c+iimfrg0NEivt20lx0HjqEo3Vo2ZlinFjRIsJZNY4wj0ELyT2ACMBVAVZeJyCkhS2U8t25XPo/MXM8Xq3IpLv3hfGwpSXFcfGI7bjy1Cy0aJXqU0BhTXwRaSFDVLGd10/8qCX4c47XjxaU8NH09T83ZRHJCHFeN6MgZfVvTrWUjSkuVFTsO8VZmNs9+s4XXF2Rx46mduWZkJxomBPy/kjEmwgT6rz9LREYAKiIJwM+ANaGLZbyQV1DENc8vJHPbASaf2I5fn9GTpg0TfvCaUd3SGNUtjdvGduWvn63j71+s54XvtnLp0PZMHNCGLmnJlPuDA4DiklJ25xeSc+gY+w4fp0mDeHq0avw/n2+MCT8BTSMvIqnAI8A4nMWEvgBuc5cQrfeGDBmiNtdW1Q4dK+LyZ+azJiePhy4ewMT+GQG9L3Prfv791Sa+XLcbVWiRnEDb5g1JjIuhqKSUopJS9h8+Tm5+ISXlmsgS4mI4q19rfnNmL1KticyYekVEFqnqkIBeGw3rkVghqVppqXLtiwv5esNenrpiMGN7taz2Z2QfOMrXG/aydPtBdh46RlFJKfGxMcTHxtC0YTwZTRrQumkSGU0a0KJRAgePFjF9dS5vZWbRukkSL11zEu1bNAzB2RljaiLohcRdj+R6oCM+zWGqek0NM9YpKyRVe2j6eh6duYE/nXsCVwzrUKfHXrTtANe+uJDGSXF8dtspJCdaX4sx9UF1CkmgU6R8CDQBZgCf+DxMmFu54xCPfbmBSYPacPlJdb9q8uAOzXjmyiFkHzjGA9Os282YcBTon38NVfWukCYxda64pJS73l1Oi0aJ3HN2nwo7yevCkI7NuW5kJ/7z9RYm9m/D0E7NPclhjKmZQK9IPhaRM0OaxNS51xdmsWpnHvee3YcmDeI9zfLL8T1Ib5zIIzPXe5rDGFN9VRYSEckXkTzgNpxickxE8ny2mzB17HgJ/5q5gaEdm3Nm31ZexyEpPpbrR3Xm2437WLL9gNdxjDHVUGUhUdXGqprifo1R1QY+z1PqKqQJvpfnbWV3fiF3TOjhWZNWeZee1J6mDeN5fNZGr6MYY6ohoKYtEZkZyDYTHgqKSnhq9mZO6Z5Wr/ojkhPjuGJYB2au3c3Og8e8jmOMCZC/pq0kEWkBpIpIMxFp7j46AoHdsVYDIjJAROaJyFIRyRSRoT777haRjSKyTkQmhCpDJPto2U72HTnOjad09jrK/7hoSDtU4d1F2V5HMcYEyN8VyQ1AJtATWOTz+BB4PIS5HgTuU9UBwB/c54hIb2Ay0Ac4HXhCRGwa2mpQVV74bivdWzZieJcWXsf5H+2aN2R45xa8vSib0tLIv1nWmEjgr4/kEVXtBNyhqp1VtZP76K+qj4lIqK5KFGcVRnDuX9npfn8O8IaqFqrqFmAjMLSC95tKZG47wKqdefxkRMd60zdS3kUntmX7/qPM37Lf6yjGmAAE1Eeiqv8CEJE4ERktIn8RkfnAf0KU63bgbyKSBfwduNvd3gbI8nldtrvNBOjtzCySE2I5b2D9/bGd3qc1DRNi+Wj5Tv8vNsZ4zm8hEZH2IjJFRN4H5gLnAl8Bp6rqj2t6YBGZISIrK3icA9wE/FxV2wE/B54te1sFH1Vh+4ebOVNEMvfs2VPTmBGloKiEaSt2cfoJrev1tO8NEmIZ0zOdz1fu+p+JHo0x9U+Vv01EZCFOs9LnwC/c5qSgUNVxVRz3JZx7VwDeBp5xv88G2vm8tC3fN3uV//yngafBmWurtnkjwYw1ueQXFjNpUP29GilzZt/WfLw8h/lb9jGiS6rXcYwxVfB3RXKyqp6jqk+UFREROasOcu0ETnW/HwNscL+fCkwWkUQR6QR0AxbUQZ6I8P7iHbRMSWRY5/rXyV7e6B5pJMXHMG3FLq+jGGP88NfZfryCzX8MURZf1wP/EJFlwJ+BKW6eVcBbwGrgM+AWVbWVGgOQV1DEnA17mNg/g9iY+tnJ7qthQhxjeqYzbeUuG71lTD0X6FxbvkL+W0hVv1HVwe7osJNUdZHPvvtVtYuq9lDVaaHOEim+WreHohJlQh/vp0MJ1Pjerdh7uJBl2Qe9jmKMqUJNCskNQU9hQm766lxaJCcwsH0zr6MEbHSPNGJjhJlrdnsdxRhThUCnSPmTiMQBqOoCEUkRkedDG80Ey/HiUr5au5uxvdLDolmrTNOGCQzp0IwZa3K9jmKMqUKgVyRxwHwR6Sci44GFOHe4mzAwf8s+8guLOa13+DRrlTmtd0vW7sona/9Rr6MYYyoR6A2JdwN3AfOBF4Afq+pjIcxlgmj66lyS4mMY2TX8htGWrR8/065KjKm3Am3aOgV4BGfE1ldAKKdHMUGkqsxYncuobmk0SAi/ack6pSbTJS2ZGdZPYky9FWjT1t+BC1X1AVW9FOdGvy9DF8sEy6qdeew8VMBpvVt6HaXGxvVqyfwt+8grKPI6ijGmAoEWkuGqurrsiaq+B5wcmkgmmL5YnUuMwNie6V5HqbFxvVtSVKLMWW9T3RhTH/lbj+RyEYmp6KY/Vd0nIl1EZGTo4pnamr46l8EdmtGiUaLXUWpsUPtmNGsYb8OAjamn/M3c1wJYIiJl65DsAZKArjhTmOwFfh3ShKbGsvYfZU1OHr85s6fXUWolNkb4Uc90Zq7ZTXFJKXGxNbn9yRgTKn7XIwEGAa8DacBY9/kO4ApVPV9VN1TxEcZDZfdfhOOw3/LG9WrJoWNFLNp2wOsoxphy/M4l7jZrTXcfJoxMX51L1/RGdEpN9jpKrZ3SPY2E2BhmrMnlpDCYdNKYaFKtNgIR+UmogpjgOnS0iPlb9of1aC1fjRLjOKlzcxsGbEw9VN3G5tv8v8TUB7PW7aakVBkfIYUEnLvct+w9wqY9h72OYozxYb2WEWr66lzSGyfSv21Tr6MEzRh3CPOM1XaXuzH1SSBL7c4SkS9FZBbQtex7EbEbEuupwuISvlq3m7G9WhITRpM0+tO2WUN6tU6xYcDG1DOBLNx9lftVgE+Aq0OWxgTF3E37OHK8JKKatcqc1iudx2Zt5MCR4zRLTvA6jjGGAK5IVHWb+9gKFPo83xb6eKYmpq/OpWFCLMO7RN7oprG9WlKqTh+QMaZ+qG4fydZQhDDBU1qqzFiTy6nd00iKD79JGv3p26YJ6Y0TbY0SY+qRahUSVZ0UqiAmOJZkHSQ3r5DxfSKvWQsgJkYY2yudOev3Ulj8PzP3GGM8UONRWyIS3vNuRKhPlueQEBfDuF6RWUjAucv9cGEx8zfv9zqKMYbaDf/9ImgpTFCUlirTVuZwavc0GifFex0nZE7umkpSfIw1bxlTT1Q5aktEHq1sFxA5NyhEiCVZB8g5VMBdp0f2xWJSfCwju6Yxc81u7puoiETOEGdjwpG/K5KrgZU4M//6PjKB46GNZqrrY7dZa2yv8F17JFDjeqWz4+Ax1uTkex3FmKjn7z6ShcBKVf2u/A4RuTckiUyNlJYq01bsivhmrTJj3GI5Y00uvTNSPE5jTHTzd0VyAbC0oh2q2in4cUxNLd5+gF15BZzVr7XXUepEeuMkBndoxifLc7yOYkzU87ceyX5VPVpXYUzNfd+sFbmjtcqb2D+Ddbn5rNtlzVvGeMnfUrsficjZIvI/bSUi0llE/igi14QunglE2Wit0d3TaJQYyKw3keHMvq2JEZi6bIfXUYyJav6atq4HRgFrRWShiHzqTtq4GXgKWKSqz4U8panSvC37yM0r5Kz+GV5HqVNpjRM5uWsqHy3LQVW9jmNM1Kryz1dV3QXcCdwpIh2B1sAxYL01edUf7yzKpnFiXERO0ujP2f0zuPOd5SzNOsjA9s28jmNMVAr4hkRV3aqqc1V1qRWR+uNwYTHTVuzirP4ZETm3lj8T+rQiITaGqct2eh3FmKgVUCERkXwRySv3yBKR90Wkc6hDmsp9uiKHY0UlXDC4rddRPNGkQTyje6Tx8fIcSkqtecsYLwR6RfIQ8CugDdAWuAP4D/AGYH0kHnpnUTadU5MZ1D56JxqYOCCDPfmFzN+8z+soxkSlQAvJ6ar6lKrmq2qeqj4NnKmqbwLWMO2RbfuOsGDLfs4f3DaqpwkZ27MlyQmxfLjUmreM8UKghaRURC4SkRj3cZHPPmtP8Mi7i3cgApMGtfE6iqcaJMRyZt/WfLx8J0cKi72OY0zUCbSQXAZcAex2H1cAl4tIA+DWEGUzVSgtVd5dlM3Irqm0btLA6ziemzy0HUeOl9id7sZ4IKBCoqqbVfVsVU11H2er6kZVPaaq39TkwCJyoYisEpFSERlSbt/dIrJRRNaJyASf7YNFZIW771GJ4vacbzbuZcfBY1HbyV7eoPbN6JreiDcWbvc6ijFRJ9BRW23dEVq7RSRXRN4Vkdr+BlsJTALmlDtWb2Ay0Ac4HXhCRMrGtf4bmAJ0cx+n1zJD2Hpp7jZaJCdw+gmtvI5SL4gIFw9px+LtB9mQa1OmGFOXAm3aeh6YCmTgjNz6yN1WY6q6RlXXVbDrHOANVS1U1S3ARmCoiLQGUtx7WRR4CTi3NhnCVdb+o8xcm8slQ9uTGBd9945U5rxBbYiPFd5cmOV1FGOiSqCFJE1Vn1fVYvfxApAWokxtAN/fBNnutjbu9+W3R51X529HgEtPau91lHoltVEip/VuyXtLdth67sbUoUALyV4RuVxEYt3H5YDfQfsiMkNEVlbwOKeqt1WwTavYXtmxp4hIpohk7tmzx1/UsFFQVMKbC7czvncrMppaJ3t5Fw1px/4jx5mxerfXUYyJGoFOFXsN8BjwMM4v7+/cbVVS1XE1yJQNtPN53hbY6W5vW8H2yo79NPA0wJAhQyJmiPLHy3M4cLSIK4d38DpKvTSqWxoZTZJ4bcE2fhwla7MY47VAR21tV9WJqpqmqumqeq6qbgtRpqnAZBFJFJFOOJ3qC1Q1B8gXkWHuaK0rgQ9DlKHeemnuVrqmN2J4lxZeR6mXYmOEy4Z14NuN+1hvne7G1Ikqr0hE5F9U0Xykqj+r6YFF5DzgXzh9LZ+IyFJVnaCqq0TkLWA1UAzcoqplDd43AS8ADYBp7iNqLNp2gOXZh/jjOX2i+k52fy4Z2p5HZ27ghe+28ufz+nodx5iI569pKzNUB1bV94H3K9l3P3B/BdszgRNClam+e3L2Jpo2jOf8QXbvSFWaJydw7oA2vLc4mzsn9KBpwwSvIxkT0fytR/JiXQUxVdu4O5/pq3P52dhuJEfRKog1dfXIjryZmcUbC7O48dQuXscxJqIFvB5JGRE5LRRBTNWenL2ZpPgYrhrR0esoYaFnqxSGd27BS99tpbik1Os4xkS0ahUSEWkKvOnOsWXqyM6Dx/hgyQ4mn9ie5snWTBOoq07uyM5DBXyxOtfrKMZEtCoLiYhcLCJzRGSZiKwC1gOvAFtFZLmIzBWRS+okaRR79pstAFw3qpPHScLLuF4tade8Ac98vdnWdDcmhPxdkdyHs6DVRJx5rTq4I7V6AGcDt7n7TYgcPHqc1xdsZ2L/DNo2a+h1nLASGyNcN7Izi7cfZMGW/V7HMSZi+Ssk/1HV+aq6TVWzVPUYgKoedLctAF4Nfczo9Z+vN3OsqIQbR1uHcU1cfGI7Uhsl8NisjV5HMSZiVVlIVPUf/j4gkNeYmtl7uJDnv93K2f0y6N6ysddxwlJSfCzXjuzM1xv2sjz7oNdxjIlIAXe2i0gbERkhIqeUPUIZzMCTX22ioKiE28Z18zpKWLt8WHtSkuJ4YtYmr6MYE5ECuiFBRP4KXIxzt3nZXeZKubVETPDsOlTAy/O2MWlQW7qkNfI6TlhrnBTPVSM68uiXG9mQm083u7ozJqgCvSI5F+ihqme6qyOeraoTQxks2j0+ayMlpcptY+1qJBiuOrkTDeJjeeIruyoxJtgCLSSbgfhQBjHfy9p/lDcWbufiE9vRrrmN1AqG5skJXDG8Ax8u3cHG3TaZozHBFGghOQosFZGn3LXSHxWRR0MZLJr97fN1xIhw65iuXkeJKDee2oUG8bE8NH2911GMiSiBTto01X2YEFu8/QBTl+3kp2O60rqJTSAQTM2TE7h2VGcenbmBlTsOcUKbJl5HMiYiBFRIbPLGuqGq/Onj1aQ3TrSJBkPkulGdeGnuVv7+xTpeuHqo13GMiQj+pki5yd8HBPIaE5iPluewZPtB7pjQw2b4DZGUpHhuOrULX63bY3e7GxMk/vpIbhORHiLSXESauisTIiIJ7raewJTQx4x8BUUl/HXaWnq3TrH1RkLsyuEdSW+cyN8+X2tzcBkTBP4KyYPA6zgLXC0DdovI/cAud9s7wOMhTRglnpq9mR0Hj/G7s3oRG2OrH4ZSg4RYbhvXjYVbD/DZyl1exzEm7Plb2Oo54Lmy5yLSCmcG4A6qeiDE2aLG1r1HePyrjZzVrzUjuqR6HScqXDykHS99t40Hpq1lTK90EuNivY5kTNiq1nokqroLuNKKSPCoKn+YuoqE2Bh+f1Zvr+NEjbjYGH53Vi+27z/KC99u9TqOMWGt2iskquoHoQgSrT5dsYs56/fwy/HdaZmS5HWcqDKqWxpjeqbz2Jcb2Xu40Os4xoStgAqJiNjkRCGQX1DEfR+t4oQ2KVwxrIPXcaLSb87sxbGiEh62mxSNqTG/hURE2gAf10GWqPPgZ+vYc7iQ+8/tS1xstS8OTRB0TW/E5cM68PqC7azcccjrOMaEJX/3kfQBPgPuqps40eO7TXt5ed42rjm5E/3bNfU6TlT7+WndaZ6cyG/fX0FJqQ0HNqa6/P0ZPAuYoqrz6iJMtDhSWMxd7y6nY4uG3DG+h9dxol6TBvH87se9WJZ9iNcXbPc6jjFhx18hWQicXxdBosmDn60l+8Ax/nZhfxok2LDT+uCcARkM79yCBz9by55863g3pjr8FZKJQIqIPFgXYaLB3E37eHHuNq4e0YkTOzb3Oo5xiQh/OvcEjhWV8MCna7yOY0xY8bdme4mqTgEO11GeiHboaBG/eGspnVKT+dUEa9Kqb7qmN+KGU7rw3pIdzFm/x+s4xoSNgIYKqeofQx0k0qkqd7+/nD35hTwyeYA1adVTt47pSue0ZO5+bwWHC4u9jmNMWLAxp3XkrcwsPl2xizsm9KBfWxulVV8lxcfytwv6sfPQMR78bK3XcYwJC9UuJCLyVCiCRLJNew5z79TVjOjSgimjOnsdx/gxuENzrh7RiZfmbmP+5n1exzGm3gvkhsROIjJQRIaIyLXApSJykYgMcqeYt6lqq3DseAm3vLqYpPgYHrpoADE2s29YuGNCd9o3b8id7y7n6HFr4jKmKv5uSHwT+AR4CPgrcAowAZgM/ANnGnlbPbESqspvP1jButx8Hr54AK2a2Fxa4aJhQhwPXtCP7fuPcv8nNorLmKr4W4avj6pWNCXtJAD3amRJ0FNFiNcXZPHe4h3cPq4bo3ukex3HVNOwzi24flRnnp6zmTE90xnbq6XXkYypl/w1bU2qaqc6y8tdFLw4kWN59kHunbqK0T3S+NmYbl7HMTX0y/Hd6dmqMXe9u9xmCDamEv7uI/E7JWogr4k2u/MKmPLSItIaJ/Kw9YuEtcS4WP45eQB5x4r59bsrbGleYypQ4+G/IvKH2hxYRC4UkVUiUioiQ3y2nyYii0Rkhft1jM++we72jSLyaH3s6C8oKuH6lxeRV1DEf64cQrPkBK8jmVrq2SqFO0/vwYw1ubz43Vav4xhT79TmPpLrannslThNZ3PKbd8LnK2qfYGfAC/77Ps3MAXo5j5Or2WGoFJV7nxnOcuyDvLwxQPonZHidSQTJNec3ImxPdO5/9M1LMs66HUcY+oVf6O28ip55AMZtTmwqq5R1XUVbF+iqjvdp6uAJBFJFJHWQIqqznX7Zl4Czq1NhmD715cbmbpsJ7+a0IMJfVp5HccEUUyM8I+L+pPeOIlbXlvMoaNFXkcypt7wd0VyEOimqinlHo2BnDrIdz6wRFULgTZAts++bHdbvfBWZhYPTV/PpIFtuHl0F6/jmBBo2jCBxy4dSG5eAXe8s8z6S4xx+SskLwGVrQH7mr8PF5EZIrKygsc5Aby3D869KzeUbargZZX+SxaRKSKSKSKZe/aEdgK+WWt3c/d7KxjVLZW/nN+Peth1Y4JkYPtm3H1GL6avzuXZb7Z4HceYeqHK+0hU9XdV7DZHV1oAAB0JSURBVPO7aqKqjqtJKBFpC7wPXKmqm9zN2UBbn5e1BXaWf6/PsZ8GngYYMmRIyP50XJp1kJtfXUzPVo359+WDSYiz6csi3dUnd2TBlv38ZdpaBrZvxuAOzbyOZIynAv6tJyIpvl9DRUSa4txNf7eqflu2XVVzgHwRGeaO1roS+DCUWfxZtyufq55fQGrjBJ6/+kQaJfq7v9NEAhHhwQv7kdG0ATe/uojcvAKvIxnjqer8+fxVua+1IiLniUg2MBz4REQ+d3fdCnQFfi8iS91H2W3hNwHPABuBTcC0YGSpiU17DnPZM/NJjIvhlWtPIr2xTX8STVKS4nn6ysEcLihmysuLKCgq8TqSMZ6RQDsMRWSxqg4SkSWqOjDEuYJqyJAhmpmZGbTPy9p/lAufnEtRSSlv3jCcrumNgvbZJrx8sWoXU15exHkD2/DQRf2tf8xEDBFZpKpD/L/S1iOptu37jjL56XkcKyrhletOsiIS5cb3acUd47vz/pIdPDVns9dxjPGENepXw8bdh7nsmXkUFpfy6nUn0au13XBo4JYfdWXtrnz++tlaurdsxJieNrmjiS41uSKJysHza3flMfnpuZSUKm9MGcYJbZp4HcnUEyLC3y7oT+/WKfzs9aWs25XvdSRj6lR1ComU+xo1Mrfu5+Kn5hEbI7wxZTg9W9mViPmhBgmx/OfKITRMiOWq5xeQc+iY15GMqTPVKSQXl/saFT5buYvLnplPi+QE3rlxhPWJmEplNG3A81efSH5BMVc9t5BDx2waFRMdAi4kZdPFR9O08S9+t5WbXl1E74wU3rlpBO2aN/Q6kqnn+mQ04cnLB7Npz2FueDmTwmIbFmwin79JG99yv64QkeU+jxUisrxuInojr6CIf3+1ibE9W/LadcNobtPBmwCN7JbK3y7sx7zN+7nj7eWUlkZlt6KJIv5Gbd3mfj0r1EHqm5SkeN65aTitUpKIi7VR0qZ6zhvYll2HCvnrZ2tJb5zI737cy+4xMRHL31xbOe7XbXUTp35p28yaskzN3XhqZ3YdOsaz32whJSme28bZkssmMvlr2vqTvw8I5DXGRCMR4Z6z+zBpUBsenrGe/9gNiyZC+WvaukBEvgH2AKXATlXdLSKdgGZAOnAG8PvQxjQmPMXECA+e34/ColLu/3QNDRJiuXxYZSszGBOe/BWSX+CsB9IUiAU6i8gcnIkWtwFHgPtCmtCYMBcXG8PDFw/gWFEJv/9wJQ0TYpk0qK3/NxoTJvz1kUzDZ4ZdEUnGKSAZqno8xNmMiRgJcTE8cdkgrnlhIXe8vYyEuBjO6ler1aqNqTeqNRxJVY8AZ1oRMab6kuKdu98Hd2jGz15fwodLd3gdyZigqPa4VlVdUH6biNwdnDjGRLbkxDheuHooQzs15+dvLuXdRdleRzKm1oJ1g8SFQfocYyJecmIcz181lBFdUrnjnWW8uXC715GMqZVgFRK708qYamiQEMszPxnCKd3SuOvdFbwyLypv1TIRIliFxOaAMKaakuJjefrKwYzrlc7vPljJk7M3eR3JmBqxKxJjPJQYF8sTlw1mYv8M/jJtLX/+dA2BLn9tTH0RrBUS3w7S5xgTdRLiYvjnxQNo1jCep+dsZv+R4/xlUl+b482EjYAKiXsn+0+Bjr7vUdWJ7tc/hyKcMdEiJka4d2Ifmicn8vCM9Rw8WsRjlw4kKT7W62jG+BXoFckHwLPARzhTpRhjgkxEuG1cN5o3SuAPH67k8mfm8/SVQ2wJA1PvBVpIClT10ZAmMcYAcMWwDrRITuD2N5cy6Ylvef7qoXRKTfY6ljGVCrQR9hERuUdEhovIoLJHSJMZE8XO7Nua168/ibyCYiY98S2ZW/d7HcmYSgVaSPoC1wN/Af7hPv4eqlDGGBjcoTnv3zyCZg0TuPSZ+Xy0bKfXkYypUKBNW+cBnW2OLWPqVocWybx70whueHkRP319CRt3H+a2sd2IibER96b+CPSKZBnOVPLGmDrWLDmBl68bygWD2/LIzA3c+MoiDhcWex3LmP8KtJC0BNaKyOciMrXsEcpgxpjvJcbF8rcL+nHP2b2ZuXY35z3+LVv3HvE6ljFA4E1b94Q0hTHGLxHh6pM70aNlY25+bTETH/uGf106iFO7p3kdzUQ5iYbpGIYMGaKZmZlexzAmaLL2H+X6lzJZn5vPr8/oyfWjOiNi/SYmeERkkaoOCeS1ATVtiUi+iOS5jwIRKRGRvNrFNMbUVLvmDXn3phFM6NOKP3+6llteW0xeQZHXsUyUCqiQqGpjVU1xH0nA+cBjoY1mjKlKcmIcT1w2iF+f0ZPPV+Vy9r++YeWOQ17HMlGoRrPCqeoHwJggZzHGVJOIcOOpXXhjyjAKi0qZ9MR3vDJvm80gbOpUoJM2TvJ5GgMMwdYgMabeOLFjcz752Uh+/tYyfvfBSuZv2c8Dk/rSKDFYE3wbU7lA/y872+f7YmArcE7Q0xhjaqxFo0ReuOpE/j17E//4Yh0rdxzisUsH0iejidfRTIQLtI/kap/H9ap6v6rurs2BReRCEVklIqUi8j8jA0SkvYgcFpE7fLYNFpEVIrJRRB4VG6ZizA/ExAi3/Kgrr10/jCOFxZz7+Lc8NXsTpaXWgGBCJ9BRWw+KSIqIxIvITBHZKyKX1/LYK4FJwJxK9j8MTCu37d/AFKCb+zi9lhmMiUjDOrfgs9tPYUzPdB6YtpZLn5nHzoPHvI5lIlSgne3jVTUPOAvIBroDv6rNgVV1jaquq2ifiJwLbAZW+WxrDaSo6lx1ehJfAs6tTQZjIlnz5ASevHwwD17QjxXZhzj9n3OYahM/mhAItJDEu1/PBF5X1ZDNaS0iycBdwH3ldrXBKWJlst1txphKiAgXDWnHp7eNokt6I372+hJuf2OJ3XNigirQQvKRiKzFGa01U0TSgAJ/bxKRGSKysoJHVR319wEPq+rh8h9XwWsrbfgVkSkikikimXv27PEX1ZiI1qFFMm/fMJzbx3Xjo+U5THh4DrPW1aqb05j/CniKFBFpBuSpaomINAQaA7GqWqtrZRH5CrhDVTPd518D7dzdTXGW9v0D8C4wS1V7uq+7BBitqjf4O4ZNkWLM95ZmHeRXby9jw+7DnD+oLX84qzdNGsb7f6OJKkGfIgVAVQ84ny2jcSZxnAr8p0YJqz7OKFXtqKodgX8Cf1bVx1Q1B8gXkWHuaK0rgQ+DfXxjIt2Adk35+GcjueVHXfhg6Q5Oe3g201fneh3LhDG/hcQdhnuDiLwPzMXp4J4FnKqqP67pgUXkPBHJBoYDn4jI5wG87SbgGWAjsIn/HdVljAlAYlwsv5rQkw9uPpnmyQlc/1Imt72xhANHbO06U31VNm2JyEJgJ/AZ8JmqbqmrYMFkTVvGVO54cSmPz9rI47M2ktIgnt+e2YtJg9rYbMJRrjpNW/4KSUIkLK9rhcQY/9bk5PGb91ewZPtBhnduwf+ddwJd0hp5Hct4JGh9JJFQRIwxgenVOoV3bxzB/eedwKqdhzjjn1/z0PT1FBSVeB3N1HM1mv3XGBOZYmKEy07qwMxfjuaMvq14dOYGznjka77eYEPoTeWqXUhEZFAoghhj6o+0xok8MnkgL187FFXlimcXcOPLi8jaf9TraKYeqskVyTNBT2GMqZdGdUvjs9tP4VcTejB7/R7GPTSbh6av59hxa+4y36tJIbGhHMZEkaT4WG75UVe+vONUJvRxmrvG/uMrPlmeYwtoGaBmhaT8HFjGmCjQukkDHr1kIG/dMJwmDRO45bXFTH56HsuzD3odzXgs4ClSwpkN/zUmuEpKldcXbOfh6evZd+Q4E/tn8KsJPWjXvKHX0UyQBO0+kkhhhcSY0MgvKOKp2Zt55pvNlJbClcM7cOuYrjRtmOB1NFNLVkjKsUJiTGjlHDrGw9PX8/aibFKS4rn1R125ckQHEuNivY5maigkhcSd/TcDOAZsVdXSmkesW1ZIjKkba3fl8cCna5m9fg9tmzXgF6d155wBbYiNsTE64SaYU6Q0AW4BLgESgD1AEtASmAc8oaqzap04xKyQGFO3vtmwlwemrWHVzjy6pCVz+7ju/Lhva2KsoISNYBaS6ThL2n6kqgfL7RsMXAGsUNVna5E35KyQGFP3SkuVz1ft4uEZ61mfe5ierRrz89O6M753S5sQMgxYH0k5VkiM8U5JqfLx8p38c8YGtuw9Qt82TfjFad0Z3SPNCko9Fswrkp6qurayaVFUdXENM9YpKyTGeK+4pJT3l+zg0S83kLX/GIPaN+WnY7sxursVlPoomIXkaVWdIiIV9YOoqo6paci6ZIXEmPrjeHEp7yzK5rEvN7DzUAF9MlK45UddOb1PK+tDqUesaascKyTG1D/Hi0v5YMkO/j17E1v2HqFLWjI3j+7KxAEZxMfaxOReC3ohEZEk4GZgJKDA18CTqlpQm6B1xQqJMfVXSany6YocHp+1kbW78mnbrAE3nNqFCwe3JSne7kPxSigKyVtAPvCKu+kSoJmqXljjlHXICokx9Z+q8uXa3Tw2ayNLth8ktVEiPxnegcuHdaBZst0pX9dCUUiWqWp/f9vqKyskxoQPVWXupn08NWczs9fvISk+hgsGt+XakZ3plJrsdbyoUZ1CEhfgZy4RkWGqOs89wEnAtzUNaIwxlRERRnRNZUTXVNbn5vPM15t5a2E2r87fztieLbl+VCeGdmpuI73qkUCvSNYAPYDt7qb2wBqgFGf0Vr+QJQwCuyIxJrztzi/glbnbeHneNg4cLaJf2yZcc3InzujbyubzCpFQNG11qGq/qm4LMJsnrJAYExmOHS/h3cXZPPfNFjbvPUJqowQuGdqeS09qT+smDbyOF1FCOvzXHcE1GjgTGBHogbxkhcSYyFJaqnyzcS8vzd3KzLW7iRFhfO+WXDm8I8M6W7NXMAS9j0REuuMUjglAM2AW8B7wi5qGNMaYmoqJEU7pnsYp3dPI2n+UV+Zv482FWUxbuYvuLRtxxfCOnDsgg8ZJ8V5HjQr+7mx/AhgGbAA+Bz5T1Z11lC1o7IrEmMhXUFTCR8t28uLcrazckUeD+FjO6teayUPbMah9M7tKqaZgTpFyMjA3nNYeqYgVEmOih6qyPPsQbyzMYurSHRw5XkLX9EZMPrEd5w1sQ4tGiV5HDAvBLCT/AR5R1ZUV7EsGLgYKVfXVmoatC1ZIjIlORwqL+WR5Dm8s3M7i7QeJjxXG927FxSe2Y2TXVJvbqwrBLCQDgN8AfYGVfL+wVTcgBXgOZ6qUwtqGDiUrJMaY9bn5vLEgi/eWZHPwaBFtmjZg0qA2nDuwDV3SGnkdr94JxfDfRsAQoDXOUrtrVHVdrVLWISskxpgyhcUlfLEql7cys/h2415KFfq3bcJ5A9twVv8MUq3pC7DZf/+HFRJjTEVy8wqYunQn7y/ZweqcPGJjhFO7p3HuwDaM790yqieNDMUVyQqcWX99HQIygf9T1X3VTlmHrJAYY/xZtyuf95fs4MOlO8g5VECjxDhOP6EVE/tnMLxLi6ib2j4UheRBoAR4zd002f2aB4xU1bNrErSuWCExxgSqtFSZt2UfHyzZwbQVu8gvLKZpw3gm9G7Fmf1aMyJKikooCsm3qnpyRdtEZIWq9q1h1jphhcQYUxMFRSXMWb+HT1fkMGPNbg4XFtOkQTzje7fkzH6tOblLKglxkVlUQjH7byMROUlV57sHGAqUDXMorkFGY4yp95LiYxnfpxXj+7SioKiErzfsZdqKHD5buYu3F2WTkhTHab1bccYJrTi5ayoNEqKzTyXQQnId8Jw7ektwmrSude8leaAmBxaRC4F7gV7AUFXN9NnXD3gKZ4hxKXCiqhaIyGDgBaAB8Clwm0bDaAFjjOeS4mM5rXdLTuvdksLiEr7ZsJdPVuTwxepdvLs4m6T4GEZ2TWVcr5aM6ZVOeuMkryPXmWqN2hKRJu57Dtb6wCK9cIrEU8AdZYVEROKAxcAVqrpMRFoAB1W1REQWALcB83AKyaOqOs3fsaxpyxgTKseLS5m/ZR8z1+xm+upcdhw8BkD/dk05rVc6Y3u1pGerxmE3RUso+kiaAPcAp7ibZgN/VNVDNU75/Wd/xQ8LyZnApap6ebnXtQZmqWpP9/klwGhVvcHfMayQGGPqgqqyLjefGatzmb5mN8uynL+52zRtwJie6ZzSPY3hXVrQKDHQxiDvhKKP5DmcO9svcp9fATwPTKp+PL+6AyoinwNpwBuq+iDQBsj2eV22u80YY+oFEaFnqxR6tkrh1jHd2J1XwJdrdzNjTS7vLMrm5XnbiIsRBnVoxindUjmlexonZDQJ+6laAi0kXVT1fJ/n94nIUn9vEpEZQKsKdv1WVT+sItNI4ETgKDBTRBbh9MuUV+nllIhMAaYAtG/f3l9UY4wJuvSUJCYPbc/koe0pLC5h0bYDfL1hL3PW7+HvX6zn71+sp1nDeEZ2S2NElxYM69yCji0ahl0zWKCF5JiIjFTVb+C/swIf8/cmVR1Xg0zZwGxV3ese61NgEPAK0NbndW2BSqe0V9WngafBadqqQQ5jjAmaxLhYRnRJZUSXVO46vSd7DxfyzYa9zNmwh6837OWjZc6vs5YpiZzUySkqJ3VuTufU5HpfWAItJDcBL5Z1tgP7gatClOlz4E4RaQgcB04FHlbVHBHJF5FhwHzgSuBfIcpgjDEhldookXMHOpNGqiqb9hxh/pZ9zNu8n3mb9zHVLSxpjRM5qVNzBrVvxsD2TemdkVLv1qmv7qitFABVraiZqXoHFjkPpxCkAQeBpao6wd13OXA3TtPVp6p6p7t9CN8P/50G/DSQ4b/W2W6MCSeqypa9R5i/xSkqC7bsJ+dQAQAJsTH0aZPCwHZOYRnQriltmzUI+lVLMKeRr3IpXVV9qJrZPGGFxBgT7nYdKmBp1gGWbD/Iku0HWb7jIAVFzpqDLZIT6J2RQu/WKf/92ik1mbhaTOUSzFFbjWucwhhjTNC0apLE6U1ac/oJrQEoKill3a58Fm8/wModh1idk8fz327leIlTXBLjYuidkcI7N44gNsSjwqosJKp6n4gkqmph2deQpjHGGBOQ+NgYTmjThBPaNPnvtqKSUjbtOczqnXms3pnHgaNFIS8iEFhn+5MichPwOHBtiPMYY4ypofjYmP/exzJpUN0dt8oGNBE5FWfNka+BRSJySlWvN8YYE30C7Ymp34OYjTHGeKbKQqKqs3HuMB8FDFHVOXWSyhhjTNgI5IrkRlU9Btwc6jDGGGPCj9/OdncdkDTgehHp6PseVb0mdNGMMcaEg0CnSPkQp8N9Bs7a7cYYYwwQeCFpqKp3hTSJMcaYsBToqK2P3QWnjDHGmB8IdIXEfCAZKASK3c2qqikhzBY0IrIH2Bbkj00F9gb5M8OBnXd0sfOOPmXn3kFV0wJ5Q7Vm/zXfE5HMQCc0iyR23tHFzjv61OTcqz01pIjcW933GGOMiVw1mWN4YtBTGGOMCVs1KSQ2XYrjaa8DeMTOO7rYeUefap97tftIRCRGVUureyBjjDGRyTrbjTHG1ErN12GMUiJyoYisEpFSdw153313i8hGEVknIhO8yhgqInK6e24bReTXXucJFRF5TkR2i8hKn23NRWS6iGxwvzbzMmMoiEg7EZklImvc/8dvc7dH9LmLSJKILBCRZe553+duj+jzLiMisSKyREQ+dp9X+7ytkFTfSmAS8IOZkEWkNzAZ6AOcDjwhIrF1Hy803HN5HDgD6A1c4p5zJHoB57+hr18DM1W1GzDTfR5pioFfqmovYBhwi/vfONLPvRAYo6r9gQHA6SIyjMg/7zK3AWt8nlf7vAMuJCIyxvdrtFLVNaq6roJd5wBvqGqhqm4BNgJD6zZdSA0FNqrqZlU9DryBc84Rx10uYX+5zecAL7rfvwicW6eh6oCq5qjqYvf7fJxfLm2I8HNXx2H3abz7UCL8vAFEpC3wY+AZn83VPu/qXJH8vdxX80NtgCyf59nutkgR6efnT0tVzQHnFy6Q7nGekHJn+h4IzCcKzt1t3lkK7Aamq2pUnDfwT+BOwHcAVbXPO9BJG31F/PBfEZkBtKpg129V9cPK3lbBtkgayRDp52dcItIIeBe4XVXzRCL+nzyqWgIMEJGmwPsicoLXmUJNRM4CdqvqIhEZXZvPqkkhiXiqOq4Gb8sG2vk8bwvsDE6ieiHSz8+fXBFprao5ItIa5y/XiCMi8ThF5FVVfc/dHBXnDqCqB0XkK5w+skg/75OBie6EvElAioi8Qg3O2zrbg2cqMFlEEkWkE9ANWOBxpmBaCHQTkU4ikoAzsGCqx5nq0lTgJ+73P8FZoyeiiHPp8SywRlUf8tkV0ecuImnulQgi0gAYB6wlws9bVe9W1baq2hHn3/OXqno5NThvuyKpJhE5D/gXkAZ8IiJLVXWCqq4SkbeA1TijX25xL5cjgqoWi8itwOdALPCcqq7yOFZIiMjrwGggVUSygXuAvwBvici1wHbgQu8ShszJwBXACre/AOA3RP65twZedEcmxgBvqerHIjKXyD7vylT7v3fANySKyBxVPaXsa+1yGmOMiRR2Z7sxxphasT4SY4wxtWKFxBhjTK1YITHGGFMrARUSEbk1UicsM8YYUzuBXpG0AhaKyFvuDLCRf6urMcaYgARUSFT1dzg32D0LXAVsEJE/i0iXEGYzpkZEpERElrpTgi8TkV+ISEibcUXkRhG5shbvL8ucUcG+Fu707odF5LFy++4XkSwROVz+fZUcZ6h7nKXuz+Y8n32XiMgKEVkuIp+JSKq7/ecisr38sY0pU63hvyLSH7gaZ/qAWThTTU9X1TtDE8+Y6hORw6rayP0+HXgN+FZV7/E2WeV8M1ewLxlnAsUTgBNU9VaffcOAbcCGyt5f7rMaAsfdG0xbA8uAsuK1E+itqntF5EHgqKre677vKmCI77GNKRNoH8nPRGQR8CDwLdBXVW8CBgPnhzCfMbWiqruBKcCt4kgSkefdv7yXiMiPwPlFKSIfiMhHIrLF7Rf8hfuaeSLS3H3d9SKy0P1r/l33FzMicq+I3OF+/5WI/FWcxZLWi8ioWp7DEVX9BiioYN+8splaA/yso6pa7D5N4vuJN8V9JLtN1ylE11xqphYCvdxPBSa5U4G8rapFAO7a7WeFLJ0xQaCqm3H+X08HbnG39QUuwZkaI8l96QnApThrr9yP8xf5QGAuUNZs9Z6qnugugrQGuLaSw8ap6lDgdpwpVhCRDBH5NNjnV10icpKIrAJWADeqarH7b/omd9tOnMXLnvUwpgkjgfaR/EFVt1Wyb01F242pZ8oGiIwEXgZQ1bU4zULd3X2zVDVfVfcAh4CP3O0rgI7u9yeIyNcisgK4DGdFzIqUzZy7qOy9qrpTVc8MytnUgqrOV9U+wInA3e5VWjxOIRmI09S1HLjbw5gmjNh9JCbiiUhnoARnOuyqRhwW+nxf6vO8lO8nOH0BuNW9orkPp3moqs8qoZqTo4rIeT4d4kOq897qcP8IPIJzJTbA3bZJnY7Tt4ARoTq2iSxWSExEE5E04EngMfcX5BycKwlEpDvQHqho6eTKNAZy3L/gLwtyXABU9X1VHeA+MmvyGW4xeqCC7Z1EJM79vgPQA9gK7AB6uz8vgNP44TrexlTKppE3kaiBOw16PM6U/i8DZetrPAE86TZNFQNXqWphNW6N+j3O8rPbcJq8Ggf6Rndo7zM1ad4Ska04HeAJInIuMF5VV7ujqy4FGrpT3j/jjrTqAuRV8FEjgV+LSBHOldbNqrrXPcZ9wBx33zacof7G+GWz/xpTD1Q1/LeGn/cK8HO3vycYn3cVNvzXVMKatoypH/IquyGxJlT18iAWkZ/jdLxXdIVjjF2RGGOMqR27IjHGGFMrVkiMMcbUihUSY4wxtWKFxBhjTK1YITHGGFMr/w8njBQdhhvzHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now use Newton-Raphson method to find the maximum\n",
      "Starting point -11; max -0.10132118364233782; iterations 50\n",
      "Starting point -11; max -0.10132118364233782; iterations 50\n",
      "No maximum found in [0, 1]\n",
      "No maximum found in [1.5, 3]\n",
      "No maximum found in [4, 5]\n",
      "No maximum found in [4.7, 6]\n",
      "No maximum found in [7, 9]\n",
      "No maximum found in [8, 10]\n",
      "No maximum found in [38, 40]\n"
     ]
    }
   ],
   "source": [
    "logLikely = lambda x, th: np.log((np.pi * (1 + (x - th)**2))**-1).sum()\n",
    "obs = np.array([1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71, -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21])\n",
    "\n",
    "thetas = np.linspace(-11, 38, 200)\n",
    "y = np.array([logLikely(obs, t) for t in thetas])\n",
    "plt.plot(thetas, y)\n",
    "plt.title(\"Log likelihood function\")\n",
    "plt.ylabel(\"y = -n*log(\\pi) - \\sum_{i=1}^{n}log(1+(x-\\\\theta)^2)\")\n",
    "plt.xlabel(\"Domain: [-11, 38]\")\n",
    "plt.show()\n",
    "\n",
    "# now, find estimates for theta using the Fixed-Point Iterative Method\n",
    "gp = lambda x, th: 2 * np.sum((x - th)/(np.pi*(1 + (x - th)**2))**2)\n",
    "gpp = lambda x, th: np.sum((2*th - 2*x)/(1 + (x - th)**2)**2)\n",
    "step = lambda x, th: gp(x, th)/gpp(x, th)\n",
    "\n",
    "print(\"Now use Newton-Raphson method to find the maximum\")\n",
    "allIntervals = [[-11, 0.5], [-1, 0.5], [0, 1], [1.5, 3], [4, 5], [4.7, 6], [7, 9], [8, 10], [38, 40]]\n",
    "for interval in allIntervals:\n",
    "    try:\n",
    "        vals = newtonsMethod(interval, interval[0], gp, step, maxIterations = 50, samples = obs)\n",
    "        print(\"Starting point {}; max {}; iterations {}\".format(-11, vals[0][0], vals[0][1]))\n",
    "    except:\n",
    "        print(\"No maximum found in [{}, {}]\".format(interval[0], interval[1]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2.1(b): Find MLE using bisection method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum -1.5533092021942139; Iterations 22; Initial Guess -1\n",
      "Maximum 1.2570688724517822; Iterations 23; Initial Guess 1\n"
     ]
    }
   ],
   "source": [
    "for intv in [[-2, 0, -1], [0, 3, 1]]:\n",
    "    vals = bisectionMethod(intv[0], intv[1], gp, maxIterations = 50, guess = intv[2], samples = obs)\n",
    "    print(\"Maximum {}; Iterations {}; Initial Guess {}\".format(vals[0], vals[1], intv[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2.1(c): Find MLE using Fixed-Point Iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged to -1.5531506088155014 after 64 iterations: L 1\n",
      "Converged to -1.5531506088155014 after 64 iterations: L 0.64\n",
      "Converged to -1.5531506088155014 after 64 iterations: L 0.25\n"
     ]
    }
   ],
   "source": [
    "vals = fixedPointMethod(-3, 0, gp, 1, samples = obs)\n",
    "print(\"Converged to {} after {} iterations: L {}\".format(vals[0], vals[1], 1))\n",
    "vals = fixedPointMethod(-3, 0, gp, 0.64, samples = obs)\n",
    "print(\"Converged to {} after {} iterations: L {}\".format(vals[0], vals[1], 0.64))\n",
    "vals = fixedPointMethod(-3, 0, gp, 0.25, samples = obs)\n",
    "print(\"Converged to {} after {} iterations: L {}\".format(vals[0], vals[1], 0.25))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2.1(d): Find MLE using Secant Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x0 -1; x1 -1.5510369668105262; test -0.04440288251240235\n",
      "x0 -1.5510369668105262; x1 -1.553446710253367; test -0.00019333307614838758\n",
      "x0 -1.553446710253367; x1 -1.5533094460158816; test 1.1677870237777285e-05\n",
      "x0 -1.5533094460158816; x1 -1.5533091724900308; test 2.3224162697631822e-08\n",
      "Theta -1.5533091724900308 in 1 iterations\n",
      "x0 3; x1 13.809263334552151; test 0.1929201298543256\n",
      "x0 13.809263334552151; x1 13.697351148454327; test -0.002018267249027645\n",
      "x0 13.697351148454327; x1 17.453306326849106; test -0.0020802502444110097\n",
      "x0 17.453306326849106; x1 20.138399421037573; test -0.0008671982151147096\n",
      "x0 20.138399421037573; x1 24.21136082765596; test -0.0005226452171447957\n",
      "x0 24.21136082765596; x1 28.275009267885945; test -0.00026102350504943536\n",
      "x0 28.275009267885945; x1 31.27537676991111; test -0.00011086702637947404\n",
      "x0 31.27537676991111; x1 31.299590107232294; test -8.875480108189952e-07\n",
      "x0 31.299590107232294; x1 31.2955491982033; test 1.7779221810250174e-07\n",
      "Theta 31.2955491982033 in 1 iterations\n"
     ]
    }
   ],
   "source": [
    "vals = secantMethod(-2, -1, gp, maxIterations = 50, samples = obs, debug = True)\n",
    "print(\"Theta {} in {} iterations\".format(vals[0], vals[1]))\n",
    "vals = secantMethod(-3, 3, gp, maxIterations = 50, samples = obs, debug = True)\n",
    "print(\"Theta {} in {} iterations\".format(vals[0], vals[1]))"
   ]
  },
  {
   "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}