EstimatingOpt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating Optimal Solutions for Combinatorial Optimization Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $S$ be a sample of size $m$ from a probability distribution bounded below by $a$, let $z$ be the minimum value in the sample. As $m$ increases, the distribution of $z$ approaches a *Weibull distribution* having $a$ as its location parameter (displacement) (Golden and Alt, 1979; Rardin and Uzsoy, 2001).\n",
    "\n",
    "The CDF of the Weibull distribution is $F_x(x_0) = prob\\{x \\leq x_0\\} = 1 - e^{-(\\frac{x_0-a}{b})^c}$, where $a$ is the location parameter, $b$ is the scale parameter, and $c$ the shape parameter.\n",
    "\n",
    "By considering $k$ samples and their corresponding minimum values $z_1, \\ldots, z_k$, we can try to estimate the parameters of the Weibull distribution. We will assume $z_1, \\ldots, z_k$ are in non-decreasing order.\n",
    "\n",
    "We can use this result to estimate the value of the optimal solution of a combinatorial optimization problem.\n",
    "\n",
    "By taking a problem instance as a probability distribution, the best solution in a sample $m$ random solutions for the instance should approach a Weibull distribution. By considering $k$ samples, we can estimate the distribution parameters, in particular, $a$ is an estimate of the optimal solution value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An example - The Traveling Salesman Problem\n",
    "Sample $n$ points on a 2D plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "n = 10\n",
    "side = 100\n",
    "x = [random.random()*side for i in range(n)]\n",
    "y = [random.random()*side for i in range(n)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEyxJREFUeJzt3X9sXWd9x/H3d04KLpvm/jBR47RrpkZmCNQGrKqIH2L9MbcFUQtVFQiJbKqWf9AoMBka7Q+ENI1WQSsgTZUiCoSJlZZikqqaGrpQxl8EHNw1pcU0lJbmJm0M1MDAGmn47o97XJxg1/E9176+j98vybrnPPece78nT/LJ8XOee09kJpKkcv1JpwuQJC0vg16SCmfQS1LhDHpJKpxBL0mFM+glqXCLBn1EfD4ijkfEY3Pazo2IhyLiyerxnKo9IuKzEXE4Ih6NiDcsZ/GSpMWdyRn9F4FrT2u7FdifmVuA/dU6wHXAlupnO3Bne8qUJLVq0aDPzG8Dvzit+QZgd7W8GxiZ0/6lbPoO0BcRF7SrWEnS0q1rcb8NmXmsWn4O2FAtDwDPztnuSNV2jJdx/vnn58UXX9xiKZK0Nh08ePBnmdm/2HatBv1LMjMjYsnfoxAR22kO73DRRRcxPj5etxRJWlMi4pkz2a7VWTfPzw7JVI/Hq/YGcOGc7TZVbX8kM3dl5lBmDvX3L/ofkiSpRa0G/f3Atmp5G7B3Tvv7q9k3VwC/nDPEI0nqgEWHbiLibuDtwPkRcQT4OHAbcG9E3Aw8A9xUbf6fwPXAYeC3wN8tQ82SpCVYNOgz870LPHXVPNsm8IG6RUmS2sdPxkpS4WrPuumUPRMNdu6b5Oj0DBv7ehkdHmRk60Cny5KkVacrg37PRIMdY4eYOXESgMb0DDvGDgEY9pJ0mq4cutm5b/KlkJ81c+IkO/dNdqgiSVq9ujLoj07PLKldktayrgz6jX29S2qXpLWsK4N+dHiQ3vU9p7T1ru9hdHiwQxVJ0urVlRdjZy+4OutGkhbXlUEPzbA32CVpcV05dCNJOnMGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwXfvJWEnqpG66+VGtM/qIuCUiHouIH0TEh6q2cyPioYh4sno8pz2lStLqMHvzo8b0DMkfbn60Z6LR6dLm1XLQR8TrgL8HLgcuBd4ZEZcAtwL7M3MLsL9al6RidNvNj+qc0f8VcCAzf5uZLwL/DbwbuAHYXW2zGxipV6IkrS7ddvOjOkH/GPDWiDgvIs4GrgcuBDZk5rFqm+eADTVrlKRVpdtuftRy0GfmE8DtwDeAB4FHgJOnbZNAzrd/RGyPiPGIGJ+ammq1DElacd1286NaF2Mz867MfGNmvg14AfgR8HxEXABQPR5fYN9dmTmUmUP9/f11ypCkFTWydYBPvvv1DPT1EsBAXy+ffPfrV+2sm1rTKyPi1Zl5PCIuojk+fwWwGdgG3FY97q1dpSStMt1086O68+i/FhHnASeAD2TmdETcBtwbETcDzwA31S1SktS6WkGfmW+dp+3nwFV1XleS1D5+BYIkFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUuFpBHxEfjogfRMRjEXF3RLwyIjZHxIGIOBwR90TEWe0qVpK0dC0HfUQMAB8EhjLzdUAP8B7gduCOzLwEeAG4uR2FSpJaU3foZh3QGxHrgLOBY8CVwH3V87uBkZrvIUmqoeWgz8wG8CngpzQD/pfAQWA6M1+sNjsCDMy3f0Rsj4jxiBifmppqtQxJ0iLWtbpjRJwD3ABsBqaBrwLXnun+mbkL2AUwNDSUrdYhaWn2TDTYuW+So9MzbOzrZXR4kJGt856PqRAtBz1wNfCTzJwCiIgx4M1AX0Ssq87qNwGN+mVKaoc9Ew12jB1i5sRJABrTM+wYOwRg2Beszhj9T4ErIuLsiAjgKuBx4GHgxmqbbcDeeiVKaped+yZfCvlZMydOsnPfZIcq0kqoM0Z/gOZF1+8Dh6rX2gV8DPhIRBwGzgPuakOdktrg6PTMktpVhjpDN2Tmx4GPn9b8FHB5ndeVtDw29vXSmCfUN/b1dqAarRQ/GSutIaPDg/Su7zmlrXd9D6PDgx2qSCuh1hm9pO4ye8HVWTdri0EvrTEjWwcM9jXGoRtJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIK13LQR8RgRDwy5+dXEfGhiDg3Ih6KiCerx3PaWbAkaWnq3Bx8MjMvy8zLgDcCvwW+DtwK7M/MLcD+al2S1CHtGrq5CvhxZj4D3ADsrtp3AyNteg9JUgvadSvB9wB3V8sbMvNYtfwcsKFN7yGtOXsmGt7fVbXVPqOPiLOAdwFfPf25zEwgF9hve0SMR8T41NRU3TKk4uyZaLBj7BCN6RkSaEzPsGPsEHsmGp0uTV2mHUM31wHfz8znq/XnI+ICgOrx+Hw7ZeauzBzKzKH+/v42lCGVZee+SWZOnDylbebESXbum+xQRepW7Qj69/KHYRuA+4Ft1fI2YG8b3kNac45OzyypXVpIraCPiFcB1wBjc5pvA66JiCeBq6t1SUu0sa93Se3SQmpdjM3M3wDnndb2c5qzcCTVMDo8yI6xQ6cM3/Su72F0eHDBfbx4q/m0a9aNpDabDegzDe7Zi7ez/zHMXryd+1pamwx6aRUb2TpwxiH9chdvDfq1ze+6kQrhxVstxKCXCuHFWy3EoJcKMTo8SO/6nlPaFrt4q7XBMXqpEEu9eKu1w6CXCrKUi7daOxy6kaTCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1Lh6t4cvC8i7ouIH0bEExHxpog4NyIeiognq8dz2lWsJGnp6p7RfwZ4MDNfA1wKPAHcCuzPzC3A/mpdktQhLX9NcUT8OfA24G8BMvN3wO8i4gbg7dVmu4FvAR+rU+RK2jPR8Pu8JRWlzvfRbwamgC9ExKXAQeAWYENmHqu2eQ7YUK/ElbNnosGOsUMv3WC5MT3DjrFDAIa9tIp5gvby6gzdrAPeANyZmVuB33DaME1mJpDz7RwR2yNiPCLGp6amapTRPjv3Tb4U8rNmTpxk577JDlUkaTGzJ2iN6RmSP5yg7ZlodLq0VaNO0B8BjmTmgWr9PprB/3xEXABQPR6fb+fM3JWZQ5k51N/fX6OM9jk6PbOkdkmd5wna4loO+sx8Dng2ImbvPHwV8DhwP7CtatsG7K1V4Qra2Ne7pHZJnecJ2uLqzrr5B+DLEfEocBnwL8BtwDUR8SRwdbXeFUaHB+ld33NKW+/6HkaHBxfYQ1KneYK2uFo3B8/MR4CheZ66qs7rdsrsxRsv6kjdY3R48JRJFOAJ2ulqBX2JRrYOGOxSF/EEbXEGvaSu5wnay/O7biSpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwte4wFRFPA78GTgIvZuZQRJwL3ANcDDwN3JSZL9QrU5LUqnac0f91Zl6WmbM3Cb8V2J+ZW4D91bokqUOWY+jmBmB3tbwbGFmG95AknaG6QZ/ANyLiYERsr9o2ZOaxavk5YMN8O0bE9ogYj4jxqampmmVIkhZSa4weeEtmNiLi1cBDEfHDuU9mZkZEzrdjZu4CdgEMDQ3Nu40kqb5aZ/SZ2agejwNfBy4Hno+ICwCqx+N1i5Qkta7loI+IV0XEn80uA38DPAbcD2yrNtsG7K1bpCSpdXWGbjYAX4+I2df5j8x8MCK+B9wbETcDzwA31S9TktSqloM+M58CLp2n/efAVXWKkiS1j5+MlaTCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYWr+xUIktRWeyYa7Nw3ydHpGTb29TI6PMjI1oFOl9XVDHpJq8aeiQY7xg4xc+IkAI3pGXaMHQIw7Gtw6EbSqrFz3+RLIT9r5sRJdu6b7FBFZTDoJa0aR6dnltSuM2PQS1o1Nvb1LqldZ8agl7RqjA4P0ru+55S23vU9jA4PdqiiMngxVtKqMXvB1Vk37WXQS1pVRrYOGOxt5tCNJBXOoJekwhn0klQ4g16SClc76COiJyImIuKBan1zRByIiMMRcU9EnFW/TElSq9pxRn8L8MSc9duBOzLzEuAF4OY2vIckqUW1gj4iNgHvAD5XrQdwJXBftcluYKTOe0iS6ql7Rv9p4KPA76v184DpzHyxWj8COCFWkjqo5aCPiHcCxzPzYIv7b4+I8YgYn5qaarUMSdIi6pzRvxl4V0Q8DXyF5pDNZ4C+iJj9xO0moDHfzpm5KzOHMnOov7+/RhmSpJfTctBn5o7M3JSZFwPvAb6Zme8DHgZurDbbBuytXaUkqWXLMY/+Y8BHIuIwzTH7u5bhPSRJZ6gtX2qWmd8CvlUtPwVc3o7XlSTV5ydjJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMK1HPQR8cqI+G5E/E9E/CAiPlG1b46IAxFxOCLuiYiz2leuJGmp6pzR/x9wZWZeClwGXBsRVwC3A3dk5iXAC8DN9cuUJLWq5aDPpv+tVtdXPwlcCdxXte8GRmpVKEmqpdYYfUT0RMQjwHHgIeDHwHRmvlhtcgQYqFeiJKmOWkGfmScz8zJgE3A58Joz3TcitkfEeESMT01N1SlDkvQy2jLrJjOngYeBNwF9EbGuemoT0Fhgn12ZOZSZQ/39/e0oQ5I0jzqzbvojoq9a7gWuAZ6gGfg3VpttA/bWLVKS1Lp1i2+yoAuA3RHRQ/M/jHsz84GIeBz4SkT8MzAB3NWGOiVJLWo56DPzUWDrPO1P0RyvlzSPPRMNdu6b5Oj0DBv7ehkdHmRkq3MWtHzqnNFLWqI9Ew12jB1i5sRJABrTM+wYOwRg2GvZ+BUI0grauW/ypZCfNXPiJDv3TXaoIq0FBr20go5OzyypXWoHg15aQRv7epfULrWDQS+toNHhQXrX95zS1ru+h9HhwQ5VpLXAi7HSCpq94OqsG63k7CuDXlphI1sHDPY1bqVnXzl0I0krbKVnXxn0krTCVnr2lUEvSStspWdfGfSStMJWevaVF2MlaYWt9Owrg16SOmAlZ185dCNJhTPoJalwBr0kFc6gl6TCGfSSVLjIzE7XQERMAc8scbfzgZ8tQzmrnce99qzVY/e4F/cXmdm/2EarIuhbERHjmTnU6TpWmse99qzVY/e428ehG0kqnEEvSYXr5qDf1ekCOsTjXnvW6rF73G3StWP0kqQz081n9JKkM9CVQR8R10bEZEQcjohbO13PcomICyPi4Yh4PCJ+EBG3VO3nRsRDEfFk9XhOp2tdDhHRExETEfFAtb45Ig5U/X5PRJzV6RrbLSL6IuK+iPhhRDwREW9aC/0dER+u/o4/FhF3R8QrS+3viPh8RByPiMfmtM3bx9H02erP4NGIeEMr79l1QR8RPcC/AdcBrwXeGxGv7WxVy+ZF4B8z87XAFcAHqmO9FdifmVuA/dV6iW4BnpizfjtwR2ZeArwA3NyRqpbXZ4AHM/M1wKU0j7/o/o6IAeCDwFBmvg7oAd5Duf39ReDa09oW6uPrgC3Vz3bgzlbesOuCHrgcOJyZT2Xm74CvADd0uKZlkZnHMvP71fKvaf6jH6B5vLurzXYDI52pcPlExCbgHcDnqvUArgTuqzYp7rgj4s+BtwF3AWTm7zJzmjXQ3zS/Mr03ItYBZwPHKLS/M/PbwC9Oa16oj28AvpRN3wH6IuKCpb5nNwb9APDsnPUjVVvRIuJiYCtwANiQmceqp54DNnSorOX0aeCjwO+r9fOA6cx8sVovsd83A1PAF6ohq89FxKsovL8zswF8CvgpzYD/JXCQ8vt7roX6uC15141Bv+ZExJ8CXwM+lJm/mvtcNqdNFTV1KiLeCRzPzIOdrmWFrQPeANyZmVuB33DaME2h/X0OzTPXzcBG4FX88dDGmrEcfdyNQd8ALpyzvqlqK1JErKcZ8l/OzLGq+fnZX9+qx+Odqm+ZvBl4V0Q8TXNo7kqaY9d91a/2UGa/HwGOZOaBav0+msFfen9fDfwkM6cy8wQwRvPvQOn9PddCfdyWvOvGoP8esKW6In8WzYs293e4pmVRjUvfBTyRmf8656n7gW3V8jZg70rXtpwyc0dmbsrMi2n27zcz833Aw8CN1WYlHvdzwLMRMXuH6KuAxym8v2kO2VwREWdXf+dnj7vo/j7NQn18P/D+avbNFcAv5wzxnLnM7Lof4HrgR8CPgX/qdD3LeJxvofkr3KPAI9XP9TTHq/cDTwL/BZzb6VqX8c/g7cAD1fJfAt8FDgNfBV7R6fqW4XgvA8arPt8DnLMW+hv4BPBD4DHg34FXlNrfwN00r0WcoPlb3M0L9TEQNGcZ/hg4RHNm0pLf00/GSlLhunHoRpK0BAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mF+388mNC7T6ucxwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.scatter(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions for the distance between two vertices and cost of a solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def distance(u, v):\n",
    "    return ((x[u]-x[v])**2+(y[u]-y[v])**2)**(1/2)\n",
    "\n",
    "def cost(s):\n",
    "    return sum(distance(s[u], s[(u+1)%n]) for u in range(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for generating a random solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "def randSol():\n",
    "    s = list(range(n))\n",
    "    random.shuffle(s)\n",
    "    \n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for obtaining the best solution from a sample of $m$ random solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimum(m):\n",
    "    return min(cost(randSol()) for i in range(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for obtaining the best solutions of $k$ samples of $m$ solutions each:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minima(k, m):\n",
    "    return sorted([minimum(m) for i in range(k)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets plot a histogram of $z_1, \\ldots, z_k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEK9JREFUeJzt3X2MZXV9x/H3R0C01QrKlOIudomssdiW1UwIrU2rkFYemi4mSiFWtxa7psFUrUmLpokP0QRTlZY+0KyFulgUiGLZKLVFJLEkFZxFRGBFp7qE3SK7yoMSIsnCt3/Mb/GyzMOdp70zP9+v5GbO+Z3fvff7zZl87pkz596bqkKS1K9njLoASdLyMuglqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTt01AUAHHXUUbVu3bpRlyFJq8r27dt/UFVjc81bEUG/bt06JiYmRl2GJK0qSe4ZZp6nbiSpcwa9JHVuzqBP8qwktyT5RpI7k7y/jR+X5OYkk0muSvLMNn54W59s29ctbwuSpNkMc0T/GHBKVZ0IbABOS3Iy8GHgoqo6HngQOK/NPw94sI1f1OZJkkZkzqCvKY+01cParYBTgM+08a3AWW15Y1unbT81SZasYknSvAx1jj7JIUluA/YA1wP/CzxUVfvalF3Amra8BrgXoG1/GHjBUhYtSRreUEFfVY9X1QZgLXAS8NLFPnGSzUkmkkzs3bt3sQ8nSZrBvK66qaqHgBuB3wCOSLL/Ovy1wO62vBs4FqBtfx7ww2kea0tVjVfV+NjYnNf7S5IWaJirbsaSHNGWnw38LrCDqcB/XZu2Cbi2LW9r67TtXy6/mFaSRmaYd8YeA2xNcghTLwxXV9Xnk9wFXJnkg8DXgUvb/EuBTyaZBB4AzlmGuqWDZt0FXxjJ8+688MyRPK/6M2fQV9XtwMunGf8uU+frDxz/CfD6JalOkrRovjNWkjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktS5Yb5hStII+M1WWioe0UtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUufmDPokxya5McldSe5M8vY2/r4ku5Pc1m5nDNzn3Ukmk9yd5DXL2YAkaXbDfATCPuBdVXVrkucC25Nc37ZdVFUfGZyc5ATgHOBlwAuBLyV5SVU9vpSFS5KGM+cRfVXdV1W3tuUfAzuANbPcZSNwZVU9VlXfAyaBk5aiWEnS/M3rHH2SdcDLgZvb0NuS3J7ksiRHtrE1wL0Dd9vF7C8MkqRlNHTQJ3kO8FngHVX1I+AS4MXABuA+4KPzeeIkm5NMJJnYu3fvfO4qSZqHoYI+yWFMhfwVVXUNQFXdX1WPV9UTwMf56emZ3cCxA3df28aeoqq2VNV4VY2PjY0tpgdJ0iyGueomwKXAjqr62MD4MQPTXgvc0Za3AeckOTzJccB64JalK1mSNB/DXHXzSuCNwDeT3NbG3gOcm2QDUMBO4K0AVXVnkquBu5i6Yud8r7iRpNGZM+ir6iYg02y6bpb7fAj40CLqkiQtEd8ZK0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXNzBn2SY5PcmOSuJHcmeXsbf36S65N8p/08so0nycVJJpPcnuQVy92EJGlmwxzR7wPeVVUnACcD5yc5AbgAuKGq1gM3tHWA04H17bYZuGTJq5YkDW3OoK+q+6rq1rb8Y2AHsAbYCGxt07YCZ7XljcDlNeWrwBFJjlnyyiVJQ5nXOfok64CXAzcDR1fVfW3T94Gj2/Ia4N6Bu+1qY5KkERg66JM8B/gs8I6q+tHgtqoqoObzxEk2J5lIMrF379753FWSNA9DBX2Sw5gK+Suq6po2fP/+UzLt5542vhs4duDua9vYU1TVlqoar6rxsbGxhdYvSZrDMFfdBLgU2FFVHxvYtA3Y1JY3AdcOjL+pXX1zMvDwwCkeSdJBdugQc14JvBH4ZpLb2th7gAuBq5OcB9wDnN22XQecAUwCjwJvXtKKJUnzMmfQV9VNQGbYfOo08ws4f5F1SZKWiO+MlaTOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1Lk5gz7JZUn2JLljYOx9SXYnua3dzhjY9u4kk0nuTvKa5SpckjScQ4eY8wngH4DLDxi/qKo+MjiQ5ATgHOBlwAuBLyV5SVU9vgS16mfYugu+MOoSpFVrziP6qvoK8MCQj7cRuLKqHquq7wGTwEmLqE+StEiLOUf/tiS3t1M7R7axNcC9A3N2tTFJ0ogsNOgvAV4MbADuAz463wdIsjnJRJKJvXv3LrAMSdJcFhT0VXV/VT1eVU8AH+enp2d2A8cOTF3bxqZ7jC1VNV5V42NjYwspQ5I0hAUFfZJjBlZfC+y/ImcbcE6Sw5McB6wHbllciZKkxZjzqpsknwZeBRyVZBfwXuBVSTYABewE3gpQVXcmuRq4C9gHnO8VN5I0WqmqUdfA+Ph4TUxMjLoMrWBeXvmzYeeFZ466hFUlyfaqGp9rnu+MlaTOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1Lk5gz7JZUn2JLljYOz5Sa5P8p3288g2niQXJ5lMcnuSVyxn8ZKkuQ1zRP8J4LQDxi4Abqiq9cANbR3gdGB9u20GLlmaMiVJCzVn0FfVV4AHDhjeCGxty1uBswbGL68pXwWOSHLMUhUrSZq/hZ6jP7qq7mvL3weObstrgHsH5u1qY5KkEVn0P2OrqoCa7/2SbE4ykWRi7969iy1DkjSDhQb9/ftPybSfe9r4buDYgXlr29jTVNWWqhqvqvGxsbEFliFJmstCg34bsKktbwKuHRh/U7v65mTg4YFTPJKkETh0rglJPg28CjgqyS7gvcCFwNVJzgPuAc5u068DzgAmgUeBNy9DzZKkeZgz6Kvq3Bk2nTrN3ALOX2xRkqSl4ztjJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5w4ddQFaXdZd8IVRlyBpnjyil6TOGfSS1DmDXpI6t6hz9El2Aj8GHgf2VdV4kucDVwHrgJ3A2VX14OLKlCQt1FIc0b+6qjZU1XhbvwC4oarWAze0dUnSiCzHqZuNwNa2vBU4axmeQ5I0pMUGfQH/lWR7ks1t7Oiquq8tfx84epHPIUlahMVeR/9bVbU7yS8C1yf51uDGqqokNd0d2wvDZoAXvehFiyxDkjSTRR3RV9Xu9nMP8DngJOD+JMcAtJ97Zrjvlqoar6rxsbGxxZQhSZrFgoM+yc8nee7+ZeD3gDuAbcCmNm0TcO1ii5QkLdxiTt0cDXwuyf7H+VRVfTHJ14Crk5wH3AOcvfgyJUkLteCgr6rvAidOM/5D4NTFFCVJWjq+M1aSOmfQS1LnDHpJ6pxBL0md84tHJK0Yo/pim50XnjmS5z1YPKKXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzXke/Co3qWmNJq5NH9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnB9qJuln3ig/KPBgfDH5sh3RJzktyd1JJpNcsFzPI0ma3bIEfZJDgH8ETgdOAM5NcsJyPJckaXbLdermJGCyqr4LkORKYCNw11I/kZ/NLkmzW65TN2uAewfWd7UxSdJBNrJ/xibZDGxuq48kuXuGqUcBPzg4VR009rQ69NZTb/1ABz3lw08bmk9PvzzMpOUK+t3AsQPra9vYk6pqC7BlrgdKMlFV40tb3mjZ0+rQW0+99QP2NKzlOnXzNWB9kuOSPBM4B9i2TM8lSZrFshzRV9W+JG8D/hM4BLisqu5cjueSJM1u2c7RV9V1wHVL8FBznt5Zhexpdeitp976AXsaSqpqqR9TkrSC+Fk3ktS5kQd9kmcluSXJN5LcmeT9bfyK9hEKdyS5LMlhbTxJLm4frXB7kleMtoOnmqmfge0XJ3lkYP3wJFe1fm5Osu5g1zyXWfZRknwoybeT7Ejy5wPjK3Yfwaw9nZrk1iS3JbkpyfFtfMXvJ5h6V3qSryf5fFs/rtU72ep/ZhtfFf3AtD2tymwYdGBPA+PLkg8jD3rgMeCUqjoR2ACcluRk4ArgpcCvAc8G3tLmnw6sb7fNwCUHveLZzdQPScaBIw+Yfx7wYFUdD1wEPP2q2tGbqac/Zuoy2pdW1a8AV7b5K30fwcw9XQK8oao2AJ8C/rrNXw37CeDtwI6B9Q8DF7W6H2SqD1g9/cDTe1qt2TDowJ6WNR9GHvQ1Zf8r2GHtVlV1XdtWwC1MXYsPUx+lcHnb9FXgiCTHHPzKpzdTP5n6/J+/Af7ygLtsBLa25c8ApybJQSl2SDP1BPwZ8IGqeqLN29PmrOh9BLP2VMAvtPHnAf/Xllf8fkqyFjgT+Je2HuAUpuqFqfrPassrvh94ek8wdaHHasyG/abrabnzYeRBD0/+GXMbsAe4vqpuHth2GPBG4IttaMV/vMIM/bwN2FZV9x0w/cl+qmof8DDwgoNZ7zBm6OnFwB8mmUjyH0nWt+krfh/BjD29BbguyS6mfu8ubNNXw376W6aC4om2/gLgoVYvPHU/rIZ+4Ok9PWk1ZkMzXU/Lmg8rIuir6vH2p/Ja4KQkvzqw+Z+Ar1TVf4+muvmbpp/fBl4P/P1oK1u4GfbR4cBP2rv4Pg5cNsoa52uGnt4JnFFVa4F/BT42yhqHleT3gT1VtX3UtSyVIXpaddkwXU9JXsgy58OKCPr9quoh4EbgNIAk7wXGgL8YmDbnxyusFAP9vBo4HphMshP4uSSTbdqT/SQ5lKnTBT88+NUO54B9tAu4pm36HPDrbXnV7CN4Sk+nAycO/EV5FfCbbXml76dXAn/Qfr+uZOqUzd8xdfpi//tlBvfDSu8Hpukpyb/Bqs6G6fbTnSxzPow86JOMJTmiLT8b+F3gW0neArwGOHf/OeBmG/Cm9h/2k4GHp/lzZ2Rm6Gd7Vf1SVa2rqnXAo+2fKzDVz6a2/Drgy7XC3tww0z4C/p2pFzGA3wG+3ZZX9D6CGXvaATwvyUvatP1jsML3U1W9u6rWtt+vc5iq7w1MvYC9rk3bBFzblld0PzBjT3+0WrMBZuzpyOXOh5XwVYLHAFvbPyOeAVxdVZ9Psg+4B/if9r+Ha6rqA0y92/YMYBJ4FHjzaMqe0bT9zDL/UuCT7RX8AaZ2/koz0z66CbgiyTuBR/jp1Q8rfR/BzD39KfDZJE8wdZXKn7T5q2E/TeevgCuTfBD4OlN9wOrtB+CfWZ3ZsBBLsp98Z6wkdW7kp24kScvLoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXP/D0JDejj8oFtPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "323.0778898792385"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = 1000\n",
    "m = 100\n",
    "z = minima(k, m)\n",
    "plt.hist(z)\n",
    "plt.show()\n",
    "z[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating the Weibull Distribution Parameters\n",
    "\n",
    "A point estimate for the optimal value $a$ given by Golden and Alt (1979) is $\\tilde{a} = 2z_1 - z_2$. An estimate given by Zanakis (1979) is $\\hat{a} = \\frac{z_1z_n-z_2^2}{z_1+z_n-2z_2}$. Cooke (1979) proposes $a^* = 2z_1-(e-1)\\sum_{i=1}^{k}\\frac{z_i}{e^i}$.\n",
    "\n",
    "According to Zanakis and Evans (1981) the order of accuracy is $a^*$, $\\hat{a}$, $\\tilde{a}$.\n",
    "\n",
    "A $1-e^{-k}$% confidence interval for $a$ was proposed by Golden and Alt (1979): $[z_1-\\tilde{b}, z_1]$, where $\\tilde{b} = z_{\\lfloor 0.63k \\rfloor+1} - \\overline{a}$; $\\overline{a}$ can be $\\tilde{a}, \\hat{a}, a^*$.\n",
    "\n",
    "Golden and Alt (1979) also give an estimate for $c$, $\\tilde{c} = \\frac{\\ln(-\\ln(0.5))}{\\ln(z_M-\\overline{a})-\\ln{\\overline{b}}}$, where $z_M$ is the median of $z_1, \\ldots, z_k$.\n",
    "\n",
    "In fact, Golden and Alt (1979) recommend using Maximum Likelihood Estimation (MLE) to refine $\\overline{a}$, $\\tilde{b}$, $\\tilde{c}$. According to Zanakis and Evans (1981), this is not justifiable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>c</th>\n",
       "      <th>linf</th>\n",
       "      <th>lsup</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>Goden and Alt (1979)</td>\n",
       "      <td>318.939104</td>\n",
       "      <td>76.903089</td>\n",
       "      <td>5.763308</td>\n",
       "      <td>246.174801</td>\n",
       "      <td>323.07789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Zanakis (1979)</td>\n",
       "      <td>322.917944</td>\n",
       "      <td>72.924249</td>\n",
       "      <td>5.455432</td>\n",
       "      <td>250.153641</td>\n",
       "      <td>323.07789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Cooke (1979)</td>\n",
       "      <td>320.473675</td>\n",
       "      <td>75.368517</td>\n",
       "      <td>5.644567</td>\n",
       "      <td>247.709373</td>\n",
       "      <td>323.07789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>MLE</td>\n",
       "      <td>305.617009</td>\n",
       "      <td>90.884685</td>\n",
       "      <td>6.219493</td>\n",
       "      <td>232.193205</td>\n",
       "      <td>323.07789</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                               a          b         c        linf       lsup\n",
       "Goden and Alt (1979)  318.939104  76.903089  5.763308  246.174801  323.07789\n",
       "Zanakis (1979)        322.917944  72.924249  5.455432  250.153641  323.07789\n",
       "Cooke (1979)          320.473675  75.368517  5.644567  247.709373  323.07789\n",
       "MLE                   305.617009  90.884685  6.219493  232.193205  323.07789"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import weibull_min\n",
    "from math import log as ln\n",
    "from statistics import median\n",
    "import pandas as pd\n",
    "\n",
    "def estimates(a):\n",
    "    b = z[math.floor(0.63*k)+1]-a\n",
    "    c = ln(-ln(0.5))/(ln(median(z)-a)-ln(b))\n",
    "    linf = z[0]-b\n",
    "    lsup = z[0]\n",
    "    return {'a':a, 'b':b, 'c':c, 'linf':linf, 'lsup':lsup}\n",
    "\n",
    "# Golden and Alt, 1979\n",
    "golden = estimates(2*z[0]-z[1])\n",
    "\n",
    "# Zanakis (1979)\n",
    "zanakis = estimates((z[0]*z[-1]-z[1]**2)/(z[0]+z[-1]-2*z[1]))\n",
    "\n",
    "# Cooke 1979.\n",
    "cooke = estimates(2*z[0]-(math.e-1)*(sum(math.exp(-(i+1))*z[i] for i in range(k))))\n",
    "\n",
    "# Scipy Maximum Likelihood Estimate (MLE)\n",
    "cMLE, aMLE, bMLE = weibull_min.fit(z, cooke['c'], loc = cooke['a'], scale = cooke['b']) # shape, location, scale\n",
    "MLE = {'a':aMLE, 'b':bMLE, 'c':cMLE, 'linf':z[0]-bMLE, 'lsup':z[0]}\n",
    "\n",
    "df = pd.DataFrame([golden, zanakis, cooke, MLE], \n",
    "                  index = ['Goden and Alt (1979)', 'Zanakis (1979)', 'Cooke (1979)', 'MLE'])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNX5+PHPM5OVsJOwBgwQCEsS9n3HhR3E5StYXKotP+1XrVqtWFvbWq3aWrF1rTtaFfyCCBaVCgkqIoEAYQ1gWJQga9izZ+b8/riXNISEDCHJnUme9+s1L2bOPffOM8NNntxzzj1HjDEopZRSLqcDUEop5R80ISillAI0ISillLJpQlBKKQVoQlBKKWXThKCUUgrQhKCUUsqmCUEppRSgCUEppZQtyOkALkZkZKSJiYlxOgyllAoo69atO2qMiaqoXkAlhJiYGFJTU50OQymlAoqIfO9LPW0yUkopBWhCUEopZdOEoJRSCgiwPoSyFBYWkpmZSV5entOhqBLCwsKIjo4mODjY6VCUUj4K+ISQmZlJgwYNiImJQUScDkcBxhiysrLIzMykffv2ToejlPJRwDcZ5eXl0axZM00GfkREaNasmV61KRVgAj4hAJoM/JD+nygVeGpFQlAqUHi9TkegVPk0IVSBgwcPMm3aNDp27EifPn0YP348O3fuvOjjrFixgokTJ1Y6jgMHDhTvn5WVxahRo6hfvz533XXXOfXmzZtHYmIi3bt356GHHiouv+++++jZsyc9e/akc+fONG7cuHjbQw89RHx8PPHx8cybN6+4fNq0aXz33XeVjrmuMAbeew/atYMRI2DHDqcjUup8mhAukTGGqVOnMnLkSHbt2sW6det48sknOXToUI3H8uyzz/Lzn/8csEb5/OlPf+KZZ545p05WVhYPPvggy5cvZ+vWrRw8eJDly5cDMHv2bNLS0khLS+Puu+/mmmuuAWDJkiWsX7+etLQ0UlJSeOaZZzh16hQAd955J3/5y19q8FMGnk2brCQwYwZERVmve/SAp56CoiKno1PqvzQhXKLk5GSCg4O54447ist69OjBsGHDMMbw4IMPEh8fT0JCQvFf1uWVl7R27Vp69erFrl27yM7O5rbbbqN///706tWLRYsWlRnLggULGDt2LAAREREMHTqUsLCwc+rs3r2bTp06ERVlTWtyxRVXsGDBgvOO9cEHHzB9+nQAtm3bxvDhwwkKCiIiIoLExEQ+//xzAIYNG8ayZcso0t9sZcrMhMGDIT0dXnsN1q2znk+cCA8/DLfc4nSESv1XwA87LeneeyEtrWqP2bMnPPdc+du3bNlCnz59ytz20UcfkZaWxsaNGzl69Cj9+vVj+PDhrFq1qszys1atWsXdd9/NokWLaNeuHb/5zW8YPXo0b775JidOnKB///5cccUVREREFO+zZ88emjRpQmho6AU/T2xsLDt27GDv3r1ER0fz8ccfU1BQcE6d77//nj179jB69GjASnB//OMf+dWvfkVOTg7Jycl069YNAJfLRWxsLBs3biz3e6jLHnrIugpISYEOHayyli1h/nyYNQuefhoefNA6z5Ryml4hVKOVK1cyffp03G43LVq0YMSIEaxdu7bccoD09HRmzpzJJ598Qrt27QD4z3/+w1NPPUXPnj0ZOXIkeXl5/PDDD+e814EDB4r/6r+QJk2a8PLLL3PDDTcwbNgwYmJicLvd59SZO3cu1113XXH5VVddxfjx4xk8eDDTp09n0KBB5+zTvHlzfvzxx0v6rmqjb76B99+3fuGfTQYlzZoFjRvDo4/WfGxKlaVWXSFc6C/56tK9e3fmz59fZcdr1aoVeXl5bNiwgdatWwNWE9OCBQuIi4srd7/w8HCfx/1PmjSJSZMmAfDqq6+WmRBefPHFc8oeeeQRHnnkEQBuvPFGOnfuXLwtLy+P8PBwn967rvB64Ze/hDZtrF/8ZWnc2EoWjzxiXUEMGFCzMSpVml4hXKLRo0eTn5/Pq6++Wly2adMmvv76a4YNG8a8efPweDwcOXKEr776iv79+5dbDtC4cWOWLFnCww8/zIoVKwAYM2YMzz//PMYYADZs2HBeHJ07d2bv3r0+xXz48GEAjh8/zksvvcTPfvaz4m3bt2/n+PHjDBo0qLjM4/GQlZVV/Nk2bdrEVVddVbx9586dxMfH+/TedcVbb1n9BX/5C5Ro2TvPPfdAZCT87nc1F5tS5TLGBMyjT58+prRt27adV1bT9u/fb66//nrToUMH061bNzN+/Hizc+dO4/V6zQMPPGC6d+9u4uPjzdy5c40xptzy5ORkM2HCBGOMMd9//73p1q2bWb16tcnJyTEzZ8408fHxplu3bsV1Shs9erT57rvvil9fdtllpkmTJiYiIsK0adPGbN261RhjzLRp00zXrl1N165dzQcffHDOMX7/+9+bhx566Jyy3Nzc4voDBgwwGzZsKN528OBB069fvzLj8Yf/Gyfk5hrTvLkxQ4YY4/VWXP+ZZ4wBY778svpjU3UTkGp8+B0rxv6rMxD07dvXlF4gJz09na5duzoUkX9ZuHAh69at4/HHH6+x95w9ezYNGzbk9ttvP29bXf2/mTcPpk2DL76AK66ouH5uLnTsCImJYA/eUqpKicg6Y0zfiuppk1EtMnXqVGp6idHGjRtzi46dPMecOdC2LdiDtCoUHm4NP122DI4fr97YlLoQTQi1TMn+gJrw05/+lKCgWjU24ZIcOABLl8JNN4HrIn66rr4aPB5YsqT6YlOqIpoQlKpC771njTC6+eaL269fP2jVCsq551CpGqEJQakqYozVXDRwIFxghHCZXC6YPBk++wx01nDlFE0ISlWRDRtgy5bKT0cxZQpkZ0NSUtXGpZSvNCEoVUXmzIHQULjhhsrtP3o01K8PH39ctXEp5StNCFVARJgxY0bx66KiIqKiooqnon777bfPm4IaICYmhoSEhOIpp++5554ai1lVrcJCa5qKyZOhSZPKHSM0FMaPh8WLdd0E5QwdHlIFIiIi2LJlC7m5uYSHh/PFF1/Qpk0bn/ZNTk4mMjKymiNU1W3lSjh6FOwJYittyhT48ENrKosSN4srVSP0CqGKjB8/niX2mMGSU0eruuHTTyE4GK688tKOM348BAXpaCPljNp1heDE/Ne2adOm8dhjjzFx4kQ2bdrEbbfdxtdff13hfqNGjSqeXO6WW27hvvvuu+SQVdXanHmywjoLF9enz0Ave07kwImy6yREN4LkZDh8GDp3hk6drE6DEho3hpEj4d//thbQUaom+ZQQRGQs8HfADbxujHmq1PZQ4B2gD5AF3GCM2SsiVwJPASFAAfCgMSbJ3mcF0ArItQ9zlTHm8CV/IockJiayd+9ePvjgA8aPH+/zftpkFPh+zBR27XQz9YaCcutIQQHcfTe88MK5G2bMgLffhhIzzo4cCb/9LRw7Bk2bVk/MSpWlwoQgIm7gReBKIBNYKyKLjTHbSlS7HThujIkVkWnA08ANwFFgkjHmRxGJB5YCJRvXf2KMOXdyokvhxPzXJUyePJkHHniAFStWFM8Oqmq/lcnBAAwbXfaqcUEH9tPuzp/CujVw//3WuNTvvoMVK6wE0bSpde6KADBkiLXf6tVWE5JSNcWXK4T+QIYxZjeAiMwFpgAlE8IU4A/28/nACyIixpiS8zRvBcJFJNQYk3/Jkfuh2267jcaNG5OQkFA8dbWq/VYmB9G6rZeYjucPDXKdPEHs5CtxnT5l9RZff721ITERrr3W6niYPRsuu8xKFlh3LbvdsGqVJgRVs3xJCG2AfSVeZwKll/IormOMKRKRk0AzrCuEs64F1pdKBm+JiAdYADxuAmnq1TJER0eXO3T07bff5uMSA8xXr14NnNuHkJiYyDvvvFP9gaoqU5APKd8EMfm6grN/4J+j5VOPEXT4ILs+/oLYSWXMdvfMM7BvH/zqV9CuHVx3HRERVtfVN99Uf/xKlVQjncoi0h2rGemqEsU/McbsF5EGWAnhJqx+iNL7zgRmAsVLSvqbM2fOnFc2cuRIRo4cCcCtt97Krbfeel4dXxe0Uf5r/Ro3uTnC0FHnNxfVW7eGpu+9RdZtd5Dbq5z1pl0uePdd+OEH+N//hQkTIDycIUPg9det+xuCg6v5Qyhl82XY6X6gbYnX0XZZmXVEJAhohNW5jIhEAwuBm40xu87uYIzZb/97Gngfq2nqPMaYV40xfY0xfX1ZM1ipmvR1UjAhoYZ+g0slhMJCWs+6l8KWrTn0wG8ufJCwMHj6aWv00Zw5AAweDDk5sGlTNQWuVBl8SQhrgU4i0l5EQoBpwOJSdRYDZ2dwuQ5IMsYYEWkMLAFmGWOKL4BFJEhEIu3nwcBEYMulfRSlat7KFUH0HVBEvXrnlke+8TLh27dx4I9P4a3foOIDjRhhLar8179CURGDB1vF2mykalKFCcEYUwTchTVCKB340BizVUQeE5HJdrU3gGYikgHcD5xdVvwuIBZ4VETS7EdzIBRYKiKbgDSsK4zXKvshArzroVaqC/8n+/cJezLcDC01ush9LIsWzz7FqSvGcmrsRN8OJgKzZsHu3TB/Pm3bWovsrFpVDYErVQ6f+hCMMZ8Cn5Yqe7TE8zzg+jL2exwobz3HchpVL05YWBhZWVk0a9YMKatXT9U4YwxZWVmEhYU5HUq1WrvK+vEZOPTchND0g3dw5eZw8KFHKbOnuTyTJ0OXLtYdaTfcwODBoglB1aiAv1M5OjqazMxMjhw54nQoqoSwsDCio6OdDqNapa4OoklTLx07lxhuWlhI0zmvc2boCPK7dLu4A7pc8Otfw223wdKlDB48lnnzrEFIbdtWvLtSlyrgE0JwcDDt27d3OgxVB6WmBNF7gOeci4CGS5cQcmA/Pz7+18od9Cc/gd/9Dp55hiFPjwWsZqPKTqmt1MXQye2UqoQD+4Uf97noO/Dc5qLIN/9JfrsYTl8+pnIHDgmBn/0MkpJIbH6QevW0H0HVHE0ISlVC6mrr4rpkQgjbnEbE2m/JuvXn58xNdNGuvRaMIXjJx/TvryONVM3RhKBUJaSuDqJBI0OnLv/tP4h885946kVw/IYZF9jTB/Hx1kyoH33EwIHWvQgF5c+bp1SV0YSgVCWsW+2md78iXPZPkPtYFo0WL+DE9Tfibdjo0g4uAtdcA8nJ9Ot4jMJCSE+/9JiVqogmBKUu0uGDwg973ec0FzX6dDGuggKOTb+pat7k2muhqIiBRz8Bqn6ZD6XKEvCjjJSqaf/tP/AUlzVaNJ+82M7kdUsodz9fFtop1rITca2jiVg+l7Cwm1n2dQG9L88r3pwQfYlXIUqVQa8QlLpI61LcRNQ3xHW3EkLQgR+JSFnFySnXXtyNaBciwqlxE2nwdTKJnU6wI/0SOqmV8pEmBKUu0rqUIHr1KyLIvr5utORjxBhOTLqmSt/n5LjJuPLzmdZwCTu2uqgDs4Eoh2lCUOoiZB0Vdn93bv9B40ULyI1PpKBjpyp9r5y+AyiMjOKKUx9z6qSLgz/q1CyqemlCUOoibFhrNd30GWA1FwV/v5d6aes4Mfnaqn8zt5tTYyYQl/Ef3BSxfas2G6nqpQlBqYuQlhpEcIiha7yVEBp/8hEAJydNrZb3yx4ygpDc0/RhHTu3aUJQ1UsTglIXYdN6N90SPISEWq8bLZpPdt8BFEZXz2p+ZwYNBWBqkyS2a0JQ1UwTglI+KiyAbZvd9OhjXR2EfreD8O3bqu3qAMATGUVe5y5cHvQlO7bqj6uqXnqGKeWj9K1uCvKFHn2sDuUGX3wGwMlxk6r1fbMHDiXxxEoO/uDlzOlqfStVx2lCUMpHG1OtJpseva0rhIbLlpLbPYGiVm2q9X3PDB5GaGG21Y+g9yOoaqQJQSkfbVzvplUbL81bGtzHj1FvXQqnrhhb7e+bPcBaYHkUyezQfgRVjTQhKOWjjeuC6NHbbi5K+g/i9XL6inHV/r5WP0JXrghK1qGnqlppQlDKBwd/FA4dcNGjr9Vc1GDZUgqbtyA3sWeNvP+ZwUMZ7P2G3VuKKq6sVCVpQlDKBxvXWfNU9OjtQQoKaPDlck6Pvori+a+rWfagoYR7c2i0YwNFmhNUNdGEoJQPNq5zExpqiOvmod6ab3GfPlX5ZTIrIXugdT/C0MIV7N+nP7aqeuiZpZQPNq530y3RQ3AINFz2Od7QUE4PH1Vj7+9p2ozjl3VnJCvYvVN/bFX10DNLqQrk50H6Fjc9+3rAGBos+5zswcMx9SJqNI68oUMYwjfs2e6tuLJSlaAJQakKbNvspqhQ6NG7iNBd3xH6/R5O1WBz0VmFA/tTj1y863U9TVU9NCEoVYFN662hnom9PdRfsQyA06OvrPE4cnr1BaDp9rU1/t6qbtCEoFQFtqRZN6RFNjc0+DKJ/A6xFLa9rMbjKGx3GafCIok5lIpXW41UNdCEoFQFNqcFkdCrCMnLI2L1N5weMdqZQEQ4GNOHvp41/PCDMyGo2s2nhCAiY0Vkh4hkiMisMraHisg8e3uKiMTY5VeKyDoR2Wz/O7rEPn3s8gwR+YdIVS1Gq1TVOXQIfsx0Ed/DQ73UFFx5uZwZ7lBCAHJ69KEr6exMPeVYDKr2qjAhiIgbeBEYB3QDpotIt1LVbgeOG2NigdnA03b5UWCSMSYBuAV4t8Q+LwM/BzrZj+qfFEapi7RmjfVvQi8PDb5KwhscTLa9RoETgkb0wYXh+LJUx2JQtZcvVwj9gQxjzG5jTAEwF5hSqs4UYI79fD5wuYiIMWaDMeZHu3wrEG5fTbQCGhpjVhtjDPAOcPUlfxqlqtiaNeB2G7omeKj/VRI5fQfgjajvWDyuob0AcKeucSwGVXv5khDaAPtKvM60y8qsY4wpAk4CzUrVuRZYb4zJt+tnVnBMpRyXkgKxcV4aZh8mfOtmR5uLADxNmrIvLJbI3SmOxqFqpxrpVBaR7ljNSP+vEvvOFJFUEUk9cuRI1QenVDm8Xli7FhJ6FVF/5QoATjucEAB+aN2HzifWYIzTkajaxpeEsB9oW+J1tF1WZh0RCQIaAVn262hgIXCzMWZXifrRFRwTAGPMq8aYvsaYvlFRUT6Eq1TV+O47OHEC4nt6qP/lcoqaNiMvPtHpsDjetR+tzY8cTM2suLJSF8GXhLAW6CQi7UUkBJgGLC5VZzFWpzHAdUCSMcaISGNgCTDLGPPN2crGmAPAKREZaI8uuhlYdImfRakqVdyh3KOI+l8lc2bYqBqb3fSCBln9CIc+0X4EVbUqPLvtPoG7gKVAOvChMWariDwmIpPtam8AzUQkA7gfODs09S4gFnhURNLsR3N72y+A14EMYBfwWVV9KKWqQkoK1K8P3TybCT5yuEYns7uQhlfGU0AwRd9oP4KqWkG+VDLGfAp8Wqrs0RLP84Dry9jvceDxco6ZCsRfTLBK1aQ1a6BvX2j0TTKAdYXgB5q0DmVrUA8apOsVgqpafnD9q5T/ycuDtDQYMADqf/MVebGdKWrV2umwiu1pPoC2h1LB43E6FFWLaEJQqgxpaVBYCAN6FVAvZRXZQ0Y4HdI5TnfuSz3vGcyOnU6HomoRTQhKleFsh/KQ4DW4c7I5M2S4swGVEjLA6lg+uWKDw5Go2kQTglJlWLMGWreG5puXY0SKl7D0F5HDu5FPCKe/0oSgqo4mBKXKsHYt9O8PLF9ObnwPPE2aOB3SOTp3D2YL8chGTQiq6mhCUKqU48dh504Y0jMbVq8me6h/9R8AtG0Lm929aLI3Db1lWVUVTQhKlZJqTyQ6OmQlFBZyxs86lMG6P+5Ay15E5GVBpt6xrKqGJgSlSjnbodztwHIICSG7/0BnAypHTpzVscwGbTZSVUMTglKlrFkDcXEQ9s1yGDQIE17P6ZDKFNovES9C0VpNCKpqaEJQqgRjrCkrRvU4Zv3lPdr52U3L0yGxPjvpTM4qTQiqamhCUKqEzExr2czJDZKt7HD55U6HVK64ONhAL4I2aUJQVUMTglIlrF1r/dv7RBJERNhjT/1T585WQqh39AfIynI6HFULaEJQqoQ1ayA4GKK2JMHw4dYLP9WgAexrZncsp6U5G4yqFTQhKFXCmjUwuusBXDu2wyj/mN30QvK66kgjVXU0IShl83isexBuaG5Nd+3PHcpntUqIZL9EYzQhqCqgCUEp244dcPo0DClIgsaNoWdPp0OqUFwcrDO98KRqQlCXThOCUrazN6TF7EmGESPA7XY2IB906WJ1LLszdkBOjtPhqACnCUEp25o10D1iLyH7dgdEcxFYVwhp9ES8XtiyxelwVIDThKCULSUFbm4bOP0HAO3awfaQHtaLTZucDUYFPE0ISgG5udbv0yuCkiEqCrp3dzokn7hcEBoXQ467viYEdck0ISiFNWqzqMjQ9UCSNdxUxOmQfNa5i4vtwQmwcaPToagApwlBKazmok58R3jW/oBpLjorLg7W5PfAbNqkayOoS6IJQSmshHBdkyTrRQAmhI0mETlxQtdGUJdEE4JSWCOMJkYkQ5s2EBvrdDgXJS4ONpFovdB+BHUJNCGoOu/IEdi7x0vP40nW7KYB1H8AVkLYTIL1QvsR1CXQhKDqvJQUiGcL9bKPBlxzEUDDhhDRsiFH6rfXKwR1STQhqDpvzRq4Quz+gwCY0K4scXGQHpKoCUFdEk0Iqs5LSYHJDZKsvoN27ZwOp1I6d4Y1uYnWhEy5uU6HowKUTwlBRMaKyA4RyRCRWWVsDxWRefb2FBGJscubiUiyiJwRkRdK7bPCPmaa/WheFR9IqYvh9cK6lCIG5H4ZkM1FZ8XFwercROsDbdvmdDgqQFWYEETEDbwIjAO6AdNFpFuparcDx40xscBs4Gm7PA/4HfBAOYf/iTGmp/04XJkPoNSlyMiADifXE154KuATwkZ0Cgt1aXy5QugPZBhjdhtjCoC5wJRSdaYAc+zn84HLRUSMMdnGmJVYiUEpv5OSAqOx+w9GjnQ0lksRFwe76UBhSD1NCKrSfEkIbYB9JV5n2mVl1jHGFAEngWY+HPstu7nodyIBNtZP1QopKXClOwkTHw8tWjgdTqXFxIAryM3BZvE69FRVmpOdyj8xxiQAw+zHTWVVEpGZIpIqIqlHjhyp0QBV7Zf6TT5DzEokgJuLwFr6uWNH2B7aw7pC0CksVCX4khD2A21LvI62y8qsIyJBQCMg60IHNcbst/89DbyP1TRVVr1XjTF9jTF9o6KifAhXKd/k5EC9TasJ8+ZaN6QFuLg4WJOXCFlZcOCA0+GoAORLQlgLdBKR9iISAkwDFpeqsxi4xX5+HZBkTPl/oohIkIhE2s+DgYmAru6halRqKozwJmFcLhg+3OlwLllcHCRn2VNYaLORqoQKE4LdJ3AXsBRIBz40xmwVkcdEZLJd7Q2gmYhkAPcDxUNTRWQv8Cxwq4hk2iOUQoGlIrIJSMO6wnit6j6WUhX79lu4nOUU9ehjraEc4Dp3hnWF9hQWmzc7G4wKSEG+VDLGfAp8Wqrs0RLP84Dry9k3ppzD9vEtRKWqR9pXp/gVqwka95DToVSJuDg4QRNyI6MJ15FGqhL0TmVVJxkDQau+IggPXHGF0+FUibg469+DzXUKC1U5mhBUnbR3L/Q5sYzC4HAYNMjpcKpEVJTV8rUzNBG2b4eCAqdDUgFGE4Kqk779Fq5gGbl9h0FYmNPhVAkR6yohtSARCguteY2UugiaEFSdtHXZAeLZSsSU2tFcdFbnzpB81O5Y1mYjdZE0Iag6yb1iufXvlYF//0FJcXHw5aE4THCwjjRSF00TgqpzcnMhdu8ycsKaQs+eTodTpeLioIhgctt30ysEddE0Iag6J3WtYbRZxok+l4Ordv0IdOli/Xu4eYImBHXRfLoPQalAtDnzZJnlX7+/l2HsZ+eoIeXWCVSxsVaO2xmWSMz+f8GxY9C0qdNhqQBRu/48UsoH9b/6EgDv2JHOBlINwsKgfXt7pBFoP4K6KJoQVJ1iDMRmLONgRAyFl8U4HU616NIFlh+xE4I2G6mLoAlB1Snfby9iWGEymQm1a7hpSV27wqrdLTHNmmlCUBdFE4KqU44uSKUBZzATA3v9gwvp0gXy8oW8zonaZKQuiiYEVafUX5FEIUGETx3qdCjV5uxIo4PN7YTg9TobkAoYmhBUnRK3exnpTQZiGjZ0OpRqczYh7AxNtFYB2r3b2YBUwNBhp6rOOL7tKMML1/NZj98T7XQwl6ii4bJNmjYg6WQXxgDfJ63iVFjZqw0mRDeqhuhUoNIrBFVnnJ5rDTd1Ta69/QdntY/18nVWAsblImybLkaofKMJQdUZDb9azhEiiZqc4HQo1a59rJfteyLI7xBLuCYE5SNNCKpu8Hrp8v1y1kdeQXBo7T/t23f0cDzLxenYeMLSNSEo39T+nwylAG/qFqKKDnGgd+2a3bQ87WOtkUX7IhMJydyH6+QJhyNSgUATgqoTcj9cAUDQlJGOxlFT2nf0AJAeat2xHJa+1clwVIDQhKDqhKYrv2AjicSOLnu0TW3TKtoQEmpYm29N7x2uzUbKB5oQVK3nOnmCDvtXkRI1joj6TkdTM9xuiOngZcOBNhQ1baYjjZRPNCGoWi90WTJBeDg8aIzTodSo9h097Mlwk9ctnrCtOoWFqpgmBFXrmQX/IYumNL26j9Oh1KiYWC/797k40zmesJ3pUFTkdEjKz2lCULWbx0PLtV/wuYyj10DjdDQ1qn1HL16vsC8yEVd+PqG7M5wOSfk5TQiqVgvfuJ4GeUdJazeO+g2cjqZmdehkjTTaGtQDQPsRVIU0IahaLfyzpRThJu/y2j9dRWntO3pxuQxrT3fFGxysN6ipCmlCULVa6Gf/4RuG0H1U7Z3dtDyhYdA2xsuOjHDyO8XpFBaqQj4lBBEZKyI7RCRDRGaVsT1URObZ21NEJMYubyYiySJyRkReKLVPHxHZbO/zDxGRqvhASp0VdGA/kd9v4jPXeHr1q5sdqrFxXjJ2uKyRRpoQVAUqTAgi4gZeBMYB3YDpItKtVLXbgePGmFhgNvC0XZ4H/A54oIxDvwz8HOhkP8ZW5gMoVZ4GSV8A8F3XsdSLcDgYh8TGedi318WZTvEEHz6I++gRp0NSfszeC54DAAAZM0lEQVSXK4T+QIYxZrcxpgCYC0wpVWcKMMd+Ph+4XETEGJNtjFmJlRiKiUgroKExZrUxxgDvAFdfygdRqrR6Sz9nDzE0G9XZ6VAcE9vZGmm0q5HVsazNRupCfEkIbYB9JV5n2mVl1jHGFAEngWYVHDOzgmMqVXnZ2TRcuYJPmES/IR6no3FMbJz12Td4rSkswrZucjIc5ef8vlNZRGaKSKqIpB45ope7ykeff05QYR6fBF1Nzz51NyG0a+8lKNiwJTOSgrbtCN+U5nRIyo/5khD2A21LvI62y8qsIyJBQCMgq4JjllzFsKxjAmCMedUY09cY0zcqqm5MTKaqwMKFHHc341TvwYSFOx2Mc4KDreGnGTvd5Cb0JHyzJgRVPl8Swlqgk4i0F5EQYBqwuFSdxcAt9vPrgCS7b6BMxpgDwCkRGWiPLroZWHTR0StVloICvJ/8m4WeyfQf4XQwzuvY2UPGDje58T0I/X6Pro2gylVhQrD7BO4ClgLpwIfGmK0i8piITLarvQE0E5EM4H6geGiqiOwFngVuFZHMEiOUfgG8DmQAu4DPquYjqTovORnXqZMsZCpDRtTN4aYlxcZ5+XGfi+Od7KmwtR9BlSPIl0rGmE+BT0uVPVrieR5wfTn7xpRTngrE+xqoUj5buJC8oAhSG1zO4wmFTkfjuE52x/LW8N50BcI3bSR78HBng1J+ye87lZW6KF4vZtEi/uMaR+8RQbj0DKdjZ2s5zW0HmlPQJlr7EVS5fLpCUCpgrF6NHDzIXKYyWJuLAGjTzktYmCnuR9CEoMqjfz+p2mXhQjzuYJYwgcHDNSGAtXpah04eMna4yE3oSeieXbhOn3I6LOWHNCGo2sMYWLCA1Iaj6dCzEZHN69b6BxdizWnkJjfR7ljeoh3L6nyaEFTtsXYt7NnDaydvYEzdWi2zQh3jPBw57OJwWzshaLORKoMmBFV7zJ2LJyiE+d6pmhBKiY2zOpbTs1pS2LI14Zs3OhyR8keaEFTt4PHAvHlsjh5HUURjhgxxOiD/0qmLNfR0Z7qb3IQehOkVgiqDJgRVO6xcCT/+yBtnpjF6NISEOB2Qf2nR0tCkqZftW6x+hNDdGbjOnHY6LOVnNCGo2uGDD/CG1ePNo5OYNMnpYPyPCHSN95C+xZrTSIwhbOtmp8NSfkYTggp8hYUwfz7bYieTKxGaEMrRNcFDxk4XJ+OsjuV6G9c7HJHyN5oQVOBbtgyysngjZzoDB0LLlk4H5J+6xnsoKhS2H29l3bG8IdXpkJSf0YSgAt/cuXgbNual3WOYUnotP1Wsa7w90miLm5ze/am3bq3DESl/owlBBbbsbFi4kO3drqGAUK7WhVjL1aadl/oNDOmb3eT06UfIgf2QmVnxjqrO0ISgAtuCBXD6NK8X3kJcHMTFOR2Q/3K5oEt3D9u3WlcIAKxe7WxQyq9oQlCB7a238HSI5fm0YXp14IOu8R52bnNzJi4Bb2gYfPut0yEpP6IJQQWuXbtgxQq29LmVIo9o/4EPuiZ4yMsT9uwLs+Y10oSgStCEoALX22+DCK/k3kKLFjBggNMB+b+u8dYdy1bHcj9Ytw7y8x2OSvkLTQgqMHk8MGcOnsuv4l8ropk8GV0MxwcxHa21EbZvsTqWKSiADRucDkv5Cf0RUoFp+XLYt4818bdx5gzccIPTAQUGtxs6d7PuWM7p1c8q1I5lZdOEoALTW29BkyY8t3syLVvCyJFOBxQ4usZbI40KmreCdu20H0EV04SgAs/Ro7BwIfnX3sjHn4cxbZr1l6/yTdcED2dOC5nfu2DQIE0IqpgmBBV4Xn8d8vP59LI7KSiA6dOdDiiwlOxYZtAg2LcP9u93OCrlDzQhqMBSVAQvvwyjRvHSl93p2BH69XM6qMAS29lLaKhh8wY7IYBeJShAE4IKNJ98Aj/8wPEZd5OUBDfeaE3trHwXHALxPT2sX+uGnj0hLMxaT0LVeZoQVGB5/nlo145/nZyE16vNRZXVq18R27e4yS4MgaFDISnJ6ZCUH9CEoALH1q2QnAx33sl784Lo2RO6dnU6qMDUq5+HoiJhzRrg8sth82Y4fNjpsJTDNCGowPHCCxAayncjfkZKitVcpCqnR+8iRAzffAOMHm0VJic7GpNyniYEFRiOH4d334Ubb+Tl/4skOBhuucXpoAJXw8bQsbPX6jro3RsaNbJu9lN1mk8JQUTGisgOEckQkVllbA8VkXn29hQRiSmx7WG7fIeIjClRvldENotImojo0k3qwp5/HrKzyb/zXubMgalToXlzp4MKbL36FbFqFXgkCEaM0ISgKk4IIuIGXgTGAd2A6SLSrVS124HjxphYYDbwtL1vN2Aa0B0YC7xkH++sUcaYnsaYvpf8SVTtdeYM/P3vMGkS83cmcuwYzJzpdFCBr3d/D6dPw5YtWP0Iu3fD3r1Oh6UcFORDnf5AhjFmN4CIzAWmANtK1JkC/MF+Ph94QUTELp9rjMkH9ohIhn08HfSsyrQ58+R5ZZH/fJ5Wx46R8bN7mP3nItpeJkR2OsNmXezrkvTsWwRYI057nO1HSEqC225zMCrlJF+ajNoA+0q8zrTLyqxjjCkCTgLNKtjXAP8RkXUion/vqTJJXh6Rr77ImSHD2dpgAOtSgrj2JwU6s2kVaB1taNPGvgWhe3do0UKHn9ZxvlwhVJehxpj9ItIc+EJEthtjvipdyU4WMwHatWtX0zEqhzX5v/cJPnyQfX//JwveDyEo2DDl+kKnw6oVRGDIEKyRRiLWaKPly8EYvduvjvLl76z9QNsSr6PtsjLriEgQ0AjIutC+xpiz/x4GFmI1JZ3HGPOqMaavMaZvVFSUD+GqWqOwkKiXnyOnZx+yeg9n8fxgLh9TSLNI43RktcbQodZURj/8gJUQDh6E7dudDks5xJeEsBboJCLtRSQEq5N4cak6i4GzgwCvA5KMMcYun2aPQmoPdALWiEiEiDQAEJEI4Cpgy6V/HFWbNJ33L0L2/cDhex5g0fxQTp5wccPNBU6HVasMHWr9u3IlVscy6GijOqzChGD3CdwFLAXSgQ+NMVtF5DERmWxXewNoZnca3w/MsvfdCnyI1QH9OfC/xhgP0AJYKSIbgTXAEmPM51X70VQgk5xsms9+iuy+Azg+cixz/hlCQq8i+gz0OB1arZKQYN2CkJQEtG8PHTvCp586HZZyiE99CMaYT4FPS5U9WuJ5HnB9Ofs+ATxRqmw30ONig1V1R9SrLxJ8+BA//PMdln8eQuYPbu7/bbY2bVexoCC46iorBxgDcvXV8I9/wMmTVqZQdYqO1VB+x330CJGv/IOTYyaQ3WcAb74cymUdPIy6qsjp0GqlCRPgwAF7aeWpU6GwUK8S6ihNCMrvNP/7X3Dl5XJo1u9J+cZN+mY3t96Rr6uiVZNx46xBRUuWYK2P0LIlfPSR02EpB2hCUH4lZHcGzf71Fsem3UR+bGfeeimUyOZeJk7VoabVpXlza5GhJUsAlwumTIHPPoPcXKdDUzVME4LyH8bQ5pFf4Q2vx+H7ZrFpg5tvvw5mxu35hIY5HVztNmECrFkDR45gNRtlZ8OyZU6HpWqYJgTlP95/n/orv+TgQ49S2Lwls58Io2mkV4ea1oAJE6xO5c8+A0aNsjqUtdmoztGEoPzDsWNw333k9OzDsRk/5avlQaxLCeKOe/OJqO90cLVfr17WzBVLlgAhITBxIixebK1hreoMTQjKP8yaBceOsf+p2RQZN889GcZl7T1ce6NeHdQElwvGj4elS61BRkydaiXpr86bTUbVYpoQlPO+/BJeew3uvZe87oksnh/Mrp1u7n4oj+Bgp4OrOyZMsG4/+PZbYOxYCAuDBQucDkvVICcnt1PKWgltxgyIjYU//IHcAx5eeiaMhF5FXDlemyuqW8npxlt3g6DghrzxXgFNOhTR9qoJNHjvPdLv/S0mvN4Fj5MQrTex1QZ6haCcY4y10s3Bg/DBB1C/Pq/MDuPwIRf3P5KndyXXsPoNYOjIIpYsDKawEI7N+Cnukydp9MlCp0NTNUQTgnLOG2/A/PnwxBPQty/r1sGcf4ZwzfQC+gzQOYuccM30ArKOuFiZFET2wCHkxXam2b/ecjosVUM0IShnpKfDL39pzbD5wAMUFsLtt0OzSMP9j+gNUU4ZOqqIqOZePpobAiIcm/FT6m1IJWzLRqdDUzVAE4KqeVlZMGkS1K8P77wDLhfPPAMbN8JvHs+loTZHOyYoCCZfX8DK5CCOHBKOXzsdb2gYTd972+nQVA3QhKBqVkEBXHMNZGbCxx9D69akp8Mf/wjXXQeXj9OOZKddfUMhHo+weH4I3saNOTlpKo0X/h+uM6edDk1VM00IquYYA3fcYY1tf/NNGDSIM2fg2mutG2Off97pABXAZe299BlQxMfzgjEGsm66DXf2GRp/PN/p0FQ102GnqkqUHL5YnubPPkmLt97i0L2/5vDwCZh9J/n1L8LZsSOYVz/I5kiRdiT7i6nTCvjtffVYl+Km74C+5HZPIPK1Fzl2wwz05pDaS68QVI2I+vtfaTH7aY5dfyOH75sFwHtvhrD03yHc81Ae/QdrMvAnV04opH4Dw7w5VufyofsfJnR3Bk3nvut0aKoaaUJQ1S7qhWdp+cwTHL/2Bvb/9XlwuVizys2zj4cxekwhP71Tp6fwN+HhMP3WfJb+O4Rtm1ycvnIc2f0G0Xz2U7iyzzgdnqommhBU9TGG5n/7My2ffozjU/+HzL+9BG43mza4uee2CC7r4OVPz+boDWh+6tY78mncxMtzT4WBCAcf+SPBRw4T+dqLToemqokmBFUtJC+PtvfMpMVzf+HY//yEzNkvg9vNznQXv7ipHs0ivbz6fjYNGjodqSpPg4Yw85f5rP46mFVfBpHTpz8nx00i8pXncR894nR4qhpoQlBVzn30CO2nT6Hxx//HwYceZf8zL4Dbze7vXMy8MYKwcHj1g2yiWhinQ1UV+J8ZBbRu6+W5J8PweuHgQ4/iysulxd+edDo0VQ00IagqVf/LJDqNGUr45o18/8ocjtx1P4iQ8o2bm6bWRwRe+yCbNm01GQSCkFC4+8E8tm9189miYAo6diLr1pk0+9ebNFi6xOnwVBXThKCqhOTl0eoPs2g/4xo8jZuwa9F/ODVhCgAfzwvmzhkRNG/h5V+LztA+1utwtOpijJtSSNcED3/5Qxg/ZgoHZ/2enISetL3/ToK/3+t0eKoKaUJQl27JEjpdNYTIN17h6E9nkvHvZPK6J5KbC0/+LoxHH6hHv0FFzPnojF4ZBCCXC556PofCQuG+mRHkmjB+eOVtEOGyO25G8vKcDlFVEU0IqvLS061ltiZOxLhc7HlvIQce+wsmPJy0VDf/M6Y+H7wdyozb83lhTo7OURTA2nf08ue/55C+2c2ffhNOQdsY9s1+hfAtm2j9m/vBq1d9tYHeqawu3vr18OST1mpaDRrA3/5GxtU3YUJCyDoivPqPUOa9E0LL1obX553Rm85qiZFXFnHn/Xm8/GwYnbt4uOX/jePQvb+mxXN/AQqtiQrDwpwOU10CTQjKN3l51mR0r78Oy5dDw4bw8MNw770QFcWpbSd55/lQ3nk1lIJ8uH5GAb+clUf9Bk4HrqrS//tlPjvT3fzt8XAO7Hfxq98+jLdBQ1r96bdw4AAsWgRNmzodpqokTQiqfPn5sGKFlQjmzbOWu7zsMvjzn+EXv4BGjUhLg5d/C+/+qyG5OcKYiQX874P5xHTQJoTayOWCv76Uw+wnwnj39VB2bHPzzCt3U9iyNdH33UFh774c+P2fOX3FWCp7x6Eux+kcn/oQRGSsiOwQkQwRmVXG9lARmWdvTxGRmBLbHrbLd4jIGF+PqRzg8VjNQc89Z01RHRlpLbb+zjswZgx88QVm1242jH2YP8xuRK9e0KsXvPsujJlYyLzPTvPXl3M1GdRyQUHw4O/zePIfOWxJczN1dH3+umc6W15fhAkKIua26bS/8WpdVCcAiTEXHvUhIm5gJ3AlkAmsBaYbY7aVqPMLINEYc4eITAOmGmNuEJFuwAdAf6A1sAzobO92wWOWpW/fviY1NfXiP6U6l8djrUeQkQE7d8KmTdbqNJs2QXa2VScmBnPlVZwaNYWNzUazZlMYq1bBt99aSyCLwODBcP31cPPNkJld8WynqvbZme7iH0+H8dXyYOo3MFx3fTY/8/6TAR8/QdCJ4+R27c7JyddycvxkCtp39OmqQa8Qqp6IrDPG9K2wng8JYRDwB2PMGPv1wwDGmCdL1Flq1/lWRIKAg0AUMKtk3bP17N0ueMyy1NmEYIz18HigqIit32chhUVIYQFSkI8UFODKzcGVl4fk5OA+cxrX6VO4T5/CffwY7qxjuLOOEnTwIMGHDhBy+ACuosLiwxfWa0hWdE9+jOrB1oaD+MY1jI1Z0ezcCceO/TeMjh2tJDBqFEyYAM2b/3ebL9Nfq9pr+1YXr78QStLnwRQVCS1CjvFg67eZeOb/iDu6GoD8iMacjE0kt2s8RW2j8bZqiWnVHBo3xNugId769TEhIXRvHwUhIVb7lKoSviYEX/oQ2gD7SrzOBAaUV8cYUyQiJ4FmdvnqUvu2sZ9XdMwqsy2iH21ztlfX4X0mmAs+F8x5Dxde3JzbBNP9It7Ti3CMphwikh9pzX6Gs5827KIjGcSyi47sy2kLOwV2QliYoUVrL81bFjF6nIf2sV46xHrpEu+hWeR/Yz5UAIcyK/c9qNqnS3cvz7ycS052Lqmrg/j26wjeXncPfzpxL43YxxiW0jt7Pb02biBh41vUw7d1s4tw48F93k/GWeU9r42Cjx8hrHH1juLy+05lEZkJzLRfnhGRHQ6EEQkcdeB9fVFBbAbIsh8Vf3V5efD9buuxdlVNxOc4ja/yfIrtJPBapd/CYz8qxZ+/O7jY+JqEX8p7XeZLJV8Swn6gbYnX0XZZWXUy7SajRli/gS60b0XHBMAY8yrwqg9xVhsRSfXlcssJ/hwbaHyXyp/j8+fYQOOrDF8a6dYCnUSkvYiEANOAxaXqLAZusZ9fByQZq3NiMTDNHoXUHugErPHxmEoppWpQhVcIdp/AXcBSwA28aYzZKiKPAanGmMXAG8C7IpIBHMP6BY9d70NgG1AE/K8xxgNQ1jGr/uMppZTylU99CMaYT4FPS5U9WuJ5HnB9Ofs+ATzhyzH9mKNNVhXw59hA47tU/hyfP8cGGt9Fq3DYqVJKqbpBB/oqpZQCNCEgImEiskZENorIVhH5o13e3p6GI8OeliPELi93mo4aju89e+qPLSLypogE2+UjReSkiKTZj0cv/A7VEtvbIrKnRAw97XIRkX/Y390mEeldXbFVEN/XJWL7UUQ+tstr7LsrFadbRDaIyL/t135x7pUTm+PnXQXx+cW5d4H4/OrcO48xpk4/AAHq28+DgRRgIPAhMM0ufwW4037+C+AV+/k0YJ5D8Y23twnW9CBn4xsJ/Nvh7+5t4Loy6o8HPrP3GwikOBFfqToLgJtr+rsrFcP9wPtn39tfzr1yYnP8vKsgPr8498qLz9/OvdKPOn+FYCxn7JfB9sMAo4H5dvkc4Gr7+RT7Nfb2y0UqOa3jJcRnjPnU3mawhvJGV1cMFxvbBXaZArxj77caaCwirZyKT0QaYv0/f1xdMVRERKKBCcDr9mvBT8690rGBNRjE6fPuQvFdQI2eexXF5w/nXlnqfEKA4su6NOAw8AWwCzhhjCmyq5SccuOcaTqwbsRsVpPxGWNSSmwLBm4CPi+xyyC7meQzEbmYmS6qMrYn7Evz2SISapeVNQ1KG6rRhb47rF+0y40xp0qU1dh3Z3sO+DUUz0/SDP8590rHVszp866C+Pzi3LtAfOAf5955NCEAxhiPMaYn1l87/YEuDod0jtLxiUh8ic0vAV8ZY762X68HLjPG9ACep5r/AikntoexvsN+QFPgoeqMoRLxnTUdq9njrBr97kRkInDYGLOuOt+nMnyIzdHz7gLx+cW558P35+i5Vx5NCCUYY04AycAgrEvKs/dplJxao3g6Djl3mo6ajG+s/f6/x5pV9v4SdU6dbSYx1r0ewSISWZOxGWMO2Jfm+cBbWEkWfJsGpdrjA7C/k/7AkhJ1avq7GwJMFpG9wFysJoS/4x/n3nmxici/7Pf2h/OuzPj86Ny70PfnD+demep8QhCRKBFpbD8Px1qjIR3rl8d1drVbgEX28/Km6ajJ+LaLyM+AMVjrSHhL1G95tl1ZRPpj/R9Xyy+NC8TWyi4TrEvjLfYui4Gb7REfA4GTxpgD1RHbheKzN1+H1YmXV6J+jX13AMaYh40x0caYGKxO4iRjzE/wg3OvnNhm+MN5V0F8fnHulRefvdnxc688fj/baQ1oBcwRayEgF/ChMebfIrINmCsijwMbsKbngHKm6XAgviLge+Bb+zz6yBjzGNbJdqe9PRdrtEp1JazyYksSkSisER1pwB12/U+xRntkADnAT6sprgvGZ2+bBjxVqn5NfncX8hD+ce6V5RWcP+8u5D0/OfcuxG/PPb1TWSmlFKBNRkoppWyaEJRSSgGaEJRSStk0ISillAI0ISillLJpQlBKKQVoQlBKKWXThKCUUgqA/w9bcCCJktiyKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "ls = np.linspace(0.9*z[0], 1.1*z[-1], 100)\n",
    "\n",
    "plt.plot(ls, weibull_min(c = cooke['c'], loc = cooke['a'], scale = cooke['b']).pdf(ls), '-b', label = 'Cooke (1979)')\n",
    "plt.plot(ls, weibull_min(c = MLE['c'], loc = MLE['a'], scale = MLE['b']).pdf(ls), '-r', label = 'MLE')\n",
    "plt.hist(z, density=True, alpha=0.2)\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets obtain the optimal solution by enumeration (only for very small $n$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "323.07788987923846"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def enumeration(n):\n",
    "    i = 0\n",
    "    s = [-1]*n\n",
    "    visited = [False]*n\n",
    "    \n",
    "    while i >= 0:\n",
    "        if not s[i] == -1:\n",
    "            visited[s[i]] = False\n",
    "        s[i] += 1\n",
    "        while s[i] < n and visited[s[i]]:\n",
    "            s[i] += 1\n",
    "        \n",
    "        if s[i] == n:\n",
    "            s[i] = -1\n",
    "            i -= 1\n",
    "        else:\n",
    "            visited[s[i]] = True\n",
    "            \n",
    "            if i < n-1:\n",
    "                i += 1\n",
    "            else:\n",
    "                yield cost(s)\n",
    "                \n",
    "opt = min(enumeration(n))\n",
    "opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "P. Cooke. Statistical Inference for Bounds of Random Variables. Biometrika. 1979.\n",
    "\n",
    "B. L. Golden and F. B. Alt. Interval Estimation of a Global Optimum for Large Combinatorial Problems. Naval Research Logistics Quarterly. 1979.\n",
    "\n",
    "R. L. Rardin and R. Uzsoy. Experimental Evaluation of Heuristic Optimization Algorithms: A Tutorial. Journal of Heuristics. 2001.\n",
    "\n",
    "S. H. Zanakis. A Simulation Study of Some Simple Estimators for the Three-parameter Weibull Distribution. Journal of Statistical Computation and Simulation. 1979.\n",
    "\n",
    "S. H. Zanakis and J. R. Evans. Heuristic \"Optimization\": Why, When, and How to Use It. Interfaces. 1981."
   ]
  }
 ],
 "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.5.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}