German_Tank_Problem.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "source": [
        "!pip install lets-plot\n",
        "!pip install pymc3"
      ],
      "metadata": {
        "id": "4ZxmOCKEE2fz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "sns.set()\n",
        "\n",
        "num_tanks = 1000\n",
        "num_captured = 15\n",
        "serial_numbers = np.arange(1, num_tanks + 1)\n",
        "num_simulations = 100_000\n",
        "\n",
        "def capture_tanks(serial_numbers, n):\n",
        "     \"\"\"Capture `n` tanks, uniformly, at random.\"\"\"\n",
        "     return np.random.choice(serial_numbers, n, replace=False)\n",
        "\n",
        "simulations = [\n",
        "    capture_tanks(serial_numbers, num_captured)\n",
        "    for _ in range(num_simulations)\n",
        "]"
      ],
      "metadata": {
        "id": "TmxqozTzFZGD"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The simplest estimate for the number of tanks is simply the maximum number seen. We plot the distribution of this value first and then the average (in red)"
      ],
      "metadata": {
        "id": "DWd_UFrVOeYA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "first_estimates = [max(s) for s in simulations]\n",
        "\n",
        "plt.hist(first_estimates)\n",
        "avg_first_estimates = np.mean(first_estimates)\n",
        "plt.axvline(x=avg_first_estimates, color='red')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "id": "UlKdRPtAFeY_",
        "outputId": "5686b6a3-b9fb-43d8-8a74-37eb42182761"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f5fd98f3e50>"
            ]
          },
          "metadata": {},
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD7CAYAAACfQGjDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAU1ElEQVR4nO3db2zVZ/3/8dc5p6Nlo+X0dG05dEYcUXI24ohtIJlRYjfXBvsH9Uabk2GyOhddxJINsgawB/mjacumM4DbdNktIomZFloIxYSZTGImDRZTS4ZhYGp66J9PKaMwOjjn+t7YPPntt130cHr+wvNxaz3XOZ9zvWnHs/0czqcuY4wRAACfwZ3pDQAAsheRAABYEQkAgBWRAABYEQkAgBWRAABYEQkAgFVepjeQCpcuXVU0mrtv/ygpWSDHmc70NpKOuXILc8Wn8JmnJElXXnsjacdMxFzmcrtdKi6+7zPX4orEzMyMfv7zn+tvf/ub8vPztWLFCu3YsUPnz59XW1ubpqam5PV61dHRoSVLlkhSStbiFY2anI6EpJzfvw1z5RbmisPFi8k/ZoJSsYe4Tjd1dXUpPz9ffX196unpUWtrqyQpFAopGAyqr69PwWBQ7e3tscekYg0AkF6zRuLq1avq7u5Wa2urXC6XJOn++++X4zgaGhpSXV2dJKmurk5DQ0OanJxMyRoAIP1mPd00PDwsr9erPXv26J133tF9992n1tZWFRQUqLy8XB6PR5Lk8XhUVlamcDgsY0zS13w+X6r+DAAAFrNGIhKJaHh4WA899JBeeOEFnT59Wj/84Q/18ssvp2N/CSkpWZDpLcxZaWlhpreQEsyVW5grDvd4kn/MBKViD7NGwu/3Ky8vL3YK6JFHHlFxcbEKCgo0OjqqSCQij8ejSCSisbEx+f1+GWOSvnY7HGc6K15ESlRpaaHGx69kehtJx1y5hbnis/BGRJJ0OcN/VnOZy+12Wb+5nvU1CZ/Pp1WrVunEiROSPvrXR47jaMmSJQoEAurt7ZUk9fb2KhAIyOfzqaSkJOlrAID0c8Xz+ySGh4e1efNmTU1NKS8vTxs2bNDq1at17tw5tbW16f3331dRUZE6Ojr04IMPSlJK1uLFTxLZiblyC3PFZ+HaNZKky91HknbMRKTqJ4m4IpFriER2Yq7cwlzx+f8jUVg0XwX56X+f8oc3Iro8dS2hx94qEnfkO64BIFMK8vNU//zBtD9vz4uNKTku124CAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAVVyRqK6uVm1trRobG9XY2Ki3335bkjQwMKCGhgbV1NSopaVFjuPEHpOKNQBAesX9k8Svf/1rHTx4UAcPHtTXvvY1RaNRbdq0Se3t7err61NVVZV2794tSSlZAwCkX8KnmwYHB5Wfn6+qqipJUnNzs44ePZqyNQBA+uXFe8eNGzfKGKPKyko999xzCofDWrx4cWzd5/MpGo1qamoqJWterzfuoUpKFsR932xVWlqY6S2kBHPlFuaKwz2e5B8zQanYQ1yR2L9/v/x+vz788EPt2rVL27dv1ze/+c2kbyZZHGda0ajJ9DYSVlpaqPHxK5neRtIxV25hrvgsvBGRJF3++JiZjEWic7ndLus313GdbvL7/ZKkefPmKRgM6tSpU/L7/RoZGYndZ3JyUm63W16vNyVrAID0mzUS165d05UrH9XJGKMjR44oEAho+fLlun79uvr7+yVJBw4cUG1trSSlZA0AkH6znm5yHEfr169XJBJRNBrV0qVLFQqF5Ha71dnZqVAopJmZGVVUVKirq0uSUrIGAEg/lzEmd0/eW/CaRHZirtzCXPFZuHaNJOly95HY8eufP5i048er58XGzL0mAQC4OxEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIDVbUViz549WrZsmc6ePStJGhgYUENDg2pqatTS0iLHcWL3TcUaACC94o7Ev/71Lw0MDKiiokKSFI1GtWnTJrW3t6uvr09VVVXavXt3ytYAAOkXVyQ+/PBDbd++Xdu2bYvdNjg4qPz8fFVVVUmSmpubdfTo0ZStAQDSLy+eO7388stqaGjQAw88ELstHA5r8eLFsY99Pp+i0aimpqZSsub1euMeqqRkQdz3zValpYWZ3kJKMFduYa443ONJ/jETlIo9zBqJf/zjHxocHNTGjRuT/uSp4jjTikZNpreRsNLSQo2PX8n0NpKOuXILc8Vn4Y2IJOnyx8fMZCwSncvtdlm/uZ41EidPntS5c+f02GOPSZIuXryo73//+1q3bp1GRkZi95ucnJTb7ZbX65Xf70/6GgAg/WZ9TeKZZ57RX//6Vx0/flzHjx/XokWL9Prrr+vpp5/W9evX1d/fL0k6cOCAamtrJUnLly9P+hoAIP3iek3is7jdbnV2dioUCmlmZkYVFRXq6upK2RoAIP1cxpjcPXlvwWsS2Ym5cgtzxWfh2jWSpMvdR2LHr3/+YNKOH6+eFxtT8poE77gGAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAFZEAAFgRCQCAVVyRePbZZ9XQ0KC1a9cqGAzqzJkzkqTz58+rqalJNTU1ampq0oULF2KPScUaACC94opER0eHDh06pO7ubrW0tGjz5s2SpFAopGAwqL6+PgWDQbW3t8cek4o1AEB6xRWJwsLC2H9PT0/L5XLJcRwNDQ2prq5OklRXV6ehoSFNTk6mZA0AkH558d5xy5YtOnHihIwx+t3vfqdwOKzy8nJ5PB5JksfjUVlZmcLhsIwxSV/z+XzJnh0AMIu4I7Fr1y5JUnd3tzo7O9Xa2pqyTc1VScmCTG9hzkpLC2e/Uw5irtzCXHG4x5P8YyYoFXuIOxL/s3btWrW3t2vRokUaHR1VJBKRx+NRJBLR2NiY/H6/jDFJX7sdjjOtaNTc7mhZo7S0UOPjVzK9jaRjrtzCXPFZeCMiSbr88TEzGYtE53K7XdZvrmd9TeLq1asKh8Oxj48fP66FCxeqpKREgUBAvb29kqTe3l4FAgH5fL6UrAEA0m/WnyQ++OADtba26oMPPpDb7dbChQv1yiuvyOVyadu2bWpra9O+fftUVFSkjo6O2ONSsQYASC+XMSZ3z8tYcLopOzFXbmGu+Cxcu0aSdLn7SOz49c8fTNrx49XzYmNmTjcBAO5eRAIAYEUkAABWRAIAYEUkAABWRAIAYEUkAABWRAIAYEUkAABWRAIAYEUkAABWRAIAYEUkAABWt/1LhwAg2xUWzVdBvv2vtzv1N9OlApEAcMcpyM9L2+W6f35uQpK0+ePn63mxMS3Pmy6cbgIAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWBEJAIAVkQAAWM0aiUuXLukHP/iBampqVF9frx//+MeanJyUJA0MDKihoUE1NTVqaWmR4zixx6ViDQCQXrNGwuVy6emnn1ZfX596enr0uc99Trt371Y0GtWmTZvU3t6uvr4+VVVVaffu3ZKUkjUAQPrNGgmv16tVq1bFPl6xYoVGRkY0ODio/Px8VVVVSZKam5t19OhRSUrJGgAg/W7rNYloNKrf//73qq6uVjgc1uLFi2NrPp9P0WhUU1NTKVkDAKTfbf0+iR07dujee+/Vk08+qT//+c+p2tOclZQsyPQW5uxO/QUmzJVb7tS57lSp+HzFHYmOjg795z//0SuvvCK32y2/36+RkZHY+uTkpNxut7xeb0rWbofjTCsaNbf1mGxSWlqo8fErmd5G0jFXbsnlue7WuCX6+XK7XdZvruM63fTSSy9pcHBQe/fu1bx58yRJy5cv1/Xr19Xf3y9JOnDggGpra1O2BgBIv1l/kvj3v/+tV199VUuWLFFzc7Mk6YEHHtDevXvV2dmpUCikmZkZVVRUqKurS5LkdruTvgYASL9ZI/HFL35R77777meufeUrX1FPT0/a1gAA6cU7rgEAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGBFJAAAVkQCAGA1ayQ6OjpUXV2tZcuW6ezZs7Hbz58/r6amJtXU1KipqUkXLlxI6RoAIP1mjcRjjz2m/fv3q6Ki4hO3h0IhBYNB9fX1KRgMqr29PaVrAID0mzUSVVVV8vv9n7jNcRwNDQ2prq5OklRXV6ehoSFNTk6mZA0AkBl5iTwoHA6rvLxcHo9HkuTxeFRWVqZwOCxjTNLXfD5fMmYFkGaFRfNVkJ/QXzPIEnfkZ6+kZEGmtzBnpaWFmd5CSjBXbknGXPXPH0zCTm5Pz4uNaX/ObJCKr8OEIuH3+zU6OqpIJCKPx6NIJKKxsTH5/X4ZY5K+drscZ1rRqElktKxQWlqo8fErmd5G0jFXbknGXHdqPLNVop8vt9tl/eY6oX8CW1JSokAgoN7eXklSb2+vAoGAfD5fStYAAJkx608SO3fu1LFjxzQxMaGnnnpKXq9Xhw8f1rZt29TW1qZ9+/apqKhIHR0dscekYg0AkH6zRmLr1q3aunXrp25funSp/vCHP3zmY1KxBgBIP95xDQCwIhIAACsiAQCwIhIAACsiAQCwIhIAACsiAQCwIhIAACsiAQCwIhIAACsiAQCwIhIAACsiAQCwIhIAACsiAQCwIhIAAKuEfsc1gNxSWDRfBfm3/787v6MaRAK4CxTk56n++YNpf96eFxvT/pxILk43AQCsiAQAwIpIAACsiAQAwIpIAACsiAQAwIpIAACseJ8EkCaJvqENyCS+YoE0ydQb2iTe1IbEZeXppvPnz6upqUk1NTVqamrShQsXMr0lALgrZWUkQqGQgsGg+vr6FAwG1d7enuktAcBdKetONzmOo6GhIb3xxhuSpLq6Ou3YsUOTk5Py+XxxHcPtdqVyi2lxJ8zwWbJhrgULCpSf5NcG4r0QXlnx/KQ+7+3I1HPf6c97zwOLP/V8mZo50f+/bvU4lzHGJLqhVBgcHNQLL7ygw4cPx25bs2aNurq69PDDD2dwZwBw98nK000AgOyQdZHw+/0aHR1VJBKRJEUiEY2Njcnv92d4ZwBw98m6SJSUlCgQCKi3t1eS1Nvbq0AgEPfrEQCA5Mm61yQk6dy5c2pra9P777+voqIidXR06MEHH8z0tgDgrpOVkQAAZIesO90EAMgeRAIAYEUkAABWRAIAYEUk0qy6ulq1tbVqbGxUY2Oj3n77bUnSwMCAGhoaVFNTo5aWFjmOE3vMrdayxczMjEKhkJ544gnV19frpz/9qaRbX6wx2y/k+N///jf2eWpsbFR1dbVWrlwpKbfnkqS33npLa9euVWNjoxoaGnTs2DFJuT/XX/7yF337299WfX29nnzySQ0PD0vKvbk6OjpUXV2tZcuW6ezZs7HbE51jTjMapNU3vvEN8+67737itkgkYh5//HFz8uRJY4wxe/fuNW1tbbOuZZMdO3aYXbt2mWg0aowxZnx83BhjzLp160x3d7cxxpju7m6zbt262GNutZaNdu7caX72s58ZY3J7rmg0aqqqqmJfh2fOnDErVqwwkUgkp+eampoyK1euNO+9954x5qM9trS0GGNy7/N18uRJMzIy8qm/LxKdYy4zEok0+6xInD592nzrW9+Kfew4jlmxYsWsa9lienraVFZWmunp6U/cPjExYSorK83NmzeNMcbcvHnTVFZWGsdxbrmWjWZmZsyqVavM4OBgzs8VjUbNypUrTX9/vzHGmL///e/miSeeyPm5Tp8+bdasWRP7+NKlS+ZLX/pSTs/1//59kegcc50x664CezfYuHGjjDGqrKzUc889p3A4rMWLF8fWfT6fotGopqambrnm9Xozsf1PGR4eltfr1Z49e/TOO+/ovvvuU2trqwoKClReXi6PxyNJ8ng8KisrUzgcljHGupaN764/fvy4ysvL9fDDD2twcDCn53K5XPrVr36lZ599Vvfee6+uXr2q1157TeFwOKfn+sIXvqCJiQn985//1Je//GX19PRIUs7P9T+JzjHXGXlNIs3279+vQ4cO6c0335QxRtu3b8/0luYsEoloeHhYDz30kP74xz9q48aNWr9+va5du5bprSXNm2++qe9+97uZ3kZS3Lx5U6+++qr27dunt956S7/5zW+0YcOGnP98FRYW6pe//KV+8Ytf6Dvf+Y4cx1FRUVHOz5Vp/CSRZv+7UOG8efMUDAb1ox/9SN/73vc0MjISu8/k5KTcbre8Xq/8fr91LVv4/X7l5eWprq5OkvTII4+ouLhYBQUFsYs1ejyeT1ys0RhjXcs2o6OjOnnypDo7OyV98iKUuTjXmTNnNDY2psrKSklSZWWl5s+fr/z8/JyeS5IeffRRPfroo5KkiYkJvf7666qoqMj5uaTEv+7mOiM/SaTRtWvXdOXKFUmSMUZHjhxRIBDQ8uXLdf36dfX390uSDhw4oNraWkm65Vq28Pl8WrVqlU6cOCHpo39J4TiOlixZYr1YYy5dyPFPf/qTVq9ereLiYkm3vghlLsy1aNEiXbx4Ue+9956kj66V5jiOPv/5z+f0XJI0Pj4uSYpGo3rppZfU3NysioqKnJ9LSvzrbq4zcu2mNBoeHtb69esViUQUjUa1dOlSbd26VWVlZTp16pRCoZBmZmZUUVGhrq4u3X///ZJ0y7VsMTw8rM2bN2tqakp5eXnasGGDVq9efcuLNebKhRxramq0ZcsWff3rX4/dlutzHTp0SL/97W/lcn30G8l+8pOf6PHHH8/5ubZs2aJTp07pxo0b+upXv6rNmzcrPz8/5+bauXOnjh07pomJCRUXF8vr9erw4cMJzzGXGYkEAMCK000AACsiAQCwIhIAACsiAQCwIhIAACsiAQCwIhIAACsiAQCw+j8AD0UirrTH3AAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We see that the average estimate (red line) is less than the total number of tanks. This shows us that the `max(tank_numbers)` is a biased estimator - it will always estimate the number to be less than the actual number of tanks"
      ],
      "metadata": {
        "id": "1WOALXhTOr-D"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "So, it's clear we need to add something to this number. One intuitive solution is to add the average size of the gap between numbers to `max(tank_number)`"
      ],
      "metadata": {
        "id": "MjQxqFBOPEf3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def max_plus_avg_spacing(simulation):\n",
        "    m = max(simulation)\n",
        "    avg_spacing = (m / num_captured) - 1 \n",
        "    return m + avg_spacing\n",
        "\n",
        "new_estimates = [max_plus_avg_spacing(s) for s in simulations]\n",
        "\n",
        "plt.hist(new_estimates)\n",
        "avg_new_estimates = np.mean(new_estimates)\n",
        "plt.axvline(x=avg_new_estimates, color='red')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "id": "ChMy166yF-12",
        "outputId": "6ead151e-97bb-4754-ee45-0aa49cc09c4d"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f5ff0091590>"
            ]
          },
          "metadata": {},
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD7CAYAAACfQGjDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAVDElEQVR4nO3dbWxT9/3+8ct2IKElwXGaBJNOY0UbcotWtEQgddrQ0q6JWG7Y9iCRVSY166qtGgtqQY2AxYybTklot07A2q3VHqEhVd0CCRVhEp3UoakDsTBlQWWiMGWKyc0JUQmUtNjf34P86/+vv/Eljuu7mPfrEfHHh/O9gsNlnxMfu4wxRgAA3IY70wsAAGQvSgIAYEVJAACsKAkAgBUlAQCwoiQAAFaUBADAKi/TC0iFq1evKxrN3Ns/SkoWy3GmMrb/ZMmVHBJZslGu5Ch8+kktXOCRc+C1TC8lYW63S8XF9952FldJTE9P64UXXtBf//pX5efna/Xq1dq9e7cuXbqktrY2TU5Oyuv1qqOjQ8uXL5eklMziFY2ajJbEJ2vIBbmSQyJLNsqJHFeuSAs8uZHlNuI63NTV1aX8/Hz19fWpp6dHra2tkqRQKKRgMKi+vj4Fg0G1t7fHtknFDACQXrOWxPXr19Xd3a3W1la5XC5J0n333SfHcTQ4OKi6ujpJUl1dnQYHBzUxMZGSGQAg/WY93DQ0NCSv16v9+/fr3Xff1b333qvW1lYVFBSovLxcHo9HkuTxeFRWVqZwOCxjTNJnPp8vVd8DAIDFrCURiUQ0NDSkBx98UM8//7zOnTunH/7wh3r55ZfTsb6ElJQszvQSVFpamOklJEWu5JDIko1yIseCmSe1OZHlNmYtCb/fr7y8vNghoIcffljFxcUqKCjQyMiIIpGIPB6PIpGIRkdH5ff7ZYxJ+mwuHGcqoyeRSksLNTZ2LWP7T5ZcySGRJRvlSo4lH0e0cIFnXmdxu13WJ9eznpPw+Xxau3atTp06JWnmt48cx9Hy5csVCATU29srSert7VUgEJDP51NJSUnSZwCA9HPF83kSQ0ND2rZtmyYnJ5WXl6fNmzdr3bp1unjxotra2vTBBx+oqKhIHR0deuCBByQpJbN48UoiOXIlh0SWbJQrOZZsWD/zSuKNnkwvJWF3eiURV0nMN5REcuRKDoks2ShXcvzfkigsWqSC/PS/T/nm9C1d++DDhLa9U0nk5DuuASBTCvLzVP/ckbTvt+fFRqWicrl2EwDAipIAAFhREgAAK0oCAGBFSQAArCgJAIAVJQEAsKIkAABWlAQAwIqSAABYURIAACtKAgBgRUkAAKwoCQCAFSUBALCiJAAAVpQEAMCKkgAAWFESAAArSgIAYEVJAACsKAkAgBUlAQCwoiQAAFaUBADAipIAAFjFVRLV1dWqra1VY2OjGhsb9c4770iS+vv71dDQoJqaGrW0tMhxnNg2qZgBANIr7lcSv/rVr3TkyBEdOXJEX/va1xSNRrV161a1t7err69PVVVV2rdvnySlZAYASL+EDzcNDAwoPz9fVVVVkqTm5mYdP348ZTMAQPrlxXvHLVu2yBijyspKPfvsswqHw1q2bFls7vP5FI1GNTk5mZKZ1+uNO1RJyeK475sqpaWFmV5CUuRKDoks2SgncizwSMqOLKlYQ1wlcejQIfn9fn300Ufau3evdu3apW9+85tJX0yyOM6UolGTsf2XlhZqbOxaxvafLLmSQyJLNsqVHEs+jmjhAk8sSybLItHvp9vtsj65jutwk9/vlyQtXLhQwWBQZ8+eld/v1/DwcOw+ExMTcrvd8nq9KZkBANJv1pK4ceOGrl2baSdjjN566y0FAgGtWrVKN2/e1JkzZyRJhw8fVm1trSSlZAYASL9ZDzc5jqNNmzYpEokoGo1qxYoVCoVCcrvd6uzsVCgU0vT0tCoqKtTV1SVJKZkBANLPZYzJ3MH7FOGcRHLkSg6JLNkoV3Is2bB+5pzEGz2SZnLVP3ck7evoebExc+ckAAB3J0oCAGBFSQAArCgJAIAVJQEAsKIkAABWlAQAwIqSAABYURIAACtKAgBgRUkAAKwoCQCAFSUBALCiJAAAVpQEAMCKkgAAWFESAAArSgIAYEVJAACsKAkAgBUlAQCwoiQAAFaUBADAipIAAFhREgAAK0oCAGA1p5LYv3+/Vq5cqQsXLkiS+vv71dDQoJqaGrW0tMhxnNh9UzEDAKRX3CXxz3/+U/39/aqoqJAkRaNRbd26Ve3t7err61NVVZX27duXshkAIP3iKomPPvpIu3bt0s6dO2O3DQwMKD8/X1VVVZKk5uZmHT9+PGUzAED65cVzp5dfflkNDQ26//77Y7eFw2EtW7Ys9rXP51M0GtXk5GRKZl6vN+5QJSWL475vqpSWFmZ6CUmRKzkksmSjnMixwCMpO7KkYg2zlsTf//53DQwMaMuWLUnfeao4zpSiUZOx/ZeWFmps7FrG9p8suZJDIks2ypUcSz6OaOECTyxLJssi0e+n2+2yPrmetSROnz6tixcv6tFHH5UkXblyRd///ve1ceNGDQ8Px+43MTEht9str9crv9+f9BkAIP1mPSfx9NNP6y9/+YtOnjypkydPaunSpXr99df11FNP6ebNmzpz5owk6fDhw6qtrZUkrVq1KukzAED6xXVO4nbcbrc6OzsVCoU0PT2tiooKdXV1pWwGAEg/lzEmcwfvU4RzEsmRKzkksmSjXMmxZMP6mXMSb/RImslV/9yRtK+j58XGlJyT4B3XAAArSgIAYEVJAACsKAkAgBUlAQCwoiQAAFaUBADAipIAAFhREgAAK0oCAGBFSQAArCgJAIAVJQEAsKIkAABWlAQAwIqSAABYURIAACtKAgBgRUkAAKwoCQCAFSUBALCiJAAAVpQEAMCKkgAAWFESAAArSgIAYBVXSTzzzDNqaGjQhg0bFAwGdf78eUnSpUuX1NTUpJqaGjU1Neny5cuxbVIxAwCkV1wl0dHRoaNHj6q7u1stLS3atm2bJCkUCikYDKqvr0/BYFDt7e2xbVIxAwCkV1wlUVhYGPvz1NSUXC6XHMfR4OCg6urqJEl1dXUaHBzUxMRESmYAgPTLi/eO27dv16lTp2SM0WuvvaZwOKzy8nJ5PB5JksfjUVlZmcLhsIwxSZ/5fL5kZwcAzCLukti7d68kqbu7W52dnWptbU3Zoj6rkpLFmV6CSksLZ7/TPJArOSSyZKOcyLFg5kltNmRJxRriLolPbNiwQe3t7Vq6dKlGRkYUiUTk8XgUiUQ0Ojoqv98vY0zSZ3PhOFOKRs1coyVNaWmhxsauZWz/yZIrOSSyZKNcybHk44gWLvDEsmSyLBL9frrdLuuT61nPSVy/fl3hcDj29cmTJ7VkyRKVlJQoEAiot7dXktTb26tAICCfz5eSGQAg/WZ9JfHhhx+qtbVVH374odxut5YsWaJXXnlFLpdLO3fuVFtbmw4ePKiioiJ1dHTEtkvFDACQXi5jTOaOy6QIh5uSI1dySGTJRrmSY8mG9TOHm97okTSTq/65I2lfR8+LjZk53AQAuHtREgAAK0oCAGBFSQAArCgJAIAVJQEAsKIkAABWlAQAwIqSAABYURIAACtKAgBgRUkAAKwoCQCA1Zw/dAgAsl1h0SIV5Kfpv7cs+mS6VKAkAOScgvy8tF2u+4WL45Kkbf9vfz0vNqZlv+nC4SYAgBUlAQCwoiQAAFaUBADAipIAAFhREgAAK0oCAGBFSQAArCgJAIAVJQEAsKIkAABWs5bE1atX9YMf/EA1NTWqr6/Xj3/8Y01MTEiS+vv71dDQoJqaGrW0tMhxnNh2qZgBANJr1pJwuVx66qmn1NfXp56eHn3uc5/Tvn37FI1GtXXrVrW3t6uvr09VVVXat2+fJKVkBgBIv1lLwuv1au3atbGvV69ereHhYQ0MDCg/P19VVVWSpObmZh0/flySUjIDAKTfnM5JRKNR/f73v1d1dbXC4bCWLVsWm/l8PkWjUU1OTqZkBgBIvzl9nsTu3bt1zz336IknntCf/vSnVK3pMyspWZzpJeTMB5DkSg6JLNkoV3Jki1R8P+MuiY6ODv373//WK6+8IrfbLb/fr+Hh4dh8YmJCbrdbXq83JbO5cJwpRaNmTtskU2lpocbGrmVs/8mSKzkksmSjVOa4W8sn0e+n2+2yPrmO63DTSy+9pIGBAR04cEALFy6UJK1atUo3b97UmTNnJEmHDx9WbW1tymYAgPSb9ZXEv/71L7366qtavny5mpubJUn333+/Dhw4oM7OToVCIU1PT6uiokJdXV2SJLfbnfQZACD9Zi2JL37xi3rvvfduO/vKV76inp6etM0AAOnFO64BAFaUBADAipIAAFhREgAAK0oCAGBFSQAArCgJAIAVJQEAsKIkAABWlAQAwIqSAABYURIAACtKAgBgRUkAAKwoCQCAFSUBALCiJAAAVpQEAMCKkgAAWFESAAArSgIAYEVJAACsKAkAgBUlAQCwoiQAAFaUBADAataS6OjoUHV1tVauXKkLFy7Ebr906ZKamppUU1OjpqYmXb58OaUzAED6zVoSjz76qA4dOqSKiopP3R4KhRQMBtXX16dgMKj29vaUzgAA6TdrSVRVVcnv93/qNsdxNDg4qLq6OklSXV2dBgcHNTExkZIZACAz8hLZKBwOq7y8XB6PR5Lk8XhUVlamcDgsY0zSZz6fLxlZAaRZYdEiFeTb/5spLS1M42qQiIRKItuVlCzO9BJy5sGfKzkksmRK/XNH0r7Pnhcb077PbJCKx0VCJeH3+zUyMqJIJCKPx6NIJKLR0VH5/X4ZY5I+myvHmVI0ahKJlhSlpYUaG7uWsf0nS67kkMiSKfOpzHJBoo8Lt9tlfXKd0K/AlpSUKBAIqLe3V5LU29urQCAgn8+XkhkAIDNmfSWxZ88enThxQuPj43ryySfl9Xp17Ngx7dy5U21tbTp48KCKiorU0dER2yYVMwBA+s1aEjt27NCOHTv+6/YVK1bojTfeuO02qZgBANKPd1wDAKwoCQCAFSUBALCiJAAAVpQEAMCKkgAAWFESAAArSgIAYEVJAACsKAkAgBUlAQCwoiQAAFaUBADAipIAAFhREgAAK0oCAGCV0GdcA5hfCosWqSCfH3fMHY8a4C5QkJ+n+ueOpH2/PS82pn2fSC4ONwEArCgJAIAVJQEAsKIkAABWlAQAwIqSAABYURIAACveJwGkSTLf0FZaWpiUvweYDSUBpEmm3tAm8aY2JC4rDzddunRJTU1NqqmpUVNTky5fvpzpJQHAXSkrSyIUCikYDKqvr0/BYFDt7e2ZXhIA3JWy7nCT4zgaHBzU7373O0lSXV2ddu/erYmJCfl8vrj+Drfblcolzps1JEOu5JD+f5bFiwuUn6GL3ZUVL8rIfjO571zf74L7l/3X/jKVOdGf1ztt5zLGmEQXlAoDAwN6/vnndezYsdht69evV1dXlx566KEMrgwA7j5ZebgJAJAdsq4k/H6/RkZGFIlEJEmRSESjo6Py+/0ZXhkA3H2yriRKSkoUCATU29srSert7VUgEIj7fAQAIHmy7pyEJF28eFFtbW364IMPVFRUpI6ODj3wwAOZXhYA3HWysiQAANkh6w43AQCyByUBALCiJAAAVpQEAMCKkpij6upq1dbWqrGxUY2NjXrnnXckSf39/WpoaFBNTY1aWlrkOE5smzvNMml6elqhUEiPP/646uvr9dOf/lTSnS+wmG0XX/zPf/4T+7dobGxUdXW11qxZM+tasy3HJ95++21t2LBBjY2Namho0IkTJyTNzyx//vOf9e1vf1v19fV64oknNDQ0JCn7s3R0dKi6ulorV67UhQsX4lpbtmf6TAzm5Bvf+IZ57733PnVbJBIxjz32mDl9+rQxxpgDBw6Ytra2WWeZtnv3brN3714TjUaNMcaMjY0ZY4zZuHGj6e7uNsYY093dbTZu3Bjb5k6zbLBnzx7zs5/9zBgz/3JEo1FTVVUVe3ydP3/erF692kQikXmXZXJy0qxZs8a8//77xpiZdbW0tBhjsv/f5fTp02Z4ePi/ftYTXXc2ZPosKIk5ul1JnDt3znzrW9+Kfe04jlm9evWss0yampoylZWVZmpq6lO3j4+Pm8rKSnPr1i1jjDG3bt0ylZWVxnGcO86ywfT0tFm7dq0ZGBiYlzmi0ahZs2aNOXPmjDHGmL/97W/m8ccfn5dZzp07Z9avXx/7+urVq+ZLX/rSvMryv3/WE113tmVKRNZdBXY+2LJli4wxqqys1LPPPqtwOKxly5bF5j6fT9FoVJOTk3eceb3eTCxfkjQ0NCSv16v9+/fr3Xff1b333qvW1lYVFBSovLxcHo9HkuTxeFRWVqZwOCxjjHWWDe+IP3nypMrLy/XQQw9pYGBg3uVwuVz65S9/qWeeeUb33HOPrl+/rt/85jcKh8PzLssXvvAFjY+P6x//+Ie+/OUvq6enR5LmZRYp8XVnc6Z4cU5ijg4dOqSjR4/qzTfflDFGu3btyvSSEhKJRDQ0NKQHH3xQf/jDH7RlyxZt2rRJN27cyPTSEvbmm2/qu9/9bqaXkbBbt27p1Vdf1cGDB/X222/r17/+tTZv3jwv/00KCwv1i1/8Qj//+c/1ne98R47jqKioaF5mudvxSmKOPrnQ4MKFCxUMBvWjH/1I3/ve9zQ8PBy7z8TEhNxut7xer/x+v3WWSX6/X3l5eaqrq5MkPfzwwyouLlZBQUHsAosej+dTF1g0xlhnmTYyMqLTp0+rs7NT0qcvFDlfcpw/f16jo6OqrKyUJFVWVmrRokXKz8+fd1kk6ZFHHtEjjzwiSRofH9frr7+uioqKeZkl0cdTNmeKF68k5uDGjRu6du2aJMkYo7feekuBQECrVq3SzZs3debMGUnS4cOHVVtbK0l3nGWSz+fT2rVrderUKUkzv4HhOI6WL19uvcBiNl988Y9//KPWrVun4uJiSXe+UGS25li6dKmuXLmi999/X9LMNcwcx9HnP//5eZdFksbGxiRJ0WhUL730kpqbm1VRUTEvsyT6eMrmTPHi2k1zMDQ0pE2bNikSiSgajWrFihXasWOHysrKdPbsWYVCIU1PT6uiokJdXV267777JOmOs0waGhrStm3bNDk5qby8PG3evFnr1q274wUWs/XiizU1Ndq+fbu+/vWvx26bjzmOHj2q3/72t3K5Zj4p7Cc/+Ykee+yxeZll+/btOnv2rD7++GN99atf1bZt25Sfn5/1Wfbs2aMTJ05ofHxcxcXF8nq9OnbsWMLrzoZMnwUlAQCw4nATAMCKkgAAWFESAAArSgIAYEVJAACsKAkAgBUlAQCwoiQAAFb/Aw3MPRA3C1+uAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "avg_new_estimates"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Kbml0n9PPiFq",
        "outputId": "3f036fde-a4c3-4cd2-b364-d76f4252b3db"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "999.9626666666666"
            ]
          },
          "metadata": {},
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "As we can see, the new estimate is almost exactly equals to 1000. This new estimate is an unbiased estimate of the actual number of tanks.\n",
        "\n",
        "You can read the derivation of how we can get this estimator and the proof that this is the UMVE from here: https://laurentlessard.com/bookproofs/tank-counting/"
      ],
      "metadata": {
        "id": "K-eh-wH_PnGf"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# also, the new estimator has a higher variance - this is the bias variance tradeoff\n",
        "\n",
        "print(np.std(first_estimates))\n",
        "print(np.std(new_estimates))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VUX2VniZQSlE",
        "outputId": "8cfd6a33-827a-4ae5-8209-588d64106f32"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "58.294202231010935\n",
            "62.180482379744994\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# let's take an example of a sequence of captured tank numbers\n",
        "# and use it to get the estimate using the new estimator\n",
        "captured = [499, 505, 190, 427, 185, 572, 818, 721, \n",
        "            912, 302, 765, 231, 547, 410, 884]\n",
        "\n",
        "print(max_plus_avg_spacing(captured))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1zv0Z7nGGLC0",
        "outputId": "b3fef1be-cc87-4803-84cd-650b6fbca73d"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "971.8\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, we will use PyMC to calculate the poterior probability distribution of the total number of tanks. Let's think of the data generating process:\n",
        "\n",
        "- Let's call the actual maximum numbers of tanks produced `N`. This `N` would be any random number between `max(tanks_captured)` to some large number (in this case, we're taking the max to be 2000). This can be generated using a `DicreteUniform` prior\n",
        "- Then, we get the produce tanks numbered: `1, 2, 3, .... N`\n",
        "- Lastly, we capture some random tanks from this sequence. This can also be modeled using a `DicreteUniform` process with mininum being `1` and maximum being `N`"
      ],
      "metadata": {
        "id": "vY1qtEq1Q1lc"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 75
        },
        "id": "K0BcQzEGEp7p",
        "outputId": "af0fcdc2-b2de-49c8-f435-88546d67cc91"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "    <div>\n",
              "      <progress value='11000' class='' max='11000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [11000/11000 00:01&lt;00:00 Sampling chain 0, 0 divergences]\n",
              "    </div>\n",
              "    "
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "\n",
              "    <div>\n",
              "      <progress value='11000' class='' max='11000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [11000/11000 00:01&lt;00:00 Sampling chain 1, 0 divergences]\n",
              "    </div>\n",
              "    "
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:pymc3:The number of effective samples is smaller than 10% for some parameters.\n"
          ]
        }
      ],
      "source": [
        "import pymc3 as pm\n",
        "\n",
        "with pm.Model():\n",
        "    num_tanks = pm.DiscreteUniform(\n",
        "        \"num_tanks\", \n",
        "\t    lower=max(captured), \n",
        "\t    upper=2000\n",
        "    )\n",
        "    likelihood = pm.DiscreteUniform(\n",
        "        \"observed\", \n",
        "\t    lower=1,\n",
        "\t    upper=num_tanks, \n",
        "\t    observed=captured\n",
        "    )\n",
        "    posterior = pm.sample(10000, return_inferencedata=True, tune=1000)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import arviz as az\n",
        "az.plot_posterior(posterior)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 299
        },
        "id": "jdQNeNmdE3GV",
        "outputId": "8104c63b-dece-4f49-9147-051cddd218d1"
      },
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f5fc57149d0>"
            ]
          },
          "metadata": {},
          "execution_count": 19
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAcgUlEQVR4nO3dd1hUV8IG8BdmGFGKgAMG1KiYVbEjKGtEXUcRVBDsxr4xuomxl6iRxG6iPhZWfWKBxW6+L7ZosC5sYtCgIJZ1o8a1IApRBAxFGGDmfH/4OYaoNAeZ8by/v7j33HvuO3fwneE6xUIIIUBERFKxrOoARET0+rH8iYgkxPInIpIQy5+ISEIsfyIiCbH8iYgkxPInIpIQy5/oNWvSpAmSkpKqOgZJjuVPVAqNRoPTp09XdQwio2L5ExFJiOVPJk+j0SAiIgJBQUHw8vLClClToNVqsW/fPrz33nvFtv39JZXZs2dj/vz5+OCDD+Dp6YkhQ4YgLS0NS5YsQbt27RAQEICff/65xGPPnDkTKSkp+PDDD+Hp6YnNmzcDACZNmoSOHTvCy8sLw4YNw/Xr1w37zJ49GwsWLMC4cePg6emJgQMH4s6dOy+cPyEhAV26dMGZM2cghMDSpUvRoUMHtG3bFkFBQfjll19e5dQRvRTLn8zCkSNHEB4ejujoaFy7dg379u0r835TpkxBXFwcVCoVBg8ejObNmyMuLg7+/v744osvStx/xYoVcHNzw4YNG3D+/HmMHTsWANC5c2ccO3YMP/30E5o1a4YZM2YU2+/w4cOYMGEC4uPj8fbbb2P16tXPzX3y5ElMnz4da9euhY+PD2JjY5GQkIBjx47h3LlzWLNmDRwcHMp4hojKh+VPZmHEiBGoXbs2HBwc0LVrV1y5cqVM+/n5+aFFixaoVq0a/Pz8UK1aNYSEhEChUKBXr15lnuePBgwYAFtbW6hUKkycOBFXr15Fdna2Ybx79+5o1aoVlEol+vTp89xxjh49innz5mHz5s1o1aoVAECpVCI3Nxc3b96EEAKNGjWCi4tLhfIRlUZZ1QGIysLZ2dnwc/Xq1fHgwYMy7VerVi3Dz9bW1lCr1cWWHz9+XO4sOp0Oq1evxtGjR5GRkQFLyyfPoTIzM2FnZwcApR5n69atCA4ORuPGjQ3rOnTogGHDhmHhwoW4d+8eevTogVmzZsHW1rbcGYlKw2f+ZLaqV6+O/Px8w3JaWtprOe6hQ4cQHR2NyMhInDt3DjExMQCA8nw6elhYGKKjo7F169Zi60eOHIl9+/bh8OHDuH37NsLDw42anegplj+ZraZNm+L69eu4cuUKtFot1q5dWynHUavVSE5ONizn5uZCpVLB0dEReXl5WLVqVbnndHFxwZYtW7Bt2zbs2rULAHDp0iVcvHgRhYWFqF69OlQqleGvCiJj428Wma2GDRvi448/xujRo9GjRw94eXlVynHGjRuHr776Ct7e3oiIiEBISAjc3NzQqVMn9O7dG23atKnQvG5ubtiyZQs2b96Mb775Brm5uQgNDUX79u3RtWtXODg4YMyYMUa+NURPWPCbvIiI5MNn/kREEuKrfUh6KSkp6N279wvHoqKi4Obm9poTEVU+XvYhIpIQL/sQEUnIZC77pKVll77Ra+boWAOZmeV/E9DrZA4ZAeY0JnPICJhHTnPICJSc09nZrkJz8pl/CZRKRVVHKJU5ZASY05jMISNgHjnNISNQOTlZ/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRM5h2+uflFZd62mkoBpaVFJaYhInqzmUz5//vmwzJv29JdDaW1yUQnIjI7vOxDRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5l0F2djYWL56HwEA/BAb6ISJiY7Hx69evYfz4D+Dv3wV9+/bCli3hhrHCwkKEhn6CAQOC4OvrjcTEhNcdn4joOSz/Mli7dhXy8/OxZ88hbN68FceOHUZU1EHD+IIFoWjd2hOHD8dg3bpN2L//G8TG/mAYb9WqDT77bBFq1apVFfGJiJ5TZeWv0WgQHh6OoKAgtGnTBv+7NQzZWZnYHDYPcycOxMZVc/E4NwcAkHTzKtZ+OQOhkwdj5cIJuHjhnGGeqKiDGDZsAPz8OmPgwGAcOLDXMJaYmIC+fXth9+4dCAz0Q3Cwf7HSLqtTp05i2LCRsLa2hqurGwIDg4vNk5qagh49ekKhUKBOnbpo2bINbt26CQCwsrLCoEFD0bp1G1haKip6uoiIjEpZlQc/fvw4IiMjodPp0LN3EO7duYlBoybBxbUewv8+H7ExB+Hj2wMRaxfgvfenoUlzL/z36kUsWjAHu3bthaOjIxwdnbB8+Rq4udXBhQuJmDFjEjw8mqNJk6YAgIyMdOTm5uDAgSOIj49DaOgsdOr0F9jb22P79i3YuXPLS/OdO/fsQUYIUeznmzdvGJYHDRqKo0ejMHbsR0hJuYv//OffGDZspLFPFxGR0VRp+Q8fPhxqtRoA4P5Oc9ja10SdtxsBAFp6dsD1Kxdx7sz3aNrCGx4t2wEAGjfzROPGHoiLO4WePQPx7ru+hvk8Pb3Qvv2fcfHieUP5KxRKjB79AZRKJTp08EX16jVw504SWrRoiREjRmPEiNGl5vTxeRc7dmxFaOh8ZGRkICrqILTafMP4u+/6YvHiefj66x3Q6XT461/HwsOjubFOExGR0VVp+T8tfgCwUqlga+fwbNlKhQJtHjLTH+DSuVj8fOmsYUwIHby9nzwY/PTTKURGbkZy8h0IoUd+fj7c3d8xbFuzZk0olc9uprW1NfLyHpcr55QpM7B69QoMGdIP9vY10b27P/75z2MAgKys3zB9+iRMnToTfn4ByMhIR2joLDg6OqFfv4HlOyFERK9JlZZ/WTg4qeH1564YOHKSYV1LdzVsrJUoKChAaOgnCA1dgE6d/gKlUok5c6YXu0RTkm3b/oHt2yNfOn7hwgUAgL19Tcybt9iwfuPG9YZn9ikp92BpaYmePQMBAC4utdG9ew/ExZ1i+RORyTL58vfy6YqwpVNx7T/n8CePNtDpdLh44RzecW8AW1tbFBYWwsHBEQqFAj/9dApnz8ahYcNGZZp75Mj3MXLk+6Vud+/eXdja2sLW1g5nz8bh4MF9WLt2EwCgXr23AQgcP34U3bv3QGZmBqKjT6BtW2/D/gUFBYYHpKKiImi1WqhUKlhYWJT/hBARGYHJl7+DkzNGj/8MUXsjsWPzClhaWqJ5sxb45JNPUaOGDSZPnoHPP5+DwsICdOzYCb6+nY2e4erVK/j731ciJycb9erVx+efL4a7+5MHGBsbWyxZsgJffbUWK1d+gWrVrNGxYyeMGjXGsP/Qof3x66+pAIBp0yYAAL755iBcXd2MnpWIqCwsRFmvkVSyQz9cL/O2Ty/7VDZnZzukpWVX+nFehTlkBJjTmMwhI2AeOc0hI1ByTmdnuwrNyTd5ERFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCQhlj8RkYSUVR2gQiwskJtfZPRpq6kUUFpaGH1eIiJTY5blX1Ckw7WkDKPP29JdDaW1WZ4SIqJy4WUfIiIJsfyJiCTE8icikhDLn4hIQix/IiIJsfyJiCTE8icikhDLn4hIQix/IiIJsfyJiCTE8icikhDLn4hIQix/IiIJvVHlfz81GV+t/BShkwbhi7lj8e/zp1+43fHvdmPGuED88vMFw7p/HduLgf38MXz4INy48V/D+kuXLmDOnOklHjciYiMWLvzsufW+vt64ezfZsE2XLj7w8+sMP7/OGDKkH1atWoaHDx8atk9MTEDfvr3KdZuJiCrijSl/nU6HyPWL0KxVOyxcsxsDRkzAroiVSLt/r9h2Dx+k4lJCLOxrOhnWZT3KwNnYE9iyfR/69u2PjRvXAQCKioqwbt0aTJpUcvmXVbduPXDixEkcORKDpUtXICMjHWPGDC/2AEBE9Dq8MeX/4NdkZP2Wgc7dQ2BpqcCfmrZGw0bNcC4upth2+3d/hd79R0OhfPa5/ZkZaajztjtsbGzg7e2DlJQnDxhbt26Fr29nuLq6GTWrUqmEu3sjLFjwBRwcHPH11zuMOj8RUWnemPJ/EQGBX+8lGZYvJsRCqbSCR8t2xbZTu7gi9d5t5ORkIyHhDBo0cMf9+78iKioK7703otLyKRQKdOrUBRcvnq+0YxARvcgb87VVLrXrwtauJr4/thedu4fgv9cu4eYvl9GoSUsAQH7+Yxw5sBXjpix+bl8bW3t07zUYs2Z8jFq11Jg5cw7CwlZi5syZOHnyX9i/fw9sbW0xbdosuLjUfuHxY2JO4PTpH8udW61WIzs7q9z7ERG9ijem/BVKJUaPD8WB3Rvxr2N7Ubf+O2jt7Qul0goAcPzQLrT9swZO6heXt2f7Lhg5pD9srJU4fToWVlZW8PDwwCefzML27f+D2NiTWL9+DRYs+OKF+2s0fvj880XF1vn6epeaOy0tDXZ29uW8tUREr+aNKX8AcKvbEONnfmlYXvvlDHh36AYA+O+Vi3iUmY6fvo8CAORkZ2H7pi/RNWAANAEDDPtotfnYuHE9Vq78O5KSkuDiUhs2Nrbw8GiG7dsjjZpXr9fj1Kkf4e3d3qjzEhGV5o0q/5S7t+Bcuw6EXo/TPxxG1m+ZaPdudwDA36YtgU6nM2wbtnQq+gz8AE1beBWbY8uWCPTqFQi12hlC2OLOnSRkZKQjMTEBbm51jJKzqKgId+8m4x//2ISMjHQMGTLMKPMSEZXVG1X+5+L+hbOxx6DT6dDwneb429RFUFo9uexjY1v80oqlpSWq17BFNevqhnV37txGfPwZbNz45Bm+i4sLhg8fhREjBsHR0emll3zKKjr6OH788XsIIaBWO8Pb2wcREduhVju/0rxEROVlIYQQVR0CAA79cL3M2zap74RrSRlGz9DSXQ0b62ePh87OdkhLyzb6cYzJHDICzGlM5pARMI+c5pARKDmns7NdheZ8o1/qSUREL8byJyKSEMufiEhCLH8iIgmx/ImIJMTyJyKSEMufiEhCLH8iIgmx/ImIJMTyJyKSEMufiEhCLH8iIgmx/ImIJMTyJyKSEMufiEhCLH8iIgmx/ImIJMTyJyKSEMufiEhCLH8iIgmx/ImIJMTyJyKSEMufiEhCLH8iIgmx/ImIJKSs6gAmxcICuflFz5YzHxdfrqBqKgWUlhavPA8RkbGw/H+noEiHa0kZhmV7u+rIys575XlbuquhtOapJiLTwcs+REQSYvkTEUmI5U9EJCGWPxGRhFj+REQSYvkTEUmI5U9EJCGWPxGRhFj+REQSYvkTEUmI5U9EJCGWPxGRhFj+REQSKrX8b9y4gZEjR8LLywt+fn44ceKEYaygoACTJk2CRqNBkyZNcObMmWL7hoeHIzAwEJ6entBoNAgPDzf+LaAy2bz5K4wcORhduvggImJjsbGHDx9i1qypCA4OgK+vN1JTU4qNR0efwIcfvo9u3TpiwoRxxcbu3EnC7NnTEBjYHT17ajBt2gTcuXO7sm8OEb2iEsu/qKgI48ePR9euXXH27FksXLgQM2fOxK1btwzbtG3bFsuXL4ezs/Nz+wshsGzZMsTHxyM8PBw7d+5EVFSU8W8Flapu3Xr46KNJ6NCh43NjlpYW8PF5F0uWLH/hvvb29hg48D0MHz76ubGcnGz4+nbGrl17cejQcXh4NMfs2dONHZ+IjKzE8r958yYePHiA0aNHQ6FQoEOHDmjbti2+/fZbAIBKpcLo0aPh7e0NS8vnpxo7diyaN28OpVIJd3d3dOvWDYmJiZVzS6hEPXsGokOHjqhRw+a5MSenWujXbyCaNm32wn3btfNBt25+UKuff4Bv1qwFAgNDYG9fE0qlEoMHD8WdO0n47bdHRr8NRGQ85b7mL4TA9evXy30gIQQSEhLwzjvvlHtfMh8XLpxHrVq1ULOmQ1VHIaISlFj+DRs2hJOTE8LDw1FYWIjY2FjEx8cjPz+/3Adau3Yt9Ho9+vfvX+GwZNoePLiPVauWYcKEqVUdhYhKUWL5W1lZYf369fjhhx/g6+uLyMhIBAQEoHbt2uU6yI4dO3DgwAFs2rQJKpXqlQKTacrMzMS0aRPQt+8A+PkFVHUcIipFqV8s27RpU+zYscOwPGTIEISEhJT5AHv27MGmTZuwc+dOvPXWWxVLae7++MXwRlJNpTD6nBWRlZWFadM+RseOnTFq1JiqjkNEZVBq+V+9ehUNGzaEXq/Hrl278ODBA/Tr188wXlBQACEEAKCwsBBarRYqlQoWFhY4ePAgVq9ejW3btqFevXqVdytM3B+/GN5YWrqry7xtUVERdDod9Ho9dDodtFotlEolFIonDyBarRZ6vR4AUFhYAK1Wi2rVqgEAdDrd/+9fBCEEtFotFAoFlEolcnNzMH36BLRs2RoffTTR6LeRiCpHqeX/7bffYs+ePSgqKoKXlxciIyOLXboJCAjAvXv3AABjxjx51hcdHY26detizZo1ePToEQYMGGDYPigoCAsXLjT27aBSLFu2GEeOfGdY3rbtH/j003no1SsIANCt27OXgA4d+uT+io1NAAAcO3YYS5cuMIx369YRPXsGYu7c+Th58ntcufIzbt26WWz+7du/kfcvPSIzYCGePm2vYod+KPsriJrUd6qUZ9J/nNferjqysvOMPq+xtHRXo0E9R6SlZRt9bmNzdrZjTiMxh4yAeeQ0h4xAyTmdne0qNCc/3oGISEIsfyIiCbH8iYgkxPInIpIQy5+ISEIsfyIiCZX6On8yYRYWSMt8bPR3D1dTKaC0tDDqnERkWlj+ZqygSIekWxlGeS/C77V0V0NpzV8NojcZL/sQEUnIZJ7eHdoTiah9W6s6xhuhd79RCBrw16qOQUQmzGTKn0xIZXwKaeZjFOkF/y+ByESw/Ok5lfEppPZ21VHf2Yb/l0BkIkzmX2LQgL+W+VIFP9jt2bypGXlG/w/fSlOJ32vAvyiIysckyl+v10Obl1Pm7XNylOXavqLz5imKoM0r/1dWljavseTkKJH3ON8oGf84r7Hz5imKkPGoEDfvGv+L3ZvWrwUY6cNpsx7nIC+3AACgsrKEwgQfVKysdMjKMv7vk7GZQ05zyAiUnNPKSgd7e3tYWpbv9Tsm8ZHOjx49QlhYWFXHICIyS5MnT4aDg0O59jGJ8tfr9cjKyqrqGM8JCwvD5MmTqzpGicwhI8CcxmQOGQHzyGkOGYHSc1bkmb9JXPaxtLQs96PW62KquX7PHDICzGlM5pARMI+c5pARMH5Oxfz58+cbdcY3TIMGDao6QqnMISPAnMZkDhkB88hpDhkB4+c0ics+RET0evHjHYiIJMTyJyKSkDTlv2zZMmg0GjRp0gS//PKLYf2tW7cwePBg+Pv7Y/Dgwbh9+/Yrjxk748vWV0XGl+XJzMzE2LFj4e/vj6CgIEyYMAEZGc/e2HbhwgX06dMH/v7+eP/995Genl6mMWPnBIDx48ejT58+CAkJwdChQ3HlyhXDmKnc50+tW7fuuTFTOpcajQYBAQEIDg5GcHAwfvzxR5PMqdVqMW/ePPTo0QNBQUH47LPPDGOmcJ/fvXvXcA6Dg4Oh0WjQvn37ys0oJBEfHy9SUlJE165dxbVr1wzrR4wYIQ4cOCCEEOLAgQNixIgRrzxm7IwvW18VGV+WJzMzU8TFxRm2+fLLL8WcOXOEEELodDrRvXt3ER8fL4QQYv369WL27NmljlVGTiGEyMrKMvx84sQJERISYlg2lftcCCEuX74sxowZU2zM1M7li3KbYs5FixaJJUuWCL1eL4QQIi0tzTBmSvf5U4sXLxYLFiyo1IzSlP9Tvz/hDx8+FF5eXqKoqEgIIURRUZHw8vIS6enpFR4zdsaS1ldlxpJyCiHE0aNHxahRo4QQQly8eFH07t3bMJaeni7atGlT6tjryLl//37Rt29fIUTVns8/ZtRqtWLQoEEiOTm52JipncuXnVtTypmTkyO8vLxETk7Oc9uZ0n3+lFarFT4+PuLy5cuVmtEkXudfVVJTU1G7dm0oFAoAgEKhgIuLC1JTUyGEqNCYk5OTyeev7Ix6vR67d++GRqMx5HRzczOMOzk5Qa/X49GjRyWOVebrr+fOnYtTp05BCIHw8HBDTlM5n2FhYejTpw/q1q1bbL0pnssZM2ZACAEvLy9MmzYN9vb2JpUzOTkZDg4OWLduHc6cOQMbGxtMnjwZ3t7eJnWfPxUTE4PatWujefPmACrv91Kaa/70+ixatAg1atTA8OHDqzrKSy1ZsgTff/89pk6diuXLl1d1nGLOnz+Py5cvY+jQoVUdpVQ7d+7EwYMHsXfvXgghsHDhwqqO9BydTofk5GQ0a9YM+/btw4wZMzBx4kTk5JjmZ/rs3bsX/fv3r/TjSF3+rq6uuH//PnQ6HYAnvyQPHjyAq6trhcfMIX9lWrZsGZKSkrBmzRrD281dXV2RkpJi2CYjI8Pwru6Sxl6HkJAQnDlzBpmZmSZzPuPj43Hjxg1069YNGo0Gv/76K8aMGYPY2FiTO5dPb79KpcLQoUORmJhoWG8qOV1dXaFUKhEYGAgAaN26NRwdHXHr1i2Tuc+fun//PuLj4xEUFFQsf2VklLr8a9WqBQ8PD3z33XcAgO+++w4eHh5wcnKq8Jg55K8sq1atwuXLl7F+/XqoVCrD+hYtWiA/Px8JCQkAgK+//hoBAQGljlWG3NxcpKamGpZjYmJQs2ZNODg4mMz5HDduHGJjYxETE4OYmBi89dZbiIiIgK+vr0mdy8ePHyM7OxsAIITA4cOH4eHhUWqW153TyckJPj4+OHXqFIAnr45JT09H/fr1TeY+f2r//v3o0qULHB0dDesqK6M07/BdvHgxjh8/jocPH8LR0REODg6IiorCjRs3MHv2bGRlZcHe3h7Lli2Du7s7AFR4zNgZX7a+KjK+LOeaNWsQGBiIBg0awNraGgBQt25drF+/HgCQmJiIefPmQavVok6dOlixYgXUanWpY8bOuXXrVowfPx55eXmwtLREzZo1MWvWLMP1VVO5z39Po9Fgw4YNaNy4MQDTOZcbNmzAxIkTodPpoNfr0ahRI4SGhsLFxcWkckZFRSE5ORmffvopHj16BKVSiSlTpqBLly4ATOs+9/f3x9y5c9G5c+di+1RGRmnKn4iInpH6sg8RkaxY/kREEmL5ExFJiOVPRCQhlj8RkYRY/kREEmL5ExFJiOVPRCSh/wNrYcBsuNkOngAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We see that the 94% HDI is between 912 to 1112 tanks with the mean being very close to the actual maximum number of tanks produce (1000)"
      ],
      "metadata": {
        "id": "KcGxDtwkTLF5"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "HbTpC43QGdCF"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}