Gaia_exercise_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferring distances with Gaia data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your goal is to infer the distances to stars using Gaia parallaxes. Inverting the parallax will give you an inaccurate distance when the parallax uncertainty is large, so it is better to infer distance. There are two ways you could do this: 1) infer the distance using a sensible prior, as used in [Bailer-Jones et al. (2019)](https://arxiv.org/pdf/1804.10121.pdf), OR infer the prior using HBM.\n",
    "\n",
    "We're going to do both.\n",
    "\n",
    "First, lets open the data file, \"parallaxes.csv\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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>Unnamed: 0</th>\n",
       "      <th>parallax_mas</th>\n",
       "      <th>parallax_error</th>\n",
       "      <th>true_distance_kpc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1.672053</td>\n",
       "      <td>0.018477</td>\n",
       "      <td>0.591963</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.375761</td>\n",
       "      <td>0.035355</td>\n",
       "      <td>2.457926</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>0.529537</td>\n",
       "      <td>0.038041</td>\n",
       "      <td>1.920681</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>1.104286</td>\n",
       "      <td>0.018510</td>\n",
       "      <td>0.901326</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>1.011174</td>\n",
       "      <td>0.021886</td>\n",
       "      <td>0.949794</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0  parallax_mas  parallax_error  true_distance_kpc\n",
       "0           0      1.672053        0.018477           0.591963\n",
       "1           1      0.375761        0.035355           2.457926\n",
       "2           2      0.529537        0.038041           1.920681\n",
       "3           3      1.104286        0.018510           0.901326\n",
       "4           4      1.011174        0.021886           0.949794"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv(\"parallaxes.csv\")\n",
    "df = df[df.parallax_error > 0]\n",
    "df = df.iloc[:100]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could just calculate the distances using 1/parallax. \n",
    "\n",
    "*Plot true distance against 1./parallax.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '1/parallax')"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4lGXWx/HvSaMICgJKh8XGCkoJi7H3VRGxgAJWXBFdwd5WV7GtvaLYKBbUtTdAVNSFFxDDmiBFVFZEIgFUSoBQQsqc94+ZxBBSJpDJTJLf57pyOTPPPfOcjGTOPHc5t7k7IiIiAHHRDkBERGKHkoKIiBRRUhARkSJKCiIiUkRJQUREiigpiIhIESUFEREpoqQgIiJFlBRERKRIQrQDqKzmzZt7x44dox2GiEiNkp6evsbdW1TUrsYlhY4dO5KWlhbtMEREahQzywinnbqPRESkiJKCiIgUUVIQEZEiSgoiIlJESUFERIooKYiISBElBRERKaKkICIS6wryYOZjsCI94qeqcYvXRETqlFXz4cPh8OtC2HYttEmO6OmUFEREYlHeVvi/B+HLJ6FhMzjnFTiwX8RPq6QgIhJrMr6CiSNg7RLofj6c9C9o0LRaTq2kICISK7Zlw+d3wddjoUl7uOB92Oe4ag1BSUFEJBb8+BlMugY2roCUK+DYf0K9RtUehpKCiEg0bVkHn9wCC96AFp3hks+g3V+iFo6SgohIhAx8/isA3rzs0B0PusOi92DKTZCzHo6+GY68HhLqVXOU21NSEBGpbhtXwUfXw+KPoHUP6PchtOwa7agAJQURkYjJzsljY04+6RlZJHdoGrw6mDsBpt4OBdvgxHuC4wfxsfNRHDuRiIjUIukZWfzwazYBh/PGpfL2Oa04aO7t8PMM6HgknDYKmu0T7TB3oKQgIhIBqUvXEnCII8CF/hGd338HEusFk0GPCyEuNqsMRSwpmFl9YAZQL3Sed9z9jhJthgAPAytCD41293GRiklEpLqkdGrGn+N+4f6EMXSPW8r6NifQZMCTsEebaIdWrkheKWwDjnP3TWaWCMwys4/dPbVEuzfdfUQE4xARqV7520he+iyTkh4hm91YevRTdDrmAjCLdmQVilhScHcHNoXuJoZ+PFLnExGJCZlpLHvhYjoGfmFG0jH8u+nfGXfsybv8suVOb61CEe3UMrN4M5sH/A585u5zSmnW38wWmNk7ZtaujNcZZmZpZpa2evXqSIYsIrJzcjfDJ7fi404gqWATQ3JvZGj2MFbmNox2ZJUS0aTg7gXu3h1oC/Q2s5ITcScBHd39YOBz4OUyXmeMu/dy914tWrSIZMgiIpW3dDo8cyikPs23rQfw120PMT3QAwNOPbh1tKOrlGoZ/nb39cB04OQSj691922hu2OByBYKFxGpSluzgnsdTDgd4hJgyBRyT3qYLRa8OkhMiCOlU7MqOVV2Th4r1m8lPSOrSl6vLJGcfdQCyHP39WbWADgBeLBEm1buvip0tx/wfaTiERGpUt9PCq5K3rwGjrg2WKYisQHJwNuXH0bq0rWkdGoWXLS2i0queXhtaEqVvG5pIjn7qBXwspnFE7wiecvdJ5vZ3UCau08ErjKzfkA+sA4YEsF4RER23abfYcqN8N0H0PIgOPctaN19uybJHZpW6Yd24ZoHgLz8AKlL19a8pODuC4AepTw+stjtW4BbIhWDiEiVcYf5b8An/wjuinbc7XD41RCfGPFTp3RqRpxBwKu2S6o0WtEsIlKRrAyYfA389B9olwL9noIW+1fb6ZM7NKVzy8ZszMln1KAeEbtKACUFEZGyBQLBXdA+vyu48OyUh+EvQ6NSoqJx/UQa10+MaEIAJQURkdKtXgwTr4Tlc2Cf4+G0J4JbZNZySgoiIsUV5MGXo+D/HoSk3eDM5+HggVEvURHplcyFlBRERAqt/AY+vBJ+WwgHng59HoFGe0U7qmqlpCAikrcVpt8Ps0fDbi1g4Gvw577RjioqlBREpG5bNis4drBuKfS4AP76L2jQJNpRRY2SgojUTTkb4fM7IO0FaNoRLvwQOh0T5aCiT0lBROqe/30Kk6+F7FVw6Ag49tbgoLIoKYhIHbJ5TXBF8sK3ocWf4ZwJ0LZXtKOKKUoKIhKTqnRTGXf49l34+KZgt9Ext8AR10FC0q6/di2jpCAiNVqFyWPDCvjoOvjfJ9AmGfqNhr0PrMYIaxYlBRGp0bJz8tiYk096Rtb2JSACAZj7EkwdCYF8OOk+OORyiIuPWqw1QfUX8BARqSKF+wxkZm3lvHGpf2xAs/YnePm04GBy6+5wxWw4dLgSQhh0pSAiManMK4BiSu4zMOen30jOnADT7oP4esFqpj0uiHqJippEVwoiEnPKvAIooXCfAYCuCcu5aNFQ+Gwk7HsCDJ8DPS9UQqgkXSmISMwJd6ex5A5N6bp3fc7a9DoXFbyP5TSFs1+CA89QMthJSgoiEnPC3mls+X+ZmHAL5C+GbufCSfdCwz2rN9haRklBRGJOhTuNbdsE/7kH5jwPe7SF898NdhnJLotYUjCz+sAMoF7oPO+4+x0l2tQDJgDJwFpgoLsvi1RMIlJzlLnT2JIvYNI1sGE59L4Ujh8J9RpHJ8gwVOkivGoQyYHmbcBx7t4N6A6cbGYpJdpcAmS5+77A48CDEYxHRGqyLevggyvg1bMgoR5c/DH0eTimEwIEZ1GtWL+1zMHyWBOxKwV3d2BT6G5i6MdLNDsduDN0+x1gtJlZ6LkiUodt9836uw/hoxtgy1o48gY46kZIrB+94MJUOIsq4HDeuFReG5oS8T2Wd1VEp6SaWbyZzQN+Bz5z9zklmrQBlgO4ez6wAdhhRMnMhplZmpmlrV69OpIhi0gsyf4V3jwf3roQGreEYdPh+NtrREKA0mdRxbqIJgV3L3D37kBboLeZdS3RpLQ5YztcJbj7GHfv5e69WrRoEYlQRSRK+oyawREP/mf77hV3+OZVeLo3/G8qnHAnXDoNWh0crTB3SvF1FOXOoooh1bJ4zd3XA9OBk0scygTaAZhZArAHsK46YhKR6Ct1kVrWMnjlDPhwOOx1IPx9NhxxLcTXvMmShbOo2jZtUCO6jiCys49aAHnuvt7MGgAnsONA8kTgIuArYADwH40niNQdxbtXCvLz2TLjKfjlWbA4OPVRSP4bxNXswgtTrj4q2iFUSiRTbyvgZTOLJ3hF8pa7Tzazu4E0d58IjAdeMbMlBK8QBkUwHhGJMYXdK53I5JGksXT/6UfY90To+zg0aRft8OqkSM4+WgD0KOXxkcVu5wBnRyoGEQlPtObSJ7fZjZkpabSc9xSe1AhOHQsHna0SFVFU8zrpRKR2WJEOH15Jm98XQZez4JSHoJEmkkSbkoKIhFWmusrkboHp98NXo6HR3jDodejcJ7LnlLDV7BEcEdll4ZaprhI/z4RnD4PZTwb3ORg+RwkhxigpiNRx1bLAKmcDTLoaXu4bvH/RJOj3JNTfo+rPJbtE3UcidVzYZap30pKZb9Ny5i3slrcWO3QEHPtPSGpYpeeQqqMrBZE6LmILrDavYd3L57PvF0PJzKnP2fn3kN75BiWEGKcrBREpu0z1znCHhW/DxzezR042j+YN4LmCfgQsocwd1CR26EpBRKrOhkzWjz8L3ruUTY068MMZUxgbN4CAJdSY2j91na4UROqI8hao7fKitUAA0l+gYOodJOXmcVf+Bbz56ym80mQ/Xhu6H6lL15LSqZmuEmoAJQUR2TVrlsDEK+GX2axscgiDNw0m0/ciPhCc2TT82H2VDGoQdR+JyM4pyIdZjwfXHfy+CPqN5vcz3mBNQkvireaUipbt6UpBpI6o0lXLqxbAxBGwaj6p9Y8gZfh4aNySZOC1oSnqLqrBlBRE6oAq2xYyLwdmPASznoCGzbh3t3/wceAQRq2rR3Joq+TkDk2VDGowdR+J1AFVsmr5l1R4/kiY+Sh0G8y806cyft3B1VMeQ6qNkoJIHbBL20Juy4YpN8ILJwevFM5/D854mi9XFNS4/YelYuo+EqkDClctb8zJZ9SgHuF37yz5HCZdAxsy4ZDL4LjboV4jIPLlMSQ6lBRE6ohKrVresg4+vRXmvw7N94e/fQrtD9muyU4nGolpSgoiMa5ad0Vzh+8+CHYXbc2Co26Co26AhHqlNq9p+w9LxZQURGJctW2As3EVfHQ9LP4IWnWHCz6All0jdz6JSRFLCmbWDpgAtAQCwBh3H1WizTHAh8DPoYfec/e7IxWTSE1TZVNJKedKwx3mToCpt0PBNjjxbkgZDvH6zlgXRfL/ej5wvbvPNbPGQLqZfebu35VoN9Pd+0YwDpEaq7SppFV6tbDuZ5h0Ffw8AzocDv2egmb7VN3rS40TsaTg7quAVaHb2Wb2PdAGKJkURKQMlZ3hE/b4Q6AA5jwHX9wDcQnQ93HoOQTiNEu9rquW60Mz6wj0AOaUcvhQM5sPrARucPdFpTx/GDAMoH379pELVCTGVHaGT1jjD799FyxRsSId9j8ZTn0M9mgTgeilJor41wIzawS8C1zj7htLHJ4LdHD3bsBTwAelvYa7j3H3Xu7eq0WLFpENWCTGNK6fSJsmDSpMCIXjD2WuMM7P5e1HhpP/7JGQtQz6j4fBbyghyHYimhTMLJFgQnjN3d8redzdN7r7ptDtKUCimTWPZEwiNc2blx0a1nTUcktZZKbB80dx9qZXSa1/JIMSR3HElD1J/2V9hKKWmipiScHMDBgPfO/uj5XRpmWoHWbWOxSP1sqL7IRSS1nkboZPboVxJ8C2jdzRaCTXFQxnzm9xqlkkpYrkmMLhwAXAQjObF3rsVqA9gLs/BwwA/m5m+cBWYJC7ewRjEqm1dhh/KFgAz14V7CpKvphvOl/HKy8sJOC5Rc+JyIwmqdEiOftoFmAVtBkNjI5UDCJ1zZSrj4Kt6+Gz24NrD/bcB4Z8BB2PYPa0JUXdSxD841TNIilJq1NEapPvJwdXJW9eDYdfDcfcAokNgO2ntyYlxDEguS39e7bVVYJsJ6ykYGaXuPv4Yvfjgdvc/a6IRSYiYRn4/FfsUZDFmOZvBusW7X0QnPsGtO6xXTsVsJNwhHulcLyZ9QcuAZoBLwL/F7GoRCQ87vxlw6cM2zKWwJpc4o67DQ6/BuITS21eqUqpUieFlRTc/VwzGwgsBLYAg939y4hGJiLlW/8LG94ewQ1b/o+0wP6MLLiMezr0J7mMhADVVGlVarSwpqSa2X7A1QTXHCwDLjCzhhGMS0TKEgjAf8fCM4fScNV/uTPvQs7JHcni/Fba/Ux2WbjrFCYBt7v7ZcDRwI/A1xGLSkRKt+ZHeKkPTLkB2vXmh7OmMiFwMgHiNJNIqkS4Ywq9C0tUhNYRPGpmEyMXlohspyAPZj8J0x8MziY6/Rnofi4HmdG55XoNHkuVCXdMYaOZdQUOBOoXO/RjRKISkT+snBcsYPfrQjjwdDjlYfq8sJiNU6cxalAP7X4mVSrcMYU7CBasewo4FngI6BfBuEQkbyt8fieMPQ42/Q7nvALnTCB9XVL5he9EdkG4YwoDgOOBX939YqAbUPqmrSKy6zJmw3NHwKzHoftgGD4HDgx+Dyu38J3ILgp3TGGruwfMLN/Mdgd+BzpFMC6RuilnI3xxF3w9Dpq0D+6TvM+x2zWp7MY7IpURblJIM7MmwFggHdgE/DdiUYnURf+bCpOvhY0rIOUKOO42SNpth2ZamSyRZJUtShraRW13d18QiYAq0qtXL09LS4vGqUUiY/Na+PQWWPAmtOgM/UZDu79EOyqpZcws3d17VdSu3CsFM+tZ3jF3n7szwYkI4A6L3oMpN0HOejj6ZjjyekjQcJ1ET0XdR4+Wc8yB46owFpG6Y+PKYDXTxVOgdU84fSLs3SXaUYmUnxTc/djyjovIjvqMmlF2f38gAHNfhs9GBhek/fVfwfGDuPjoBCtSQkXdR2eVd7y0fZdF6qo+o2awetM21mTn4sB541J5bWjKH4lh7U8w6WpYNhM6HgmnjYJm+0Q1ZpGSKuo+Oq2cYw4oKYgA6RlZ/PBr9nY7mxVtddm2Mcx5Fv5zb7Ck9WmjoOdFYOVuTFjlBj7/FaBKqVK+irqPLq6uQERqsuILygolJsRxTJPVMP5SWDkX9j8F+j4Gu7eOTpAiYQh7O04zOxXoQrHaR+5+dznt2wETgJZAABjj7qNKtDFgFNCH4D4NQzSjSWqi4gvKDGjdyHi762xaT3oG6jeBAS9ClzOr/epApLLC3Y7zOaAhwbpH4wiWvaho8Vo+cL27zzWzxkC6mX3m7t8Va3MKsF/o5xDg2dB/RWqU5A5Nefvyw0hdupbjGy2j85xbYd5iOHggnPwANNwz2iGKhCXc2keHufuFQFZoX+ZDgXblPcHdVxV+63f3bOB7oE2JZqcDEzwoFWhiZq0q9RuIxIjkVkkM3zaezh8NgNzNcN47cNaYmEkI2Tl5rFi/VQX0pFzhJoWc0H+3mFlrIA/4U7gnCa2C7gHMKXGoDbC82P1MdkwcmNkwM0szs7TVq1eHe1qR6vPTNHgmBVKfgV5/gyu+gv1OjHZURQoHwlVZVSoS9s5rodpHDwNzCW7J+Xo4TzSzRgS38bymcKOe4odLecoOdTfcfYy793L3Xi1atAgzZJGdM/D5r4pm6lRoaxZ8OBxeOQPiEmHIlOBgcv3dIxtkJamyqoSrwjEFM4sDvnD39cC7ZjYZqO/uG8J4biLBhPBaGWsaMtm+G6otsDKsyEUiJDsnj405+aRnZJVfbO77ScFVyZvXwBHXBstUJDbYpXNHatqoKqtKuCq8UnD3AMXKXbj7tjATggHjge/d/bEymk0ELrSgFGCDu68KL3SRqhdON8uwZz4i9YFT4c3zodFeMGwanHDnLicEiFy/f2Fl1bZNG2y/oE6khHC7j6aaWf/QB324DgcuAI4zs3mhnz5mdrmZXR5qMwVYCiwhWJb7ikq8vkiVK7ebxR3m/ZvHVg+jZ84cOH4kXDoNWnWrknNHut+/cf1E2jRpoIQg5Qp3ncJ1wG5AvpnlEBwLcHcvs+PU3WdR+phB8TYODA8zBpGIK6ubZchjb/P3TaM5JDCPZQkHcm/8FdzQvi/J8YlVdu7SEpI+wKW6hZUU3L1xpAMRiQU7bGDTbnd++fhxnt5wP45xV+BiXs45ngBxO9Y22kWR7vdXeQsJR2VWNDcluMis+IrmGZEISiSaplx9VPDG6sXw4kDaL5/D9EA3bs27hFU0L5oeV9Xf5rWjmsSCcFc0DwWuJjg7aB6QAnyF9lOQ2qggD2Y9Qd60B9hCfRYdfB9/+7oDAYykhDjy8gM4kfk2X5SQRKIk3CuFq4G/AKnufqyZdQbuilxYIlGyYi5MvBJ++5ZPClK4M+8iNn/TlH+d0YWsLbmkdGrG7R8s1Ld5qbXCTQo57p5jZphZPXf/wcwOiGhkItUpdwtMvx++Gg2N9mZKl0e5Mj1YcSU+P0DWllyGH7svEJzF07h+ohKC1ErhJoXM0IrmD4DPzCwLLTKT2uLnmTDpKli3NLjPwYl3s/fvTv0FqeTlB3boJtKArdRm4c4+OjN0804zmwbsAXwSsahEqkPOBvjsDkh/kVVxLXkg6S4u7HYRyQ2akNwBXhuaQurStaR0aqarAqkzKjP7qCdwBMHaRF+6e27EohKJtMWfwEfXQfYqfu1yKcfPPYwtXo9Pi00zLfwRqUvCWtFsZiOBl4FmQHPgRTO7LZKBiUTE5jXwziXw+kCovwdc8jnvNr+cLV4P2H4Vc6UK44nUEuFeKQwGerh7DoCZPUCwWuq/IhWYSJVyh4XvwCc3Q85GOObWYBG7hCRSCrJKXTQWdmE8kVok3NpHyyi2aA2oB/xU5dGI7IQKv9FvyIR/D4T3hvLDtuZckPgI6X8aBglJwB+7pt140gFFXUfaf0DqqnCvFLYBi8zsM4JjCicCs8zsSQB3vypC8YlUqMxv9IEApL8YHEz2Apb3vp2+Mw8g3+P4ukSJipLjB6pDJHVVuEnh/dBPoelVH4pI5RV+ow8429ciWvtTcBFaxpfwp6PgtCeZOL+AfF8MVPxBr/0HpK4qNymY2RjgY+C90D7LIjGl5Df6OT/9RnLmBJh2H8TXg36jocf5YEZKpyzqJ8aVuvagJNUhkrrKgtWryzgY3PjmZOB4IBeYCnzi7vOrJ7wd9erVy9PS0qJ1eokhA5//iuycPL5flY0DBycs5/W9X2W3tQuhc1/o8wjs3mq756RnZGntgdRJZpbu7r0qalfulYK7pwKpBBetNQP+ClxvZgcB3xBMEG9VRcAilZWdk8eaTdtIIpcRCR9wefwkNq1rzPNN/sl1A2+EUvaE0toDkfKFvXjN3dcCrwOvh3ZguxHYN1KBiZSncCyhB4v5d9JY9o1byXsFRzK2waXs3mCvooTQZ9QMdQGJVELYSaE4d3czG+Hu7as6IJGKpGdk8ezU+dwe/zIXxU9lJc24KPdm5sT3oGV8fbKL7XFc6iC0iJSpooHmBWUdAvau+nBEypeekcUz457nrrixtI5fy4SCExnFYE75y36MbL0Ht32wsCgJnNWzraaVilRSRVcKewMnASVX7hgwu7wnmtkLQF/gd3fvWsrxY4APgZ9DD73n7neHEbPUVVvW0XDKVYyPn8RPgVacnTeSH+t14cWLe5PcoSlPT1uyXRIw0LRSkUqqKClMBhq5+7ySB8xsegXPfQkYDUwop81Md+9bweuIwHcfwkc30HnLWp4LnMGovDPwhPq8FkoIsOPagrN6tuWbX7I0piBSCRXNPrqknGPnVvDcGWbWcefCkrpqh4Hh7F/ho+vhh8nQqht2/rtMfGs9zUr5oC9tbYG2txSpnJ0aaK5Ch5rZfIIb9tzg7otKa2Rmw4BhAO3ba2y7piusU1Rys5rtVyd/xSdHLadj+r2Qvw1OuBMOvRLiE4AZZb62koDIrgm3IF4kzAU6uHs34CmCu7qVyt3HuHsvd+/VokWLagtQIiM7J4/0X7Lofe9n2xWau/qNbwg4tLPfGGf30vHLm2CvLnD5l8GKpvEJKlQnEmFRSwruvtHdN4VuTwESzax5tOKR6pGekcX3q7LJL3B+z85l8Ng/PtiTLMDf4j/m06R/0MOW8EvKPTDkI2j+x3KY0grViUjViVr3kZm1BH4LrXnoTTBB6S+8lktdupbihVUKP9jrr/+RRzbdTM/EJUwLdGfDcQ9zxjG9d3i+CtWJRFbEkoKZvQ4cAzQ3s0zgDiARwN2fAwYAfzezfGArMMjLK8QktUJKp2YYFCWGhgkBztz4KnvPeIoNVp+rcofzkR/Gdb5nqc9XoTqRyIpYUnD3wRUcH01wyqrUMS0aJ5FbEGDYPhu4ePUjNPhmMev+dBp//eFU1vju1E8s/wpAg8kikRPt2UdSh6RnZHH2c7NJ8m3clPQOFy/5GGvUEga/wZ4HnMJeo2ZQX1cAIlGlpCARUXzaaeHto/ZvwSG2iPsTx9Ex7je+bdmfrhc9DvX3AHQFIBILojklVWqx7Jw8VoQK02Xn5PHj8pW0+/IWXk+6F4DBuf9kQfc7ixKCiMQGXSlIlUvPyOK7VcGN+s55fjbHks4LiS/QIm89XzQbxIiVJ7GVenwzeREHtGysriKRGKKkIFXuvbmZADRjA3fGv8xp8al8H2jHsLzr+HbVPgRC7VS5VCT2KClIlVu9MYcz4mYxMnECu5HDI3ln83zBaeSRgAEJcYa7a52BSAxSUpCqtX4516+5jQOSUpkb2Jeb8oaxxNsWHa6XGMfIvl3I2pKrfZJFYpCSglSNQADSxrNlyu20J8DdBRfxUt6JBEJzGRLijXN6taN/z7Ykd2hKn1EzeP2/v2j6qUiM0ewj2XVrfoSX+sCUG5gb2IcTtz3Iq4GTOajtHx/2HnDaNGlAcoemKmonEsOUFGTnFeTBzMfg2cPh9+/44oA7OD/3FjJ9LwoCTpc2exBnwabFxw9U1E4kdqn7SCotPSOLJfO/5LRl99Fw3SL4cz/o8whN1iVR/7tU8vIDJCbE0b9nW/r3bEvq0rXbjR+oqJ1I7LKaVoOuV69enpaWFu0w6pTiu6FZfg6pL97EsPjJZNGYjcc9wD5H/7EJX3pG1g5JoDThthORqmFm6e7eq6J2ulKQchXfDe3RcS8xquF4eiYs5+38o7iv4HyGBpIZXqx9coemYX3Ih9tORKqXkoKUqvCb/Ir1W2ngW7kp4Q0uiv+MrEBLLim4hekFB6nrR6QWUlIQYPsuIoCzn5tNwOH4xAVMrTeOVqzllcDJdB34KFck7kZPdf2I1EpKCnVUySRQ2EU04NnZNEiKZ3fP5vbEV+kfP5Pf63Xg/T8/xIE9jqNHh6ZFVU+HH7tveacQkRpISaEOKj5OcN64VM7q2bZoiqjjHJP/JXfVe4kmbOaZwJmkDHyA/p1aFj0/OyePjTn5pGdk6UpBpJbROoU6qOQ6AQMM2Issnk98nGeSnmSVN+OOlk9zyCWP07NYQtDCM5HaTUmhDipcJwDBdQJdWu3OwPhpfF7vRo6Om899eYM5M/duuvY8fIcrAS08E6ndItZ9ZGYvAH2B3929aynHDRgF9AG2AEPcfW6k4pE/JHdoSueWjdmYk8+zffYk571zOTfxW1IDf+bmvEvJ8OCVwbcrN+zwXC08E6ndIjmm8BIwGphQxvFTgP1CP4cAz4b+K1Wo5IBy6tK1vJ22nEBBPq8fNI82HzzCpoBxa/4lvF5wLF7s4tFKeb3iCUXF7ERqn4glBXefYWYdy2lyOjDBg0uqU82siZm1cvdVkYqprik+oHzO819REOr32d+W81DiGNp9/RPzG6Zw2ebz+JXtv/EnJcRxVs+2pb2s9lIWqcWiOfuoDbC82P3M0GM7JAUzGwYMA2jfvn21BFcbFO//Lwg4SeRxRcKHXBH/Idk05Jq8Efzc4GR+ZWPRc+IMBvduz1mhEtciUrdEMymU1jtRaiEmdx8DjIFg7aNIBlWbNG2YVHS7uy3hwcQxHBCXyfsFh3NP3gVsSWzCyN4dWPD+QhyIjzPuOb0r5x6ixCtSV0UzKWQC7YrdbwusjFIstc4xD09j2dotNCCH6xPe5m/xn/ArTbk490ZmWk9jSBqiAAAOaUlEQVQGHtKu6GrggJaNVZxORIDoJoWJwAgze4PgAPMGjSfsuvSMLN6dm8mytVs4LO5bHkgYS/u41bySfwIP5g8iJ2437u63/dWAitOJSKFITkl9HTgGaG5mmcAdQCKAuz8HTCE4HXUJwSmpF0cqlroiPSOL/s/OZnc280DCawxKmM7SQEsGbrudOf5nAOLdydqSG+VIRSRWRXL20eAKjjtsV3VZdkF6RhY3v7uAv8Z9zT2JL9KMjTybfxpP5PdnG8GxBUNrC0SkfKp9VIMVlrdu2jCJJz6YxciEl+ibNIfvAh24JO8GVjY4AKeA+IIA8fFxDEgO7oSmriIRKYuSQg2VnpEVKm/tDEiYxdSkCTQgl4fyBjKm4FTySeC+kzprEFlEKkVJoYZKXbqW1qzmvsRxHBW/kLTA/tycdyk/eRsA/nrg3kWDyUoGIhIuJYUaJj0ji9SfVnPImvcYkvQYAPcEhpB4yDAyvswgDicpIY7Ljt4nypGKSE2kpFCDpGdk8c+x73JP3PP0ivsfK1sczuf7/IM+Bx5EcoemnNi1tbqKRGSXKCnUFAV5/PLB3XwY/wpbSeK63Mtp0Po87j3l4KImWm8gIrtKSSHGpWdksWT+LP764z2cmb2YyYFDuCvvIlbThHOttEohIiI7T0khhs39aRVfv3gjQ+M/Yh27c1netXwa+AsA8Qb9y6hiKiKys5QUYtWyL+n0zuX0TPiFN/OP4f78c9lAIwAS4oy7T++qriIRqXJKCrEmZyN8fiekjWcze3F1wT+ZVdCFxMQ47u3bhawtuRpIFpGIUVKIJf+bCpOvwTeu5MX8U3g4/2wK4hswqHc77W8gItVCSSEWbF4Ln/wDFr4FLTrzbo8XufurYL2i+IDTukkDJQQRqRZKCtHkDt++Cx/fFOw2OvpmOPJ6/rRiC/XTUsnLD6iAnYhUKyWFaNmwgvXvXEWT5Z+zuXk3drvoWdi7CwDJHerx2tAULUQTkWqnpFDdAgGY+zIFn95Ovdxt/Cv/PP7926m8ktOa5GLNtBBNRKJBSaEapGdk8d7cTJpuy+Sy9U/Q+NdUVjXpxfmbzmOZ7018IFjgTklARKJNSSHC0jOyGPjsTC6O/5gRCW+TRwIZh9/Pmv0H8ev4OcRr3EBEYoiSQgSlZ2Tx8vuTeTfpMbrFLeWzgmRuz7uYCxIPZXjHPTVuICIxR0khQh6ePJ+k1Md5NH4iG2w3RuReyeRACknxf1wVaNxARGJNXCRf3MxONrPFZrbEzP5RyvEhZrbazOaFfoZGMp7qMvWTiZzx38FcnfA+EwOHcsK2h5kcOJRubZvw+rBDlQhEJGZF7ErBzOKBp4ETgUzgazOb6O7flWj6pruPiFQc1WrbJvjPvzhhznOssj0ZknsT0wPdAUiKN0ae1kUJQURiWiS7j3oDS9x9KYCZvQGcDpRMCjVeekYWmWkfcfLP91NvUybTdz+DK38/jc00AKBNk/o8ObinEoKIxLxIdh+1AZYXu58Zeqyk/ma2wMzeMbN2pb2QmQ0zszQzS1u9enUkYt1p8/63jKXjLuL0hcNZmV3AD33eZo/+T7AllBAS400JQURqjEgmhdJ2gPES9ycBHd39YOBz4OXSXsjdx7h7L3fv1aJFiyoOcxd8N5H93jmOM+Nm8nR+P/rk3s8XmzuR3KEp7/z9MG486QDe0BiCiNQgkew+ygSKf/NvC6ws3sDd1xa7OxZ4MILxVIn0jCwW/LCYM1Y+QdOMj7E9u3D21htYUNBhu/UGmlkkIjVRJJPC18B+ZvYnYAUwCDi3eAMza+Xuq0J3+wHfRzCeXfbv1AzmTnqa2xJepQG5ZPa6mbZ9buS2zE1abyAitULEkoK755vZCOBTIB54wd0XmdndQJq7TwSuMrN+QD6wDhgSqXh21YfTvqLdtBs5N3Eh/w0cwC15l3JWo2MZHp+oqwIRqTUiunjN3acAU0o8NrLY7VuAWyIZw65K/3kNKz97khNWPIebcVvexbxWcDzxcfEqTSEitY5WNJdj0fyv4d2/c1rcj0wLdOOfeZewkubEGdojWURqJSWF0hTkwawn6Dz9QTZaPa7JvYIPAodjGAlxxt2nd+XcQ9pHO0oRkSqnpFDSirkw8Ur47Vs2/Kkvpy3py6++O0nxxtm9tFeyiNRuSgqFcrfA9Pvhq9HQaG8Y9G/27HwqT2ZkaWaRiNQZSgoAP8+ESVfBuqXQ8yI48W5o0ATQegMRqVvqdFL45scMEv9zF11XvQtNO8KFE6HT0dEOS0QkaupkUkjPyOK7aW9wws8PsRdZvOCn0r3vI/Ts1DraoYmIRFWdSwrzFy9h1WtXckHcbH7wdlyedw3fsi/X/bKFnvtEOzoRkeiqM0khfdk61qW+ymE/PsKfbROP5g3guYJ+5JNAvUTtkSwiAnUkKaRnZPH22Pt4IHEs8wL7cmvgNhYH2hCfEMfg5Lb01zRTERGgjiSF1KVrmRg4DMsL8E7gWM7u3ZFTmzTQNFMRkRLqRFJI6dSMpxIa8Fb+8SQmxOnKQESkDHUiKSR3aMprQ1O0CE1EpAJ1IimAFqGJiIQjkttxiohIDaOkICIiRZQURESkiJKCiIgUUVIQEZEiSgoiIlLE3D3aMVSKma0GMqIdRxmaA2uiHUSYFGtk1KRYoWbFq1h3TQd3b1FRoxqXFGKZmaW5e69oxxEOxRoZNSlWqFnxKtbqoe4jEREpoqQgIiJFlBSq1phoB1AJijUyalKsULPiVazVQGMKIiJSRFcKIiJSREmhkszsZDNbbGZLzOwfpRwfYmarzWxe6GdoNOIMxfKCmf1uZt+WcdzM7MnQ77LAzHpWd4zFYqko1mPMbEOx93VkdcdYLJZ2ZjbNzL43s0VmdnUpbWLivQ0z1lh6b+ub2X/NbH4o3rtKaVPPzN4MvbdzzKxj9Ucadqwx83kQNnfXT5g/QDzwE9AJSALmAweWaDMEGB3tWEOxHAX0BL4t43gf4GPAgBRgTgzHegwwOdrvaSiWVkDP0O3GwP9K+XcQE+9tmLHG0ntrQKPQ7URgDpBSos0VwHOh24OAN2M41pj5PAj3R1cKldMbWOLuS909F3gDOD3KMZXJ3WcA68ppcjowwYNSgSZm1qp6otteGLHGDHdf5e5zQ7ezge+BNiWaxcR7G2asMSP0fm0K3U0M/ZQc+DwdeDl0+x3geDOzagqxSJix1jhKCpXTBlhe7H4mpf+B9Q91GbxjZu2qJ7SdEu7vEysODV2qf2xmXaIdDECo66IHwW+JxcXce1tOrBBD762ZxZvZPOB34DN3L/O9dfd8YAPQrHqjDAojVqg5nweAkkJllfZtpOQ3g0lAR3c/GPicP77RxKJwfp9YMZfgMv1uwFPAB1GOBzNrBLwLXOPuG0seLuUpUXtvK4g1pt5bdy9w9+5AW6C3mXUt0SRm3tswYq1JnweAkkJlZQLFM31bYGXxBu6+1t23he6OBZKrKbadUeHvEyvcfWPhpbq7TwESzax5tOIxs0SCH7Kvuft7pTSJmfe2olhj7b0t5O7rgenAySUOFb23ZpYA7EGUux7LirWGfR4ASgqV9TWwn5n9ycySCA5yTSzeoES/cT+CfbixaiJwYWimTAqwwd1XRTuo0phZy8J+YzPrTfDf7tooxWLAeOB7d3+sjGYx8d6GE2uMvbctzKxJ6HYD4ATghxLNJgIXhW4PAP7joVHd6hROrDXs8wCAhGgHUJO4e76ZjQA+JTgT6QV3X2RmdwNp7j4RuMrM+gH5BL+9DIlWvGb2OsGZJc3NLBO4g+BgGO7+HDCF4CyZJcAW4OLoRBpWrAOAv5tZPrAVGBSND4KQw4ELgIWh/mSAW4H2EHPvbTixxtJ72wp42cziCSant9x9com/sfHAK2a2hODf2KAYjjVmPg/CpRXNIiJSRN1HIiJSRElBRESKKCmIiEgRJQURESmipCAiIkWUFKTGM7NmxapQ/mpmK4rdT4rgeWeZWffQ7U/NrHE5ba8zs/qRikWkqmhKqtQqZnYnsMndHynxuBH89x6ownPNAka4+7ww2mYCXUMrX0Vilq4UpNYys33N7Fsze45gfZ92Zra+2PFBZjYudHtvM3vPzNJCNfJTSnm9hmb2dqi42RtA/WLHMs2siZk1DhWVmx869wAzuxbYC5hpZp+H2o8JnWuRFdu/IPQ6d5rZN6Hz7B96vLGZvWxmC0OPnxF6/BQz+8rM5lpwj4HdIvJmSp2hpCC13YHAeHfvAawop92TwEPu3gs4BxhXSpsRQFaouNmDBCuOltQHWObu3dy9K8HKmY8TrKJ5pLufEGr3j9C5ugEnmtmBxV7jt1C844DrQo/dCax294NCz/k/M9sL+AdwvLv3BBYAO2yiI1IZKnMhtd1P7v51GO1OAA6wP8ryNzWzBu6+tVibo4CHANz9GzNbVMrrLAAeMLMHgEnu/mUZ5xtsZpcQ/BtsTTB5fRc6Vli0Lp1gkimM74zQuR3ICl0tHAjMDsWdBMwK43cVKZOSgtR2m4vdDrB92eXiA78G9A5tnlSecgfh3P17M+tF8MP8YTOb7O73FW9jZvsR/Ebf293Xm9mrJWIprKpZwB9/o1bKuQ34xN0vqCBmkbCp+0jqjNAgc5aZ7WdmccCZxQ5/DgwvvFM4q6iEGcB5oePdgB02ozGzNgQHul8BHiO4xShANsHtMAF2D93fGKqieVIY4U8l2H1VuP9zU2A2cLSZdQo9vlso4YjsNCUFqWtuBj4BviBYl7/QcODw0CDud8ClpTx3NNDMzBYA1wJppbTpBnwdqkh6E1B4lTAG+Dw00DyXYFfRtwRr7JfVxVTcXcDeZvYtMI/g+MRvwCXAm2Y2n2CS2D+M1xIpk6akiohIEV0piIhIESUFEREpoqQgIiJFlBRERKSIkoKIiBRRUhARkSJKCiIiUkRJQUREivw/6G1+Z0tdAR0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "rel_err = df.parallax_error/df.parallax_mas\n",
    "plt.errorbar(df.true_distance_kpc, 1./df.parallax_mas, yerr=rel_err, fmt=\".\")\n",
    "\n",
    "x = np.linspace(min(df.true_distance_kpc), max(df.true_distance_kpc), 100)\n",
    "plt.plot(x, x, zorder=10);\n",
    "\n",
    "plt.xlabel(\"True distance\")\n",
    "plt.ylabel(\"1/parallax\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first try to infer distances using the [Bailer-Jones prior](https://arxiv.org/pdf/1804.10121.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Draw the non-hierarchical and hierarchical PGMs for this problem. In the hierarchical version, L will be a free parameter.*\n",
    "\n",
    "*Using PyMC3, infer the distance to each star with the non-hierarchical model.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [D]\n",
      "Sampling 4 chains: 100%|██████████| 40000/40000 [00:15<00:00, 2637.10draws/s]\n"
     ]
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "\n",
    "L = 500\n",
    "logL = np.log(L)\n",
    "\n",
    "with pm.Model() as simple_model:\n",
    "        \n",
    "    # Prior over distance as a pymc3 potential object\n",
    "    D = pm.Uniform(\"D\", lower=0, upper=1000, shape=len(df),\n",
    "                   testval=np.array(1. / df.parallax_mas))\n",
    "    r = 1e3 * D\n",
    "    pm.Potential(\"lnprior\", -3*logL + 2*pm.math.log(r) - r/L)\n",
    "    \n",
    "    # The model\n",
    "    mu = 1. / D\n",
    "\n",
    "    # Likelihood\n",
    "    pm.Normal(\"parallax_obs\", mu=mu, sigma=np.array(df.parallax_error),\n",
    "              observed=np.array(df.parallax_mas))\n",
    "    \n",
    "    simple_trace = pm.sample(tune=5000, draws=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Measured distance [kpc]')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xuc3HV97/HXey/JBlkkTYLQXI0gClQuWUOoFtHiKVABi1AQxGJF0EIF66n1coqY9uGp1oNFsA0RsCjXcjk2cAALCkYKSdmFAAmRGlMXAlFC3CRLwya72c/5Y34zzM7Ozvw22bntvp+PxzyY38z39/t9MiTzme9dEYGZmRlAU60DMDOz+uGkYGZmOU4KZmaW46RgZmY5TgpmZpbjpGBmZjlOCmZmluOkYGZmOU4KZmaW01LrAEZr+vTpMW/evFqHYWbWULq6ul6JiBnlyjVcUpg3bx6dnZ21DsPMrKFI6k5Tzs1HZmaW46RgZmY5TgpmZpbjpGBmZjlOCmZmluOkYGZmOU4KZmaW46RgZmY5TgpmZhVy5jWPceY1j9U6jFFxUjAzsxwnBTMzy3FSMDOzHCcFMzPLcVIwM7McJwUzM8txUjAzsxwnBTOzBlCtOQ9OCmZmluOkYGZmOQ23R7OZWb1rtKUt8rmmYGZmORVLCpLaJP2HpKckrZH0lSJlzpO0SdKq5HF+peIxM7PyKtl8tAN4X0S8KqkVeETSfRGxoqDcbRFxcQXjMDOzlCpWU4iMV5PD1uQRlbqfmVm96e3r58Utr9HV3VPrUFKraJ+CpGZJq4CXgQciYmWRYh+S9LSkOyTNrmQ8ZmbV0NvXz3+98iprN/ayoec1zrl2RcMkhoomhYjYFRFHALOAhZIOKyhyNzAvIt4BPAjcUOw6ki6Q1Cmpc9OmTZUM2cxsj3R19/CzX/Xycu/OXNNI/8AgK9ZvrmlcaVVl9FFEbAEeBk4oeH1zROxIDr8DLBjh/KUR0RERHTNmzKhorGZme2LF+s0MFjSUt7Y0sWj+tNoENEqVHH00Q9K+yfMpwPHAzwrKHJB3eAqwtlLxmJlVw6L502hS5rmAGe2TuOn8RSyYO7WmcaVVydFHBwA3SGomk3z+JSLukbQY6IyIZcCnJZ0CDAC/Ac6rYDxmZhW3YO5U3rZ/O9v6BtinrYX2ttaGSQhQwaQQEU8DRxZ5/bK8518AvlCpGMzMaqG9rZX2ttZah7FbPKPZzKwBVGt4q5OCmVmdy45oqsbwVicFM6tL1do/oBHkj2iq9PBWJwUzszF224XHcNuFx4zZ9fJHNFV6eKuXzjYzq3P5I5quPOvIio5mclIwM2sA2RFNlR7e6uYjM2to7nsYW04KZmaW46RgZmY5TgpmZpbjpGBmZjkefWRmViFjOVehWlxTMDOzHCcFM2tojbgPcj1zUjCzhlXNheJ2V6PNo3BSMLOGVc2F4iYKJwUzq0tpmoWquVDcRDHi6CNJf5Hi/P+OiGvGMB4zs1yz0GDAOdeuGHGP42ouFFdr1RrJVKqm8JfA3kB7icdnKx2gmU08o2kWam9rZea+U8Z1QqimUvMUvh8Ri0udLOkNJd5rA5YDk5P73BERXy4oMxn4HrAA2AycGRG/TBe6mY1X2WahwXCzULWNWFOIiM+VO7lMmR3A+yLicOAI4ARJiwrKfBzoiYgDgW8CXysfspmNtXobIZNtFpo1dcqITUdWGWU7miV9VdK+ecdTJf1tufMi49XksDV5REGxU4Ebkud3AL8vSakiN7NxYaSENF6ahRptHkWa0UcnRsSW7EFE9AAnpbm4pGZJq4CXgQciYmVBkZnAC8l1B4CtgOuJZjYuNMI8ikJpkkJz0vYPgKQpZPoJyoqIXRFxBDALWCjpsIIixWoFhbUJJF0gqVNS56ZNm9Lc2swaRKP9kh6NRpxHkSYp3Aj8SNLHJf0p8ACvN/mkktQ0HgZOKHhrAzAbQFIL8EbgN0XOXxoRHRHRMWPGjNHc2szq2Fj8kr7twmPqduG5RpxHUTYpRMTXgb8F3g4cAvxN8lpJkmZk+yKS2sXxwM8Kii0D/iR5fjrw44gYVlMws/GpEX9Jj0YjdpinXTp7FZkmo0HgqZTnHADcIKmZTPL5l4i4R9JioDMilgHXAd+XtI5MDeGsUUVvZg1tIgw9bW9rpb2ttSESAqRICpLOBy4DfkymD+AqSYsj4vpS50XE08CRRV6/LO95H3DGaIM2s/FhIs1IbhRpagp/CRwZEZsBJE0DHgVKJgUzqy/ZYZ/11v7eaL+kx7s0Hc0bgN68416SYaRmZja+pKkpvAislPSvZIaLngr8R3bBvIi4ooLxmVkV9Pb1s61vgK7unqr/Yq+3mstEl6am8AvgB7w+f+BfgY28viiemTWwRpxgZZWTpqZwQ+EidZLeGRGPVyYkM6umYsNC3b4/dhqtJpSmpnCnpJnZA0nH4k5ms3Gj0hOsurp7+PZD61wDaRBpagqfBH4g6WTgKOCrpFz7yMzq354MCy02oin/ta7uHs5Y8iiDAW2tTQ0zgWsiK5sUIuJxSZ8G/g3oA94fEV6AyGwcGcthofmd1nvSNNVozS7jRantOO9m6OJ0e5FZxfQ6SUTEKZUOzswaS+E2mpd94FDaWpvoHxgctzOWx5tSNYVvVC0KMxsXCmsGPdt3ctP5i1ixfjOL5k9z01EDGDEpRMRPqhmImTW+YmsZLZg71cmggYw4+kjSPeVOTlPGzPZMvW2VWUojrgpqQ5VqPnq3pGUl3heZpbTNrAFUYtZysWvee8mxY3Jtq41SSeHUFOfvHKtAzKxyCjuAx+JXfCWuabXnPgWzCaASs5Y9E3p8SjOj2cwaXCVmLTfiVpNWXtqd18ysgVViMxtvkDM+pUoKyR7LcyLiuQrHY2YVUmrW8u7OHvYGOeNP2eajZM2jVcD9yfERZUYlmVmNlBu++uzGbQ0zvNVqI02fwuXAQmALQESsAuaVO0nSbEkPSVoraY2kS4qUOU7SVkmrksdlxa5lNpH19vXz4pbXvMqoVUWapDAQEVt349oDwGcj4u3AIuAiScXmNfw0Io5IHot34z5m45Y3wLFqS5MUVks6G2iWdJCkq4BHy50UERsj4onkeS+wFphZ+iwzy1ds2OdYaKRZ0lZdaZLCnwOHAjuAm8mslHrpaG4iaR5wJLCyyNvHSHpK0n2SDh3h/AskdUrq3LTJq3bbxDHaYZ/Vbmq67cJjvMT1OFM2KUTE9oj4UkS8M3n8r4joS3sDSXsDdwKXRsS2grefAOZGxOHAVWT2gi4Ww9KI6IiIjhkzZqS9tVnDG81aQm5qsrGQZvTRA5L2zTueKumHaS4uqZVMQrgpIu4qfD8itkXEq8nze4FWSdNTR282AbS3tTJz3yllh32OpqnJndc2kjTNR9MjYkv2ICJ6gP3KnSRJwHXA2oi4YoQy+yflkLQwiWdsGk3NJohs/0CapqZdg8F/vfIqaze6RmHFpZm8NihpTkQ8DyBpLkN3ZBvJu4BzgWckrUpe+yIwByAilgCnA5+SNAC8BpwVEWmubWYFys0w7u3rZ/vOXWzfuSv3mtcsskJpksKXgEckZRfIOxa4oNxJEfEImeW1S5W5Grg6RQxmlkKpGcZ/+I7f5tmNQxcl8JpFVqhsUoiI+yUdRWaugYDPRMQrFY/MzMZU/q5oAqa3T2LJRzpcS7Ah0i6INxn4TVL+EElExPLKhWVmYy2/eWmfthavWWRFlU0Kkr4GnAmsAQaTlwNwUjCrgnLzAPJ3Pysn27zkuQU2kjQ1hQ8CB0fEjkoHY2bpZGcjf+6Et3n3MxtTaYakrgdaKx2ImY1epZbBsIkrTU1hO7BK0o/ILHUBQER8umJRmVkq+Z3HHklkYyFNUliWPMysznj3MxtraYak3lCNQMwsvfzO5dHsfuYOZisnzeijg4D/DRwCtGVfj4j5FYzLzEaQXfgu27k8f/obaG9zt5+NjTQdzd8F/onMpjnvBb4HfL+SQZnZyAo7l7f1DdQ2IBtX0iSFKRHxI0AR0R0RlwPvq2xYZjaSwoXv9mlLOwfVrLw0f5v6JDUBP5d0MfAiKVZJNZuosnMIKtV+v7udy5WOy8aHNDWFS4G9gE8DC4CPAB+tZFBmVlraPRbMRitNTWFeRDwOvAp8DEDSGRTfWtNswvKexzYepKkpfCHla2Zm1uBGrClIOhE4CZgp6Vt5b+1DZiSSmRWRP4egUs077hewSilVU3gJ6AT6gK68xzLgDyofmllj6e3rZ723urQGN2JNISKeAp6SdHNE9ANImgrMTvZpNrNE/oSyLG91aY0oTZ/CA5L2kfRbwFPAdyVdUe4kSbMlPSRpraQ1ki4pUkaSviVpnaSnkx3ezBpO/oSyLC9QZ40oTVJ4Y0RsA04DvhsRC4DjU5w3AHw2It5OZivPiyQdUlDmROCg5HEBmZnTZg0nf0KZgP3aJ9Xd3ga9ff28uOU1N2lZSWmSQoukA4A/Bu5Je+GI2BgRTyTPe4G1wMyCYqcC34uMFcC+yb3MaubMax4rOby02PvZCWWzpk7h7Qe08+bpe+9WQih3792Vbd5yX4eVkyYpLAZ+CKyLiMclzQd+PpqbSJoHHMnwuQ0zgRfyjjcwPHGYNYTshLJ6XJzOm/FYWmmWzr4duD3veD3wobQ3kLQ3cCdwadIMNeTtYrcsco0LyDQvMWfOnLS3Nquqeh4m6s14LK1S8xQ+FxFfl3QVRb6o0+y8JqmVTEK4KSLuKlJkAzA773gWmaGwhfdaCiwF6OjoGBaL2ViqxjyDavNmPJZWqeajtcl/Oxk6TyH7KEmSgOuAtREx0milZcBHk1FIi4CtEbExbfBmY63Wbe+V7Az2ekmWRql5Cncn/93dndfeBZwLPCNpVfLaF4E5yXWXAPeSmTW9jsxe0B/bzXuZjYlibe8L5k6tygqjhZvn1NvoJZsYSjUf3U2RZqOsiDil1IUj4hGK9xnklwngojIxmlXNSG3v+U1KlWpeGikhmVVTqY7mbyT/PQ3YH7gxOf4w8MsKxmRWM8Xa3vN/wX/4OyvoHxgkGPtf8+4MtnpQqvnoJwCS/iYijs17625JyysemVmNtLe10t7WmvuyL/wFn60+j/WveXcGWz1Is5/CDEnzk6GoSHozMKOyYZnV3klXLmdb3wB/dtyBQ37BZxNDJX7NFyYks2pLkxQ+AzwsaX1yPI9kzoDZeNXb159rMlp8zxrm/NZeDAwGV551JH/9g2f8a97GrTST1+6XdBDwtuSln0XEjsqGZVZd+aOLbrvwGL790Dqe3fgckGkmOqNjNhe990DAv+ZtfEtTUyBJAk9VOBazulGrTt9KDnmt5xnXVj/SrH1kNu4VThrLX+DO8wVsInFSsAlvpFnMngFsE1GpyWslN7zJLottVmt7Ott4tJPG3Axj41mpPoX/k/y3Degg06cg4B1klsB+d2VDM6uONLOYXVuwiWLE5qOIeG9EvBfoBo6KiI5k17UjyaxVZNZwTrpyOe/+2o+HLDhXrP+g1gvjmdVKmj6Ft0XEM9mDiFgNHFG5kMxGJ+3KoqW+6Av7D7wpjU1UaYakrpV0LZm1jwL4CK8vq21WM2de89iQSWbl1iIq1XdQ2E/gdYhsokpTU/gYsAa4BLgUeBYvcW11YlvfwJAv+ktufXLEPY6zX/RQ/oveQ1Jtokozo7lP0hLg3oh4rgoxmaW2T1sLIlOFbW5uYp+2oX+l80cmLZg7lds/+busWL+ZRfOnlf2i98xlm4jKJgVJpwB/D0wC3izpCGBxuf0UzCqtt6+fV17d8fqmH1F+p9YFc6f6S96shDR9Cl8GFgIPA0TEKknzKheSWXn5exxk7RoMNr26Y8gw0j0ZVur5CDYRpelTGIiIrRWPxGwU8juNs5qbxCu9O3Oji25e+byHlZqNUpqksFrS2UCzpIMkXQU8WuG4zEbU1d3DS1tey+31OqmlibOPnsMZHbOHbIBz3+qNHlZqNkppksKfA4cCO4Cbga1kRiGVJOl6SS9LWj3C+8dJ2ippVfK4bDSB28TU1d3DGUse5aaVzwMwo30St3xiEV/9o9/htKNmDRlddOJhB6QebWRmGSX7FCQ1A1+JiL8EvjTKa/8zcDXwvRJlfhoRHxjlda3BlVqrqNw6RvnNRgFMbmnO9RUU287yxhW/9IY4ZqNQsqYQEbuABbtz4YhYDvxmd861iavc7ORF86flmo0Ew4agFs5M9kqnZqOTpvnoSUnLJJ0r6bTsY4zuf4ykpyTdJ+nQMbqmNYBnXtw6bA2ik65czrMb03cMtzaLv/ng75Qsk91JzczSSZMUfgvYDLwPODl5jEWTzxPA3Ig4HLgK+MFIBSVdIKlTUuemTZvG4NZWS7/e1sf2nbuGfflv6xvIlRmpY3jF+s25zuRdg+HOY7MxlmZGc0WWtIiIbXnP75X0j5KmR8QrRcouBZYCdHR0lJ+hZHWrq7uH7s3bc8c789Ygasn2CjNyx3C5NYlcKzDbM2lmNH8XGPZFHBF/uic3lrQ/8OuICEkLydRa/LNvnMv/pQ/QJLFo/jS6unt4/jeZZCHgsg8cWrQfoFhnspmNnTQzmu/Je94G/BHwUrmTJN0CHAdMl7SBzMzoVoCIWAKcDnxK0gDwGnBWRIp1Cqyh5f/SF7D41MNYMHcq335oXW5UUZOgZ/vOEa/hNYnMKidN89Gd+cfJl/2DKc77cJn3ryYzZNUmmOl7TyKAJR/pyH2xe6lqs/qQpqZQ6CBgzlgHYuNL/nyD7PPPnfA2zljyKIMBed0HwOiahdxvYFY5afoUehnap/Ar4K8qFpGNS89u3MYltz6ZayIaDIZscgNuFjKrB2maj9qrEYiNL8+8uJWBwaCru4fevn76+nfR3ERu7wMBU/eaNOQc1wDMaq/sPAVJ75L0huT5RyRdIWlu5UOzRtXV3cP2nbvYOTDIGUse5dmNvfTvCrZsH8g1GwWw+J41XrnUrM6kmbz2T8B2SYcDnwO6Kb2ekU1w1/zkF7nnhctb78o79sqlZvUn7X4KAZwKXBkRVwJuUrIR/Xpb34jvtSRNSOBRRmb1KM3oo15JXwA+AhybrJzaWtmwrJGd+c45PLXhmWGv7z25mRv+9GiA1Pskm1l1pakpnElmL4WPR8SvgJlk9mw2G+akK5fzjw+vY/99JjO5Zehfr08dd2Buj+T/9/RLXHLrk+5TMKszZZNCRPwqIq6IiJ8mx89HhPsUbJjsvskbel5jy2v9fPnkQ3Mdy22trzcV5ZfzNplm9SXN6KNFkh6X9KqknZJ2SfKezTZM/gY4/QOD9Gzfydv2b2fW1CncdP6iXFNRYTl3NpvVjzR9ClcDZwG3Ax3AR8nMarYJonA3tK7unqJ9AovmT6OttYn+gcFcJ/JF7z1w2PW8pIVZ/Uq1zEVErJPUnOzE9l1Jj1Y4LqsjvX39bOsbyDXzZJeqaGttGlIDWDB3Kjedv6hsJ/KCuVO5/ZO/685mszqUJilslzQJWCXp68BG4A2VDcvqRbb9fzDgnGtXcNpRs4Y1/eR/qWc7kstJW87MqitNUjiXTN/DxcBngNnAhyoZlNVetonoxS2vDUkCgmFNRGY2fqRZ+6hb0hTggIj4ShVisioo7Cc46crlbOsb4M+OO5DVL23llpXPE8CklqbcekWtLU2cdtQsnny+x5vcmI1TaVZJPRn4BjAJeLOkI4DFEXFKpYOzyinsJ8g2EX3x/w6ddLZr1yBnHz2H3953Sq7936uZmo1faZqPLgcWAg8DRMQqSfMqFpFVXKl+gnzi9dpBfgLITyhODGbjS9q1jzwvYRwpnCcgXl+PKN/7D3nTkNFF4IlnZuNdmqSwWtLZQLOkgyRdBXhIagPLzhOATE2gfXILxTbHnt4+eVhNwBPPzMa3NEnhz4FDyax/dAuwDbi03EmSrpf0sqTVI7wvSd+StE7S05KOGk3gtmfOWjiHc46ew2UfOJQly9cXLVOs9lCYUDz6yGx8STP6aDvwpeQxGv9MZjb0SOsknUhmZvRBwNFk9m04epT3sDKyo4quPOtIAO58YkNuZFFb68i/CSYlfQmFRrOXspk1nhGTgqRlpU4sN/ooIpaX6ZA+FfheslfDCkn7SjogIjaWuq6ll9+h/MfXPMaugt7k/oFBXu7dMey83ztoOpce/9YRv/DvveTYisRrZrVXqqZwDPACmSajlRRvTdgTM5PrZ21IXhuWFCRdAFwAMGfOnDEOY/zKb/8vTAiQaf7Zr33ykNdamlQyIZjZ+FaqT2F/4IvAYcCVwPuBVyLiJxHxkzG4d7EkU6y/k4hYGhEdEdExY8aMMbj1+NfV3cN3fvqLEd8/fNYbuen8RZx21KxcH0FLk1h86mFOCGYT2Ig1hWTxu/uB+yVNBj4MPCxpcURcNQb33kBmyYysWcBLY3DdCa+ru4fT/+nR4hkWaGkWl518aO7L34vTmVlWyY7mJBn8IZmEMA/4FnDXGN17GXCxpFvJdDBvdX/Cnrt55fMsXf6LEROCgDM7Zu/WInZmNv6V6mi+gUzT0X3AVyKi6NDSEuffAhwHTJe0Afgyyd7OEbEEuBc4CVgHbAc+thvxW56bVz4/bJmKfAImtxYfVWRmBqDM4J8ib0iDwH8nh/mFBERE7FPh2Irq6OiIzs7OWty6rp105XJ+/vKr9O96/X9V++RmenfsAjJDTE9fMIsPFSxZYWYTg6SuiOgoV65Un0KaiW1WQ9nlrafuNSk39DTfMW+ZzoXveYv7C8wstVQ7r1n96eruye2A1tKkYQmhuQkufM9b3F9gZqPi2kCDyp+DMDgYuWGlkGnfO+udc5wMzGzUnBQazM0rn+fc61bS+1o/ba1NNAsmtTZxwe/Np6VJNMmdyWa2+9x81EDyRxf99Oev8Mlj59M+pTXXX/D+Q/d3/4GZ7REnhQZy2+PPDzles3Eb3//462sIuv/AzPaUk0Kdyx9h9NSGoXsdnXjYATWKyszGKyeFOlY4wijf/zjkTZx9tBcHNLOx5aRQp7q6e/iHB/9zyAijliYREbS2NHHhe95S2wDNbFxyUqhD+TUEgCZlZiRf9oFD6dm+0x3JZlYxTgp1KH8OQhPwrgNLb3pjZjZWPE+hDi2aP23IHAQnBDOrFtcUaig7sqiwOWjB3KncdP4izzkws6pzUqiRru4ezrl2BTsHBpnU0sRN5y/yHgdmVnNOClXQ1d3DXU9sICC3dPWK9Zvp6x8EoH9gkBXrNzsJmFnNOSlUWFd3Dx/+TqZGAHBH5wvccsExuX6D/oFBWluaWDR/Wo0jNTNzUqioru4eFt+9JpcQAPp3BSvWb+ai9x7ofgMzqztOChVy88rn+et/Xc2ugo0OWpuVqxW438DM6k1Fh6RKOkHSc5LWSfp8kffPk7RJ0qrkcX4l46mWru6eognh8Flv5JYLjnEiMLO6VbGagqRm4NvA+4ENwOOSlkXEswVFb4uIiysVRy3c9cSGYQlhUrO47ORDnRDMrK5VsqawEFgXEesjYidwK3BqBe9XN17u3THkeOa+ba4hmFlDqGRSmAm8kHe8IXmt0IckPS3pDkmzKxhPxXR19/Dth9bR1d0DwH7tk4e8f9zB+zkhmFlDqGRHs4q8VrC9PHcDt0TEDkmfBG4A3jfsQtIFwAUAc+bU13LR+YvXtbVmJqGddtQsbu/akBtu6q0xzaxRVLKmsAHI/+U/C3gpv0BEbI6IbFvLd4AFxS4UEUsjoiMiOmbMmFGRYNMqrBWsWL85917+JLRbPrGI//kHB3PLJxa5lmBmDaOSNYXHgYMkvRl4ETgLODu/gKQDImJjcngKsLaC8eyxYktTLJo/jUktwyehebipmTWiiiWFiBiQdDHwQ6AZuD4i1khaDHRGxDLg05JOAQaA3wDnVSqePZWdiFa4NIUnoZnZeKKIwmb++tbR0RGdnZ1Vu19Xdw93PrGB2ztfoH/X65/VpGZ5RJGZNQxJXRHRUa6cZzSXkG0u2tE/OKyH/IyO2U4IZjbueJOdElas38zOgeEJYVKzPKLIzMYl1xRKyO9Ebm4Sxx28H9PbJ+eWvzYzG2+cFErwDmhmNtE4KZThoaVmNpFM6D6FwoloZmYT3YSsKWS3x7y98wUGBqPoHslmZhPRhEsKxYaZeo9kM7OMCZMUurp7WLF+My9teS03Kxkyq/Z5j2Qzs4wJkRTy1yxqaRKTWprYtWuQ5uYmTl8wy0NMzcwSEyIpZCehDQbsGgzOXDibmftO8TBTM7MCEyIpFK5k6pqBmVlxEyIpeBKamVk6EyIpgCehmZmlMaEnr5mZ2VBOCmZmluOkYGZmOU4KZmaW46RgZmY5TgpmZpajiMLNJuubpE1Ad63jGMF04JVaB5GSY62MRooVGitex7pn5kbEjHKFGi4p1DNJnRHRUes40nCsldFIsUJjxetYq8PNR2ZmluOkYGZmOU4KY2tprQMYBcdaGY0UKzRWvI61CtynYGZmOa4pmJlZjpPCKEk6QdJzktZJ+nyR98+TtEnSquRxfi3iTGK5XtLLklaP8L4kfSv5szwt6ahqx5gXS7lYj5O0Ne9zvazaMebFMlvSQ5LWSloj6ZIiZeris00Zaz19tm2S/kPSU0m8XylSZrKk25LPdqWkedWPNHWsdfN9kFpE+JHyATQDvwDmA5OAp4BDCsqcB1xd61iTWI4FjgJWj/D+ScB9ZLaqXgSsrONYjwPuqfVnmsRyAHBU8rwd+M8ifw/q4rNNGWs9fbYC9k6etwIrgUUFZf4MWJI8Pwu4rY5jrZvvg7QP1xRGZyGwLiLWR8RO4Fbg1BrHNKKIWA78pkSRU4HvRcYKYF9JB1QnuqFSxFo3ImJjRDyRPO8F1gIzC4rVxWebMta6kXxeryaHrcmjsOPzVOCG5PkdwO9LUpVCzEkZa8NxUhidmcALeccbKP4P7ENJk8EdkmZXJ7TdkvbPUy+OSarq90k6tNbBACRNF0eS+ZWYr+4+2xKxQh19tpKaJa0CXgYeiIgRP9uIGAC2AtOqG2VGilihcb4PACeF0Sr2a6Twl8HdwLyIeAfwIK//oqlHaf489eIJMtP0DweuAn5Q43iQtDdwJ3BpRGwrfLvIKTX7bMvEWlefbUTsiogjgFnAQkmHFRQWwYJuAAAFW0lEQVSpm882RayN9H0AOCmM1gYgP9PPAl7KLxARmyNiR3L4HWBBlWLbHWX/PPUiIrZlq+oRcS/QKml6reKR1ErmS/amiLirSJG6+WzLxVpvn21WRGwBHgZOKHgr99lKagHeSI2bHkeKtcG+DwAnhdF6HDhI0pslTSLTybUsv0BBu/EpZNpw69Uy4KPJSJlFwNaI2FjroIqRtH+23VjSQjJ/dzfXKBYB1wFrI+KKEYrVxWebJtY6+2xnSNo3eT4FOB74WUGxZcCfJM9PB34cSa9uNaWJtcG+DwBoqXUAjSQiBiRdDPyQzEik6yNijaTFQGdELAM+LekUYIDMr5fzahWvpFvIjCyZLmkD8GUynWFExBLgXjKjZNYB24GP1SbSVLGeDnxK0gDwGnBWLb4IEu8CzgWeSdqTAb4IzIG6+2zTxFpPn+0BwA2Smskkp3+JiHsK/o1dB3xf0joy/8bOquNY6+b7IC3PaDYzsxw3H5mZWY6TgpmZ5TgpmJlZjpOCmZnlOCmYmVmOk4LVDUnT8laT/JWkF/OOJ1Xwvo9IOiJ5/kNJ7SXK/oWktkrFUkoS53OSTso7PmIPr3lw8vluGZsordF5noLVjYjYDGS/nC8HXo2Ib+SXSSZZKSIGKxTDH5Qp8hfA9UBfJe6fwpkRsap8sXQi4jlJHcArY3VNa2yuKVjdk3SgpNWSlpBZp2d2/i9bSWdJujZ5/iZJd0nqTNa6X1TkentJuj1ZpOxWoC3vvQ2S9pXUniwO91Ry79MlfQbYD/ippAeT8kuTe61R3j4EyXUul/Rkcp+3Jq+3S7pB0jPJ6x9MXj9R0mOSnlBmr4A3jOLzaZZ0Y3K/FklbJH0zudYDkqYl5d4q6cfJn+kJ1WgfAqtvTgrWKA4BrouII4EXS5T7FvD1iOgA/hi4tkiZi4GeZJGyr5FZObTQScAvI+LwiDiMzAqY3ySzGubvRcTxSbnPJ/c6HHi/pEPyrvHrJN5rydQwAC4HNkXE7yTn/ETSfsDngd+PiKOAp4Fhm+GMoAW4GXgmIi5PXnsjsCK51mPAXyev3wJ8M1n47neTP4vZEG4+skbxi4h4PEW544GD9fry+lMlTYmI1/LKHAt8HSAinpS0psh1ngb+TtLfAXdHxL+PcL8PS/o4mX9Lv00meT2bvJddfK6LTJLJxvfB5N4B9CS1hUOAR5O4JwGPpPizQmbJh5sj4mt5rw0AtyfPbwRuljQVmB4Rdyf3rlXzl9U5JwVrFP+d93yQocsn53f8CliYbIJUSsn1XSJibdLWfhLw95LuiYiv5peRdBCZX/QLI2KLpBsLYsmujrmL1/+tqci9BdwfEeeWibmYfyezycw/5K3GWXj9GOF1s2HcfGQNJ+lk7pF0kKQm4I/y3n4QuCh7MMLonOXAOcn7hwPDNpWRNJNMR/f3gSvIbBUK0EtmW0uAfZLjbclqmOU6qQH+jUzzVXYf56nAo8B7JM1PXn9DknDSWErmz3yrMstIQ2YhwdOS52cDj0RED/CKpJOTe7RJ2ivlPWwCcVKwRvVXwP3Aj8isr591EfCupBP3WeATRc69Gpgm6WngM0BnkTKHA48rs7Lo54BsLWEp8GDS0fwEmaai1WTWyh+piSnfV4A3SVoNrCLTP/Fr4OPAbZKeIpMk3priWgBExNeTOP6ZzL/prcBRkp4A3g38bVL0HOCzyZ/7EWBG2nvYxOFVUs0ahKRHgItLDUlNaguvRMS+o7juqM+x8cs1BbPGsRm4MTt5bSxIOphMTenXY3VNa2yuKZiZWY5rCmZmluOkYGZmOU4KZmaW46RgZmY5TgpmZpbjpGBmZjn/H33dd3cTszNHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(df.true_distance_kpc,\n",
    "             np.mean(simple_trace[\"D\"], axis=0),\n",
    "             yerr=np.std(simple_trace[\"D\"], axis=0), fmt=\".\");\n",
    "plt.xlabel(\"True distance [kpc]\")\n",
    "plt.ylabel(\"Measured distance [kpc]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [D, logL]\n",
      "Sampling 4 chains: 100%|██████████| 40000/40000 [00:16<00:00, 2408.65draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as hierarchical_model:\n",
    "    \n",
    "    logL = pm.Normal(\"logL\", mu=np.log(500), sd=10)\n",
    "    L = pm.Deterministic(\"L\", pm.math.exp(logL))\n",
    "        \n",
    "    # Prior over distance as a pymc3 potential object\n",
    "    D = pm.Uniform(\"D\", lower=0, upper=1000, shape=len(df),\n",
    "                   testval=np.array(1. / df.parallax_mas))\n",
    "    r = 1e3 * D\n",
    "    pm.Potential(\n",
    "        \"lnprior\", -3*logL + 2*pm.math.log(r) - r/L)\n",
    "    \n",
    "    # The model\n",
    "    mu = 1. / D\n",
    "\n",
    "    # Likelihood\n",
    "    pm.Normal(\n",
    "        \"parallax_obs\", mu=mu, sigma=np.array(df.parallax_error),\n",
    "        observed=np.array(df.parallax_mas))\n",
    "    \n",
    "    hierarchical_trace = pm.sample(tune=5000, draws=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFBRJREFUeJzt3X+QXeV93/H3J1Co7NYWPwTDrOSKxBqnjqe1mS2QMuOh0BB+uJY7Y6a4bq1SzaiZQuLWbY1I/iCT1B15+gObSaKOArIhdZEpTQZNoppokBlPM4UgYoKNSYqKCVqhoPWAaBuIGexv/7jP1teru7vSvbt7tXver5k795znPPfe55mzez/3POdXqgpJUvf8yLgbIEkaDwNAkjrKAJCkjjIAJKmjDABJ6igDQJI6ygCQpI4yACSpowwASeqoM8fdgPmcf/75tXHjxnE3Q5JWlCeffPI7VbVuoXqndQBs3LiRgwcPjrsZkrSiJPmTk6nnEJAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FELBkCS3UmOJfnmgGX/MkklOb/NJ8ldSQ4leTrJJX11tyR5rj22LG43JEmn6mTOBP4i8CvAff2FSTYAPwW82Fd8HbCpPS4DdgKXJTkXuAOYBAp4Msneqnp11A6o267YcYAjx98YuGxi7Rp+b/tVy9wiaeVYMACq6mtJNg5YdCfwaeChvrLNwH1VVcBjSdYmuQi4EthfVa8AJNkPXAvcP1Lr1XlHjr/BCztuGLhs4/bfWebWSCvLUPsAknwYOFJVfzhr0QRwuG9+qpXNVS5JGpNTvhhckrcBvwBcM2jxgLKap3zQ+28DtgG8613vOtXmSZJO0jBXA/0x4GLgD5MArAf+IMml9H7Zb+irux54qZVfOav80UFvXlW7gF0Ak5OTA0NC3TPXWP/E2jVjaI20OpxyAFTVN4ALZuaTvABMVtV3kuwFbk2yh95O4Neq6miSh4F/k+Sc9rJrgNtHbr06Y76xfknDOZnDQO8H/gfwniRTSbbOU30f8DxwCPh14J8CtJ2/vww80R6/NLNDWJI0HidzFNDHFli+sW+6gFvmqLcb2H2K7ZOGNrF2zcAjgTw8VOo5re8IJo1iri95Dw+VerwUhCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkd5WGg6hzPD5B6DACdVpbjmj+eHyD1GAA6rXjNH2n5uA9AkjrKAJCkjjIAJKmjDABJ6ih3AmssvMOXNH4GgMbCo32k8XMISJI6ygCQpI4yACSpo9wHIDVeI0hds2AAJNkNfAg4VlXva2X/Fvg7wJvA/wJurqrjbdntwFbge8DPVdXDrfxa4PPAGcDdVbVj8bsjDc9rBKlrTmYI6IvAtbPK9gPvq6q/BvxP4HaAJO8FbgJ+or3m15KckeQM4FeB64D3Ah9rdSVJY7JgAFTV14BXZpX9blW91WYfA9a36c3Anqr6blV9GzgEXNoeh6rq+ap6E9jT6kqSxmQxdgL/Y+C/tekJ4HDfsqlWNle5JGlMRgqAJL8AvAV8aaZoQLWap3zQe25LcjDJwenp6VGaJ0max9ABkGQLvZ3DH6+qmS/zKWBDX7X1wEvzlJ+gqnZV1WRVTa5bt27Y5kmSFjBUALQjem4DPlxVr/ct2gvclOTsJBcDm4DfB54ANiW5OMlZ9HYU7x2t6ZKkUZzMYaD3A1cC5yeZAu6gd9TP2cD+JACPVdXPVNUzSR4AvkVvaOiWqvpee59bgYfpHQa6u6qeWYL+SJJO0oIBUFUfG1B8zzz1PwN8ZkD5PmDfKbVOkrRkvBSEJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkd5T2BpQV4r2CtVgaAtADvFazVyiEgSeootwC0pK7YcYAjx984oXxi7ZoxtEZSPwNAS+rI8Td4YccN426GpAEcApKkjjIAJKmjDABJ6igDQJI6ygCQpI4yACSpoxYMgCS7kxxL8s2+snOT7E/yXHs+p5UnyV1JDiV5Osklfa/Z0uo/l2TL0nRHknSyTmYL4IvAtbPKtgOPVNUm4JE2D3AdsKk9tgE7oRcYwB3AZcClwB0zoSFJGo8FA6Cqvga8Mqt4M3Bvm74X+Ehf+X3V8xiwNslFwE8D+6vqlap6FdjPiaEiSVpGw+4DuLCqjgK05wta+QRwuK/eVCubq/wESbYlOZjk4PT09JDNkyQtZLF3AmdAWc1TfmJh1a6qmqyqyXXr1i1q4yRJPzBsALzchnZoz8da+RSwoa/eeuClecolSWMybADsBWaO5NkCPNRX/ol2NNDlwGttiOhh4Jok57Sdv9e0MknSmCx4NdAk9wNXAucnmaJ3NM8O4IEkW4EXgRtb9X3A9cAh4HXgZoCqeiXJLwNPtHq/VFWzdyxLkpbRggFQVR+bY9HVA+oWcMsc77Mb2H1KrZMkLRnPBJakjjIAJKmjDABJ6igDQJI6ynsCa1F483dp5TEAtCi8+bu08jgEJEkd5RaANKSJtWvYuP13Bpb/3varxtAi6dQYANKQ5vqSHxQK0unIISBJ6igDQJI6ygCQpI4yACSpowwASeooA0CSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjhopAJL88yTPJPlmkvuT/MUkFyd5PMlzSb6c5KxW9+w2f6gt37gYHZAkDWfoAEgyAfwcMFlV7wPOAG4CPgvcWVWbgFeBre0lW4FXq+rdwJ2tniRpTEYdAjoTWJPkTOBtwFHgKuDBtvxe4CNtenObpy2/OklG/HxJ0pCGDoCqOgL8O+BFel/8rwFPAser6q1WbQqYaNMTwOH22rda/fNmv2+SbUkOJjk4PT09bPMkSQsY+n4ASc6h96v+YuA48F+A6wZUrZmXzLPsBwVVu4BdAJOTkycs13h5719p9RjlhjB/G/h2VU0DJPlN4G8Ca5Oc2X7lrwdeavWngA3AVBsyeifwygifrzHw3r/S6jHKPoAXgcuTvK2N5V8NfAv4KvDRVmcL8FCb3tvmacsPVJW/8CVpTEbZB/A4vZ25fwB8o73XLuA24FNJDtEb47+nveQe4LxW/ilg+wjtliSNaKR7AlfVHcAds4qfBy4dUPfPgRtH+TxJ0uLxpvDSIptYu2bgjeEn1q6Z80by0jgYANIim+tLflAoSOPktYAkqaMMAEnqKANAkjrKAJCkjjIAJKmjDABJ6igDQJI6ygCQpI4yACSpowwASeooA0CSOsoAkKSOMgAkqaMMAEnqKANAkjrK+wFooCt2HODI8TdOKJ9Yu2YMrZG0FAwADXTk+Bu8sOOGcTdD0hJyCEiSOmqkAEiyNsmDSf4oybNJfjLJuUn2J3muPZ/T6ibJXUkOJXk6ySWL0wVJ0jBGHQL6PPCVqvpokrOAtwE/DzxSVTuSbAe2A7cB1wGb2uMyYGd7ljrBm8XrdDN0ACR5B/BB4B8BVNWbwJtJNgNXtmr3Ao/SC4DNwH1VVcBjbevhoqo6OnTrpRXEm8XrdDPKENCPAtPAF5J8PcndSd4OXDjzpd6eL2j1J4DDfa+famU/JMm2JAeTHJyenh6heZKk+YwSAGcClwA7q+oDwJ/RG+6ZSwaU1QkFVbuqarKqJtetWzdC8yRJ8xklAKaAqap6vM0/SC8QXk5yEUB7PtZXf0Pf69cDL43w+ZKkEQwdAFX1p8DhJO9pRVcD3wL2Alta2RbgoTa9F/hEOxrocuA1x/8laXxGPQroZ4EvtSOAngduphcqDyTZCrwI3Njq7gOuBw4Br7e6kqQxGSkAquopYHLAoqsH1C3gllE+T5K0eDwTWJI6ygCQpI4yACSpo7waqDRmXiJC42IASGPmJSI0Lg4BSVJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkd5YlgHXfFjgMcOf7GCeUTa9eMoTWSlpMB0HFHjr/BCztuGHczJI2BQ0CS1FEGgCR1lAEgSR1lAEhSRxkAktRRIwdAkjOSfD3Jb7f5i5M8nuS5JF9OclYrP7vNH2rLN4762ZKk4S3GYaCfBJ4F3tHmPwvcWVV7kvxHYCuwsz2/WlXvTnJTq/f3FuHzpVXJO4VpqY0UAEnWAzcAnwE+lSTAVcDfb1XuBX6RXgBsbtMADwK/kiRVVaO0QVqtvFOYltqoQ0CfAz4NfL/Nnwccr6q32vwUMNGmJ4DDAG35a62+JGkMhg6AJB8CjlXVk/3FA6rWSSzrf99tSQ4mOTg9PT1s8yRJCxhlC+AK4MNJXgD20Bv6+RywNsnM0NJ64KU2PQVsAGjL3wm8MvtNq2pXVU1W1eS6detGaJ4kaT5DB0BV3V5V66tqI3ATcKCqPg58Ffhoq7YFeKhN723ztOUHHP+XpPFZivMAbqO3Q/gQvTH+e1r5PcB5rfxTwPYl+GxJ0klalKuBVtWjwKNt+nng0gF1/hy4cTE+T5I0Os8ElqSO8n4AHeGNXyTNZgB0hDd+kTSbQ0CS1FFuAUgrjNcI0mIxAKQVxmsEabE4BCRJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSR3kYqLRKeH6ATpUBIK0Snh+gU+UQkCR1lFsAq4xX/ZR0sgyAVcarfko6WQ4BSVJHGQCS1FEGgCR11NABkGRDkq8meTbJM0k+2crPTbI/yXPt+ZxWniR3JTmU5OkklyxWJyRJp26ULYC3gH9RVX8VuBy4Jcl7ge3AI1W1CXikzQNcB2xqj23AzhE+W5I0oqGPAqqqo8DRNv1/kjwLTACbgStbtXuBR4HbWvl9VVXAY0nWJrmovY+kJTLXGcIzyzxLuLsW5TDQJBuBDwCPAxfOfKlX1dEkF7RqE8DhvpdNtTIDQFpC833Be5Zwt428EzjJXwL+K/DPqup/z1d1QFkNeL9tSQ4mOTg9PT1q8yRJcxhpCyDJX6D35f+lqvrNVvzyzNBOkouAY618CtjQ9/L1wEuz37OqdgG7ACYnJ08ICM19ti94xq+kkzd0ACQJcA/wbFX9h75Fe4EtwI72/FBf+a1J9gCXAa85/j8cz/aVtBhG2QK4AviHwDeSPNXKfp7eF/8DSbYCLwI3tmX7gOuBQ8DrwM0jfLYkaUSjHAX03xk8rg9w9YD6Bdwy7OdJkhaXZwJLUkcZAJLUUV4OWuowbyPZbQaA1GHeRrLbHAKSpI5yC+A05u0dJS0lA+A05glfkpaSQ0CS1FEGgCR1lENAkk7g4aHdYABIOoGHh3aDQ0CS1FEGgCR1lAEgSR3lPoDTgCd8aaVw5/DqYgCcBjzhSyuFO4dXF4eAJKmj3AKQNDKHhlYmA2AZOdav1WquL/krdhwwGE5jBsAycqxfXWMwnN6WPQCSXAt8HjgDuLuqdix3G5aav/Sl+bkz+fSwrAGQ5AzgV4GfAqaAJ5LsrapvLWc7lpq/9KXhuC9heS33FsClwKGqeh4gyR5gM3BaB8B8v+j9o5QWj0NGy2u5A2ACONw3PwVctsxtmNN8X/SDftHP90cpafGcajAsptUcMssdABlQVj9UIdkGbGuz/zfJHy95qxbwJ0BuP+X65wPfWaImne7sezetyr6fwv//6dT/v3IylZY7AKaADX3z64GX+itU1S5g13I2aikkOVhVk+NuxzjYd/veRSux/8t9JvATwKYkFyc5C7gJ2LvMbZAkscxbAFX1VpJbgYfpHQa6u6qeWc42SJJ6lv08gKraB+xb7s8dgxU/jDUC+95NXe47rMD+p6oWriVJWnW8GqgkdZQBMIIkZyT5epLfbvMXJ3k8yXNJvtx2dJPk7DZ/qC3fOM52L4YBff9ikm8neao93t/Kk+Su1venk1wy3paPLskLSb7R+nmwlZ2bZH9b9/uTnNPKV1X/5+j7LyY50rfur++rf3vr+x8n+enxtXx0SdYmeTDJHyV5NslPrvT1bgCM5pPAs33znwXurKpNwKvA1la+FXi1qt4N3NnqrXSz+w7wr6rq/e3xVCu7DtjUHtuAncvYxqX0t1o/Zw772w480tb9I20eVmf/Z/cden/3M+t+H0CS99I70u8ngGuBX2uXg1mpPg98pap+HPjr9P7+V/R6NwCGlGQ9cANwd5sPcBXwYKtyL/CRNr25zdOWX93qr0iz+76AzcB91fMYsDbJRUvawPHoX8ez130X+j/IZmBPVX23qr4NHKJ3OZgVJ8k7gA8C9wBU1ZtVdZwVvt4NgOF9Dvg08P02fx5wvKreavNT9C59AX2XwGjLX2v1V6rZfZ/xmba5e2eSs1vZoMt/TLCyFfC7SZ5sZ64DXFhVRwHa8wWtfLX1f1DfAW5t6373zDAIq6vvPwpMA19oQ593J3k7K3y9GwBDSPIh4FhVPdlfPKBqncSyFWWOvgPcDvw48DeAc4HbZl4y4G1WZN/7XFFVl9DbzL8lyQfnqbva+j+o7zuBHwPeDxwF/n2ru5r6fiZwCbCzqj4A/Bk/GO4ZZEX03QAYzhXAh5O8AOyhN/TzOXqbeTPnVvRf5uL/XwKjLX8n8MpyNngRndD3JP+pqo62zd3vAl/gB5v6C17+Y6Wpqpfa8zHgt+j19eWZTfz2fKxVX1X9H9T3qnq5qr5XVd8Hfp3Vue6ngKmqerzNP0gvEFb0ejcAhlBVt1fV+qraSG8n14Gq+jjwVeCjrdoW4KE2vbfN05YfqBV6AsYcff8Hff8EoTcO+s32kr3AJ9pREZcDr81sMq9ESd6e5C/PTAPX0Otr/zqeve5XRf/n6vusse2/yw+v+5vaUXAX09sh+vvL2ebFUlV/ChxO8p5WdDW9y9iv6PXuLSEX123AniT/Gvg6bYdRe/6NJIfo/fK/aUztW0pfSrKO3qbvU8DPtPJ9wPX0dgC+Dtw8nuYtmguB32r78M8E/nNVfSXJE8ADSbYCLwI3tvqrqf9z9f030jvst4AXgH8CUFXPJHmA3hflW8AtVfW9sbR8cfwsvb/zs4Dn6a3LH2EFr3fPBJakjnIISJI6ygCQpI4yACSpowwASeooA0CSOsoAkKSOMgAkqaMMAEnqqP8HCWh8b7Pv/pYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(hierarchical_trace[\"L\"], 50, histtype=\"step\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bonus exercises\n",
    "\n",
    "Choose one of the following:\n",
    "\n",
    "*1) Redo this analysis with emcee.*\n",
    "\n",
    "*2) Redo this analysis, in PyMC3, where you don't know, or don't trust, the uncertainties on the parallaxes.*\n",
    "\n",
    "*3) Reproduce the file \"parallaxes.csv\". How did I sample from the prior distribution?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lnprob(r_kpc, parallax_mas, parallax_mas_error):\n",
    "    lnlike = -.5*np.sum((1./r_kpc - parallax_mas)**2/parallax_mas_error**2)\n",
    "    \n",
    "    if r_kpc > 0:\n",
    "        lnprior = log_bj_prior(r_kpc, 500)\n",
    "    else:\n",
    "        lnprior = -np.inf\n",
    "    return lnlike + lnprior\n",
    "\n",
    "print(lnprob(1, 1, .1))\n",
    "print(lnprob(1, 5, .1))\n",
    "print(lnprob(5, 1./5, .001))\n",
    "print(lnprob(5, 10, .001))\n",
    "\n",
    "import emcee\n",
    "import corner\n",
    "\n",
    "ndim, nwalkers = 1, 100\n",
    "p0 = [np.random.rand(ndim) for i in range(nwalkers)]\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[df.parallax_mas[0], df.parallax_error[0]])\n",
    "sampler.run_mcmc(p0, 1000, progress=True);\n",
    "\n",
    "samples = sampler.chain[:, 50:, :].reshape((-1, ndim))\n",
    "print(\"true = \", df.true_distance_kpc[0])\n",
    "print(\"inferred = \", np.median(samples, axis=0))"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}