Gaia_exercise.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": 87,
   "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": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv(\"parallaxes.csv\")\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": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '1/parallax')"
      ]
     },
     "execution_count": 99,
     "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",
    "\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": [
    "*1. Draw the non-hierarchical PGM for this problem.*\n",
    "\n",
    "*2. Draw the hierarchical PGM for this problem*\n",
    "\n",
    "*3. Using PyMC3, infer the distance to each star with the non-hierarchical model. To define the prior with PyMC3, you'll need to use a ['Potential'](https://pymc-devs.github.io/pymc/modelbuilding.html)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "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:14<00:00, 2670.02draws/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": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"392pt\" height=\"171pt\"\n",
       " viewBox=\"0.00 0.00 392.00 171.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 167)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-167 388,-167 388,4 -4,4\"/>\n",
       "<g id=\"clust1\" class=\"cluster\">\n",
       "<title>cluster100</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M20,-8C20,-8 364,-8 364,-8 370,-8 376,-14 376,-20 376,-20 376,-143 376,-143 376,-149 370,-155 364,-155 364,-155 20,-155 20,-155 14,-155 8,-149 8,-143 8,-143 8,-20 8,-20 8,-14 14,-8 20,-8\"/>\n",
       "<text text-anchor=\"middle\" x=\"357\" y=\"-15.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">100</text>\n",
       "</g>\n",
       "<!-- parallax_obs -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>parallax_obs</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"110\" cy=\"-57\" rx=\"94.4839\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"110\" y=\"-53.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">parallax_obs ~ Normal</text>\n",
       "</g>\n",
       "<!-- D -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>D</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"202\" cy=\"-129\" rx=\"57.6901\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"202\" y=\"-125.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">D ~ Uniform</text>\n",
       "</g>\n",
       "<!-- D&#45;&gt;parallax_obs -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>D&#45;&gt;parallax_obs</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M180.1973,-111.937C168.4285,-102.7266 153.7056,-91.2044 140.7868,-81.094\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"142.6391,-78.0993 132.607,-74.6924 138.325,-83.6118 142.6391,-78.0993\"/>\n",
       "</g>\n",
       "<!-- lnprior -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>lnprior</title>\n",
       "<polygon fill=\"none\" stroke=\"#000000\" points=\"367.5,-75 222.5,-75 222.5,-39 367.5,-39 367.5,-75\"/>\n",
       "<text text-anchor=\"middle\" x=\"295\" y=\"-53.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">lnprior ~ Deterministic</text>\n",
       "</g>\n",
       "<!-- D&#45;&gt;lnprior -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>D&#45;&gt;lnprior</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M224.0397,-111.937C235.7602,-102.8631 250.379,-91.5453 263.297,-81.5443\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"265.7281,-84.0885 271.4927,-75.1992 261.4429,-78.5534 265.7281,-84.0885\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1c2af93e48>"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.model_to_graphviz(simple_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the resulting measured distances against the true distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Measured distance [kpc]')"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucXHV9//HXezcbFjRITILQhCRGIgpYIFlxqYposQI/BeVSbtJiRbSFCtaf1sujiNT6q9afFsEK4WJRrnL52cAPsKIgUkjKLgYIRGtMXQhECWFDNiabZLOf/nHOTGZnZ2fOJju33ffz8ZgHc+Z855xPhmQ+870rIjAzMwNoqXcAZmbWOJwUzMwsz0nBzMzynBTMzCzPScHMzPKcFMzMLM9JwczM8pwUzMwsz0nBzMzyJtU7gNGaPn16zJ07t95hmJk1le7u7hcjYkalck2XFObOnUtXV1e9wzAzayqSerKUc/ORmZnlOSmYmVmek4KZmeU5KZiZWZ6TgpmZ5TkpmJlZnpOCmZnlOSmYmVmek4KZWZWceuUjnHrlI/UOY1ScFMzMLM9JwczM8pwUzMwsz0nBzMzynBTMzCzPScHMzPKcFMzMLM9JwcysCdRqzoOTgpmZ5TkpmJlZXtPt0Wxm1uiabWmLQq4pmJlZnpOCmZnlVS0pSGqX9J+SHpf0lKQvlihztqS1kpalj3OqFY+ZmVVWzT6FLcC7ImKjpDbgIUn3RMSSonK3RMT5VYzDzKwu+vq3saF/gO6eXhbOmVrvcDKpWk0hEhvTw7b0EdW6n5lZo+jr38Z/v7iRFWv6WN27mTOvXkJ3T2+9w8qkqn0KklolLQNeAH4UEUtLFDtJ0hOSbpO03wjXOVdSl6SutWvXVjNkM7Nd0t3Tyy9+28cLfVvzv4K3DQyyZNW6usaVVVWTQkRsj4hDgVnA4ZIOLipyJzA3Iv4QuA+4boTrLIqIjojomDFjRjVDNjPbJUtWrWOwqE2kbVILnfOm1SegUarJ6KOIWA88ABxT9Pq6iNiSHl4FLKxFPGZm1dI5bxotSp4LmDFlMjec09k0fQpV62iWNAPYFhHrJe0OHA18pajMvhGxJj08HlhRrXjMzGph4ZypvGGfKWzoH2DP9klMaW9rmoQA1R19tC9wnaRWkhrJ9yPiLkmXAF0RsRj4uKTjgQHgJeDsKsZjZlYTU9rbmNLeVu8wdkrVkkJEPAEcVuL1iwqefxb4bLViMDMbL2o1vNUzms3MGlxuRFMthrc6KZhZQ6rV/gHNoHBEU7WHtzopmJk1uMIRTdUe3uqls83MxtgtHz0CGLsltAtHNF162mFV7VNwUjAzawK5EU3VHt7q5iMza2ruexhbTgpmZpbnpGBmZnlOCmZmluekYGZmeR59ZGZWJbmhqc3ENQUzM8tzUjCzptbXv43n1m9umu0uG52Tgpk1rVouFLezmm0ehZOCmTWtWi4UN1E4KZhZ06rlQnETxYijjyT9TYb3/z4irhzDeMzMgGybytRyobiJolxN4VPAK4EpZR6frHaAZjbxjKavYEp7GzP32n3cJ4RbPnpETYa4lpun8L2IuKTcmyW9osy5duBBYLf0PrdFxBeKyuwGfBdYCKwDTo2I32QL3czGq1J9BeP9S79RjFhTiIhPV3pzhTJbgHdFxCHAocAxkjqLynwY6I2I/YFvAF+pHLKZjbVGGyHjvoL6qdjRLOnLkvYqOJ4q6UuV3heJjelhW/qIomInANelz28D/liSMkVuZuNCqYSU6yuYNXV3bjins6lrCc02jyLL6KNjI2J97iAieoHjslxcUqukZcALwI8iYmlRkZnAs+l1B4CXAf8kMLNx0VfQDPMoimVJCq1p2z8AknYn6SeoKCK2R8ShwCzgcEkHFxUpVSsork0g6VxJXZK61q5dm+XWZmZ114zzKLIsiHc98GNJ3yH5wv4LdjT5ZBIR6yU9ABwDLC84tRrYD1gtaRLwKuClEu9fBCwC6OjoGJY0zKx5ZRl6Wk4jLzqX6xsZjObpG6lYU4iIrwJfAt4IHAj8ffpaWZJm5Poi0trF0cAviootBv48fX4y8JOI8Je+2QTRjM0ro9GMfSNZl85eRtJkNAg8nvE9+wLXSWolST7fj4i7JF0CdEXEYuAa4HuSVpLUEE4bVfRm1tQmwtDTKe1tTGlva5o/V8WkIOkc4CLgJyR9AJdJuiQiri33voh4AjisxOsXFTzvB04ZbdBmNj40Y/PKeJelpvAp4LCIWAcgaRrwMFA2KZiZVeJlKhpPlqSwGugrOO4jHUZqZs0jNxeg0Tpmm615ZbzLkhSeA5ZK+jeS0UcnAP+ZWzAvIr5exfjMrAZ2dQTQrmi0JDXRZZmn8GvgB+yYP/BvwBp2LIpnZk1svI8AstHJUlO4rniROklvjohHqxOSmdXSRBgBVE/NVhPKUlO4XdLM3IGkI3Ens9m4Ue3F57p7evnW/StdA2kSWWoKHwN+IOl9wALgy2Rc+8jMGt9YjwAq7NDu7unllCseZjCgva2laSZwTWQVk0JEPCrp48C/A/3AuyPCCxCZjSM7OwKo1Iimwk7rXWmaarZml/Gi3HacdzJ0cbo9SFYxvUYSEXF8tYMzs+aS67QeDDjz6iVc9N6DaG9rYdvAoCenNYlyNYWv1SwKMxsXimsGvZu2csM5nSxZtY7OedPcdNQERkwKEfHTWgZiZs2v1LIVC+dMdTJoIiOOPpJ0V6U3ZyljZrum0bbKLKcZVwW1oco1H71N0uIy50WylLaZNYFazVq++4Ijq3Ztq75ySeGEDO/fOlaBmFn1FHcAj9Wv+Houj2HV4T4FswmgGrOWq5VorL6yzGg2syZXjVnLzbj/sFWWdec1M2ti1di3wBvkjE+ZkkK6x/LsiPhlleMxsyopN2t5Z2YPe4Oc8ali81G65tEy4N70+NAKo5LMrE4qDV99es2GMR3eOqW9jZl77e6EMI5k6VO4GDgcWA8QEcuAuZXeJGk/SfdLWiHpKUkXlChzlKSXJS1LHxeVupbZRNbXv43n1m/2KqNWE1mSwkBEvLwT1x4APhkRbwQ6gfMklZrX8LOIODR9XLIT9zEbt7wBjtValqSwXNIZQKuk+ZIuAx6u9KaIWBMRj6XP+4AVwMzy7zKzQh7hY7WWJSn8NXAQsAW4kWSl1AtHcxNJc4HDgKUlTh8h6XFJ90g6aDTXNRvvRjuUNGtTUzMtnWG1lWU/hU3A59PHqEl6JXA7cGFEbCg6/RgwJyI2SjqOZC/o+SWucS5wLsDs2bN3JgyzpjSaET71mEzmPQ/Gnyyjj34kaa+C46mSfpjl4pLaSBLCDRFxR/H5iNgQERvT53cDbZKmlyi3KCI6IqJjxowZWW5tNm5kHeEzmqYmd17bSLI0H02PiPW5g4joBfau9CZJAq4BVkTE10cos09aDkmHp/G40dRsFHJNQVmamrYPBv/94kZWrHHntZWWZfLaoKTZEfEMgKQ5DN2RbSRvBc4CnpS0LH3tc8BsgIi4AjgZ+EtJA8Bm4LSIyHJtMytSqampr38bm7ZuZ9PW7fnXxmodJBs/siSFzwMPScotkHckaft+ORHxEMny2uXKXA5cniEGM8ug3KzlDf0Dw17z8hRWrGLzUUTcCywAbgG+DyyMiEx9CmbWOC497bB885KAGVMme2VTGybrgni7AS+l5Q+UREQ8WL2wzGysFTYv7dk+acQahU1sFZOCpK8ApwJPAYPpywE4KZjVQKVhn4Ub3VSSa17yUFIbSZaawvuBAyJiS7WDMbNschPPPn3MG7zRjY2pLENSVwFt1Q7EzEbPy2DYWMtSU9gELJP0Y5KlLgCIiI9XLSozy8Qb3dhYy5IUFqcPM2sw3ujGxlqWtY+uq0UgZpZdYedyubkJZqOVZfTRfOD/AAcC7bnXI2JeFeMysxEUL3w3b/ormNKerdvPo46skiwdzd8Bvk2yac47ge8C36tmUGY2suLO5VIzlc12VpaksHtE/BhQRPRExMXAu6oblpmNpHjhuz3bs85BNassy9+mfkktwK8knQ88R4ZVUs0mqtwcgmo11exs53K147LxIUtN4UJgD+DjwELgg8CfVTMoMysv6x4LZqOVpaYwNyIeBTYCHwKQdAqlt9Y0m7C8vaWNB1lqCp/N+JqZmTW5EWsKko4FjgNmSvpmwak9SUYimVkJhXMIqtW8434Bq5ZyNYXngS6gH+gueCwG3lP90MyaS1//NlZ5q0trciPWFCLiceBxSTdGxDYASVOB/dJ9ms0sVTihLMdbXVozytKn8CNJe0p6NfA48B1JX6/0Jkn7Sbpf0gpJT0m6oEQZSfqmpJWSnpC0YCf+DGZ1VzihLMcL1FkzypIUXhURG4ATge9ExELg6AzvGwA+GRFvBDqB8yQdWFTmWGB++jiXZOa0WdMpnFAmYO8G3Oqyr38bz63f7CYtKytLUpgkaV/gT4G7sl44ItZExGPp8z5gBTCzqNgJwHcjsQTYK72XWd2ceuUjZYeXljqfm1A2a+ruvHHfKbx2+it3KiFUuvfOyjVvua/DKsmSFC4BfgisjIhHJc0DfjWam0iaCxzG8LkNM4FnC45XMzxxmDWF3ISyrIvT1ZI347GssiydfStwa8HxKuCkrDeQ9ErgduDCtBlqyOlStyxxjXNJmpeYPXt21lub1VQjDxP1ZjyWVbl5Cp+OiK9KuowSX9RZdl6T1EaSEG6IiDtKFFkN7FdwPItkKGzxvRYBiwA6OjqGxWJm5XkzHsuqXE1hRfrfrp25sCQB1wArImKk0UqLgfMl3Qy8BXg5ItbszP3MxkotJp/V497ejMeyKDdP4c70vzu789pbgbOAJyUtS1/7HDA7ve4VwN0ks6ZXkuwF/aGdvJfZmCjewCY3gqgWK4yOdG+zWirXfHQnJZqNciLi+HIXjoiHKN1nUFgmgPMqxGhWM6U6ZBfOmTrkF3y1fs2PdG+zWirXfPS19L8nAvsA16fHpwO/qWJMZnVTqkO28Bf86VctYdvAIMHY/5p3Z7A1gnLNRz8FkPT3EXFkwak7JT1Y9cjM6qBUh+y37l855Bd8rvo81r/m3RlsjSDLfgozJM1Lh6Ii6bXAjOqGZVY/uQ7Zv/vBk2zoH+Cvjtp/yC/4XGKoxq95dwZbvWVJCp8AHpC0Kj2eSzpnwGy86uvflm8yuuSup5j96j0YGAwuPe2wfLLwr3kbj7JMXrtX0nzgDelLv4iILdUNy6y2CkcX3fLRI/jW/St5es0vgaSZaGAw8ttf+te8jWdZagqkSeDxKsdi1jCKO31rVSuo5pDXRp5xbY0jy9pHZuNe8QqihQvceb6ATSROCjbhjbSCaG6BOycEm0jKTV4ru+FNbllss3rb1dnGo5005mYYG8/K9Sn83/S/7UAHSZ+CgD8kWQL7bdUNzaw2Rpo0Vs81kMzqZcTmo4h4Z0S8E+gBFkRER7rr2mEkaxWZNZ3jLn2Qt33lJ0M2mSnVf+BNaWyiytKn8IaIeDJ3EBHLgUOrF5LZ6GTdZrLcF31x/4E3pbGJKsuQ1BWSriZZ+yiAD7JjWW2zujn1ykeGTDKrtBZRub6D4n4Cr0NkE1WWmsKHgKeAC4ALgafxEtfWIDb0Dwz5or/g5p+PuMdx7oseKn/Re0iqTVRZZjT3S7oCuDsiflmDmMwy27N9EiKpwra2trBn+8h/pRfOmcqtH/sjlqxaR+e8aRW/6D1z2SaiiklB0vHAPwGTgddKOhS4pNJ+CmbV1te/jRc3btmx6UcM3/6jeLjqwjlT/SVvVkaWPoUvAIcDDwBExDJJc6sXklllhXsc5GwfDNZu3DJkGOmuDCv1fASbiLL0KQxExMtVj8RsFAo7jXNaW8SLfVvzo4tuXPqMh5WajVKWpLBc0hlAq6T5ki4DHq5yXGYj6u7p5fn1m/N7vU6e1MIZb5nNKR37DdkA557lazys1GyUsiSFvwYOArYANwIvk4xCKkvStZJekLR8hPNHSXpZ0rL0cdFoAreJqbunl1OueJgblj4DwIwpk7npI518+QNv4sQFs4aMLjr24H0zjzYys0TZPgVJrcAXI+JTwOdHee1/BS4HvlumzM8i4r2jvK41uXJrFVVax6iw2SiA3Sa15vsKSm1nef2S33hDHLNRKFtTiIjtwMKduXBEPAi8tDPvtYmr0uzkwrkGLWLYENTimcle6dRsdLI0H/1c0mJJZ0k6MfcYo/sfIelxSfdIOmiMrmlN4Ok1G4ZNMjvu0gd5ek35juGFc6Yy+9V7sGf7JL70/jdx9wVHlr1Pbic1M8smS1J4NbAOeBfwvvQxFk0+jwFzIuIQ4DLgByMVlHSupC5JXWvXrh2DW1s99fVvY/O27fz3ixuHfPFv6B/IPx+pY7i7p5dnXtrEhv4BLrnrKY8oMhtjWWY0V2VJi4jYUPD8bkn/Iml6RLxYouwiYBFAR0fH8BlK1jS6e3pZsaaPAF7o28rpVy3hpo8ky0hMyrULMXLHcKW9D1wrMNs1WWY0fwcY9kUcEX+xKzeWtA/wu4gISYeT1Fo8ZnCcW7Jq3ZC/TIU1gmde2gQkm3Zc9N6DSvYDeKE6s+rKMqP5roLn7cAHgOcrvUnSTcBRwHRJq0lmRrcBRMQVwMnAX0oaADYDp0WUWKfAxpXOedPyaxXBji/2whpAi6B309aS7y81wsjMxk6W5qPbC4/TL/v7Mrzv9ArnLycZsmoTyMI5U3njvlNYu3ELf3LgPpy4YFb+iz1rDcAL1ZlVT5aaQrH5wOyxDsTGl8L5BoXPj7v0QVau/T37z3gF//CBN+XLj6YG4H4Ds+rJ0qfQx9A+hd8Cf1u1iGzcKlzE7he/7Ru2SJ1rAGb1l6X5aEotArHx5cnnXmZgMOju6aWvfxsr1/6ej13fle83GAy4/bHVTgBmDSZLTeGtwLKI+L2kDwILgEsjoqfq0VlT6u7pZdPW7QCccsXD+USwtm8rrS2wfTA5vq17NScV9Cm4Wcis/rJMXvs2sEnSIcCngR7Kr2dkE9yVP/11/nnx8taDgzueb9/ulUvNGk3W/RQCOIGkhnAp4CYlG9HvNvSXPZ+bouZ5BmaNJ8vooz5JnwU+CByZrpzaVt2wrJmd+ubZPL76yZLndmtr4aL3HkTvpq2Z9kk2s9rKUlM4lWQvhQ9HxG+BmSR7NpsNc9ylD/IvD6xknz13Y7dJQ/96vX3+dG44p5Mz3jKb///E81xw88+9dpFZg6mYFCLitxHx9Yj4WXr8TES4T8GGyQ05Xd27mfWbt/GF9x2UX+a6va2FC49+PQvnTB1SzttkmjWWiklBUqekRyVtlLRV0nZJ3rPZhilerK5301besM8UZk3dnRvO6cw3FZVa1M7MGkOWPoXLgdOAW4EO4M9IZjXbBFG8G1p3Ty9LVq0b1ifQOW8a7W0tbBsYzHcin/fO/Yddz4vamTWuTMtcRMRKSa3pTmzfkfRwleOyBtLXv40N/QP5Zp7c3IP2tpYhNYCFc6ZywzmdJRNGoYVzpnLrx/6oYjkzq70sSWGTpMnAMklfBdYAr6huWNYoCpemOPPqJZy4YFbZ/QwWzpma6Us+azkzq60sSeEskr6H84FPAPsBJ1UzKKu/XBPRc+s3D0kCgmFNRGY2fmRZ+6hH0u7AvhHxxRrEZDVQ3E9w3KUPsqF/gL86an+WP/8yNy19hgAmT2rJ73/QNqmFExfM4sQFs9z0YzZOZVn76H3A14DJwGslHQpcEhHHVzs4q57ifoJcE9Hn/t/QSWfbtw8yY8pkJk9qzS9pnUsopTqRzay5ZWk+uhg4HHgAICKWSZpbtYis6sr1ExQSSe3g2x/sGFIjKEworimYjS9Z1z7yvIRxpHiegNixHlGhdx/4miGjiwBPPDMb57IkheWSzgBaJc2XdBngIalNLDdPAJKawEF/8KqS5aZP2W1YTcATz8zGtyxJ4a+Bg0jWP7oJ2ABcWOlNkq6V9IKk5SOcl6RvSlop6QlJC0YTuO283DyBT73nAG44p5OLFi+nROtRydpDcULx6COz8SXL6KNNwOfTx2j8K8ls6JHWSTqWZGb0fOAtJPs2vGWU97AKcqOKLj3tMCD5pX9r17MMDEb+tYHtw1PC5HSkUbHR7KVsZs1nxKQgaXG5N1YafRQRD1bokD4B+G66V8MSSXtJ2jci1pS7rmVX2KH8p1c+wvai3uQzr17CkfNnDHvf2+dPzy9eV8rdFxxZlXjNrP7K1RSOAJ4laTJaSunWhF0xM71+zur0tWFJQdK5wLkAs2fPHuMwxq/C9v/ihABJn0DxhjitLSqbEMxsfCvXp7AP8DngYOBS4N3AixHx04j46Rjcu1SSKdW0TUQsioiOiOiYMWP4L1sbrrunl6t+9usRz+eGm5765tn5/xGtLeLvTzjYCcFsAhuxppAufncvcK+k3YDTgQckXRIRl43BvVeTLJmRMwt4fgyuO+F19/Ry8rcfLp1hgdYWOO3NszlxwSwWzpnKAftM8QxlMwMqdDSnyeB/kSSEucA3gTvG6N6LgfMl3UzSwfyy+xN23Y1Ln2HRg78eMSGIJCH8wwfelH/Ni9OZWU65jubrSJqO7gG+GBElh5aWef9NwFHAdEmrgS+Q7u0cEVcAdwPHASuBTcCHdiJ+K3Dj0meGLVNRSCR7JJcaVWRmBqBk8E+JE9Ig8Pv0sLCQgIiIPascW0kdHR3R1dVVj1s3tOMufZBfvbCRbQXDS6fs1krflu1AMsT05IWzOCltMjKziUVSd0R0VCpXrk8hy8Q2q6Pc8tZT95icH3pa6IjXTeej73id+wvMLLNMO69Z4+nu6c3vgDapRcMSQmsLfPQdr3N/gZmNimsDTapwDsLgYOSXnoAdnclOBmY2Wk4KTebGpc9w1jVL6du8jfa2FloFk9taOPft85jUIlrkzmQz23luPmoihaOLfvarF/nYkfOYsntbvr/g3Qft4/4DM9slTgpN5JZHnxly/NSaDXzvwzvWEHT/gZntKieFBlc4wujx1UP3Ojr24H3rFJWZjVdOCg2seIRRoT858DWc8RYvDmhmY8tJoUF19/Tyz/f915ARRpNaRETQNqmFj77jdfUN0MzGJSeFBlRYQwBoUTIj+aL3HkTvpq3uSDazqnFSaECFcxBagLfuX37TGzOzseJ5Cg2oc960IXMQnBDMrFZcU6ij3Mii4uaghXOmcsM5nZ5zYGY156RQJ909vZx59RK2DgwyeVILN5zTOSwxOBmYWa05KdRAd08vdzy2moD80tVLVq2jf9sgkOyVvGTVOicBM6s7J4Uq6+7p5fSrkhoBwG1dz3LTuUfk+w22DQzSNqmFznnT6hypmZmTQlV19/RyyZ1P5RMCwLbtwZJV6zjvnfu738DMGo6TQpXcuPQZ/u7flrO9aKODtlblawXuNzCzRlPVIamSjpH0S0krJX2mxPmzJa2VtCx9nFPNeGqlu6e3ZEI4ZNaruOncI5wIzKxhVa2mIKkV+BbwbmA18KikxRHxdFHRWyLi/GrFUQ93PLZ6WEKY3Couet9BTghm1tCqWVM4HFgZEasiYitwM3BCFe/XMF7o2zLkeOZe7a4hmFlTqGZSmAk8W3C8On2t2EmSnpB0m6T9qhhP1XT39PKt+1fS3dMLwN5Tdhty/qgD9nZCMLOmUM2OZpV4rWh7ee4EboqILZI+BlwHvGvYhaRzgXMBZs9urOWiCxeva29LJqGduGAWt3avzg839daYZtYsqllTWA0U/vKfBTxfWCAi1kVErq3lKmBhqQtFxKKI6IiIjhkzZlQl2KyKawVLVq3LnyuchHbTRzr53+85gJs+0ulagpk1jWrWFB4F5kt6LfAccBpwRmEBSftGxJr08HhgRRXj2WWllqbonDeNyZOGT0LzcFMza0ZVSwoRMSDpfOCHQCtwbUQ8JekSoCsiFgMfl3Q8MAC8BJxdrXh2VW4iWvHSFJ6EZmbjiSKKm/kbW0dHR3R1ddXsft09vdz+2Gpu7XqWbdt3fFaTW+URRWbWNCR1R0RHpXKe0VxGrrloy7bBYT3kp3Ts54RgZuOON9kpY8mqdWwdGJ4QJrfKI4rMbFxyTaGMwk7k1hZx1AF7M33Kbvnlr83MxhsnhTK8A5qZTTROChV4aKmZTSQTuk+heCKamdlENyFrCrntMW/tepaBwSi5R7KZ2UQ04ZJCqWGm3iPZzCwxYZJCd08vS1at4/n1m/OzkiFZtc97JJuZJSZEUihcs2hSi5g8qYXt2wdpbW3h5IWzPMTUzCw1IZJCbhLaYMD2weDUw/dj5l67e5ipmVmRCZEUilcydc3AzKy0CZEUPAnNzCybCZEUwJPQzMyymNCT18zMbCgnBTMzy3NSMDOzPCcFMzPLc1IwM7M8JwUzM8tTRPFmk41N0lqgp95xjGA68GK9g8jIsVZHM8UKzRWvY901cyJiRqVCTZcUGpmkrojoqHccWTjW6mimWKG54nWsteHmIzMzy3NSMDOzPCeFsbWo3gGMgmOtjmaKFZorXsdaA+5TMDOzPNcUzMwsz0lhlCQdI+mXklZK+kyJ82dLWitpWfo4px5xprFcK+kFSctHOC9J30z/LE9IWlDrGAtiqRTrUZJeLvhcL6p1jAWx7CfpfkkrJD0l6YISZRris80YayN9tu2S/lPS42m8XyxRZjdJt6Sf7VJJc2sfaeZYG+b7ILOI8CPjA2gFfg3MAyYDjwMHFpU5G7i83rGmsRwJLACWj3D+OOAekq2qO4GlDRzrUcBd9f5M01j2BRakz6cA/1Xi70FDfLYZY22kz1bAK9PnbcBSoLOozF8BV6TPTwNuaeBYG+b7IOvDNYXRORxYGRGrImIrcDNwQp1jGlFEPAi8VKbICcB3I7EE2EvSvrWJbqgMsTaMiFgTEY+lz/uAFcDMomIN8dlmjLVhpJ/XxvSwLX0Ud3yeAFyXPr8N+GNJqlGIeRljbTpOCqMzE3i24Hg1pf+BnZQ2Gdwmab/ahLZTsv55GsURaVX9HkkH1TsYgLTp4jCSX4mFGu6zLRMrNNBnK6lV0jLgBeBHETHiZxsRA8DLwLTaRpnIECs0z/cB4KQwWqV+jRT/MrgTmBsRfwjcx45fNI0oy5+nUTxGMk3/EOAy4Ad1jgdJrwRuBy4rrIxQAAAFdUlEQVSMiA3Fp0u8pW6fbYVYG+qzjYjtEXEoMAs4XNLBRUUa5rPNEGszfR8ATgqjtRoozPSzgOcLC0TEuojYkh5eBSysUWw7o+Kfp1FExIZcVT0i7gbaJE2vVzyS2ki+ZG+IiDtKFGmYz7ZSrI322eZExHrgAeCYolP5z1bSJOBV1LnpcaRYm+z7AHBSGK1HgfmSXitpMkkn1+LCAkXtxseTtOE2qsXAn6UjZTqBlyNiTb2DKkXSPrl2Y0mHk/zdXVenWARcA6yIiK+PUKwhPtsssTbYZztD0l7p892Bo4FfFBVbDPx5+vxk4CeR9urWUpZYm+z7AIBJ9Q6gmUTEgKTzgR+SjES6NiKeknQJ0BURi4GPSzoeGCD59XJ2veKVdBPJyJLpklYDXyDpDCMirgDuJhklsxLYBHyoPpFmivVk4C8lDQCbgdPq8UWQeitwFvBk2p4M8DlgNjTcZ5sl1kb6bPcFrpPUSpKcvh8RdxX9G7sG+J6klST/xk5r4Fgb5vsgK89oNjOzPDcfmZlZnpOCmZnlOSmYmVmek4KZmeU5KZiZWZ6TgjUMSdMKVpP8raTnCo4nV/G+D0k6NH3+Q0lTypT9G0nt1YqlnDTOX0o6ruD40F285gHp57t+bKK0Zud5CtYwImIdkPtyvhjYGBFfKyyTTrJSRAxWKYb3VCjyN8C1QH817p/BqRGxrHKxbCLil5I6gBfH6prW3FxTsIYnaX9JyyVdQbJOz36Fv2wlnSbp6vT5ayTdIakrXeu+s8T19pB0a7pI2c1Ae8G51ZL2kjQlXRzu8fTeJ0v6BLA38DNJ96XlF6X3ekoF+xCk17lY0s/T+7w+fX2KpOskPZm+/v709WMlPSLpMSV7BbxiFJ9Pq6Tr0/tNkrRe0jfSa/1I0rS03Osl/ST9Mz2mOu1DYI3NScGaxYHANRFxGPBcmXLfBL4aER3AnwJXlyhzPtCbLlL2FZKVQ4sdB/wmIg6JiINJVsD8BslqmG+PiKPTcp9J73UI8G5JBxZc43dpvFeT1DAALgbWRsSb0vf8VNLewGeAP46IBcATwLDNcEYwCbgReDIiLk5fexWwJL3WI8Dfpa/fBHwjXfjuj9I/i9kQbj6yZvHriHg0Q7mjgQO0Y3n9qZJ2j4jNBWWOBL4KEBE/l/RUies8AfyjpH8E7oyI/xjhfqdL+jDJv6U/IEleT6fncovPdZMkmVx870/vHUBvWls4EHg4jXsy8FCGPyskSz7cGBFfKXhtALg1fX49cKOkqcD0iLgzvXe9mr+swTkpWLP4fcHzQYYun1zY8Svg8HQTpHLKru8SESvStvbjgH+SdFdEfLmwjKT5JL/oD4+I9ZKuL4oltzrmdnb8W1OJewu4NyLOqhBzKf9BssnMPxesxll8/RjhdbNh3HxkTSftZO6VNF9SC/CBgtP3AeflDkYYnfMgcGZ6/hBg2KYykmaSdHR/D/g6yVahAH0k21oC7Jkeb0hXw6zUSQ3w7yTNV7l9nKcCDwPvkDQvff0VacLJYhHJn/lmJctIQ7KQ4Inp8zOAhyKiF3hR0vvSe7RL2iPjPWwCcVKwZvW3wL3Aj0nW1885D3hr2on7NPCREu+9HJgm6QngE0BXiTKHAI8qWVn000CulrAIuC/taH6MpKloOcla+SM1MRX6IvAaScuBZST9E78DPgzcIulxkiTx+gzXAiAivprG8a8k/6ZfBhZIegx4G/CltOiZwCfTP/dDwIys97CJw6ukmjUJSQ8B55cbkprWFl6MiL1Gcd1Rv8fGL9cUzJrHOuD63OS1sSDpAJKa0u/G6prW3FxTMDOzPNcUzMwsz0nBzMzynBTMzCzPScHMzPKcFMzMLM9JwczM8v4HzXOUMtB+VJoAAAAASUVORK5CYII=\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": "markdown",
   "metadata": {},
   "source": [
    "*4. Using PyMC3, infer the distance to each star with the hierarchical model.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "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:17<00:00, 2280.97draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as hierarchical_model:\n",
    "    \n",
    "    # Hyperprior over L\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": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEoFJREFUeJzt3X+MXeed1/H3h5QEtyx10kyjME5xoBaoi1A3MknB0io0sM2Pqi5SIlIW6i2WzEopFBbY2MsfWcEWueJHuhW7QabJNi2lbghbxVqs7VpxqxWIZONsorRptmSUZpNxTOxV4vAj3u66/fLHfYbejueX587MnZnn/ZKse85znnvvM4+u7+ee55zznFQVkqT+/JFxN0CSNB4GgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTbxl3AxZy5ZVX1vbt28fdDEnaUJ588snfq6qJxeqt6wDYvn07J06cGHczJGlDSfK7S6nnEJAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ1aNACSPJDkdJJvDpX9yyS/k+SZJF9JsnVo24EkU0m+neQDQ+U3t7KpJPtX/k+RJF2MpewBfA64eVbZMeDPV9VfAP4HcAAgyXuAO4Efbc/55SSXJLkE+CXgFuA9wEdaXUnSmCx6JXBV/WaS7bPKfmNo9THg9ra8GzhcVd8FvpNkCri+bZuqqhcAkhxudb81UuulZdh18Dgnz55bcv3JrVv4b/vfv4otksZjJaaC+DvAl9vyJINAmDHdygBenlV+wwq8t3TRTp49x4sHb1ty/e37/8sqtkYan5ECIMk/Bc4DX5wpmqNaMfdQU83zmvuAfQDvete7RmmeOjffL/3JrVvG0Bpp/Vl2ACTZA3wQuKmqZr7Mp4FrhqptA15py/OV/5CqOgQcAti5c+ecISEtxcX+0pd6s6zTQJPcDNwNfKiq3hzadAS4M8llSa4FdgC/BTwB7EhybZJLGRwoPjJa0yVJo1h0DyDJl4AbgSuTTAP3MDjr5zLgWBKAx6rqp6vq2SQPMTi4ex64q6q+117n48BXgUuAB6rq2VX4eyRJS7SUs4A+Mkfx/QvU/yTwyTnKjwJHL6p1kqRV45XAktSpdX1HMGkxC53T79k+0sIMAG1onukjLZ9DQJLUKQNAkjplAEhSpzwGIC1icuuWOecDcpI4bXQGgLSI+b7knSROG51DQJLUKQNAkjrlEJA2BKd2llaeAaANwQu+pJXnEJAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ3yOgCtK17wJa0dA0Drihd8SWvHISBJ6pR7ANIyeZ8AbXQGgLRM3idAG51DQJLUKQNAkjq1aAAkeSDJ6STfHCq7IsmxJM+3x8tbeZJ8JslUkmeSXDf0nD2t/vNJ9qzOnyNJWqql7AF8Drh5Vtl+4NGq2gE82tYBbgF2tH/7gPtgEBjAPcANwPXAPTOhIUkaj0UDoKp+E3htVvFu4MG2/CDw4aHyz9fAY8DWJFcDHwCOVdVrVfU6cIwLQ0WStIaWewzgqqo6BdAe39nKJ4GXh+pNt7L5yi+QZF+SE0lOnDlzZpnNkyQtZqUPAmeOslqg/MLCqkNVtbOqdk5MTKxo4yRJP7DcAHi1De3QHk+38mngmqF624BXFiiXJI3Jci8EOwLsAQ62x0eGyj+e5DCDA75vVNWpJF8F/sXQgd+fAA4sv9na6Jz0TRq/RQMgyZeAG4Erk0wzOJvnIPBQkr3AS8AdrfpR4FZgCngT+BhAVb2W5J8DT7R6/6yqZh9YVkec9E0av0UDoKo+Ms+mm+aoW8Bd87zOA8ADF9U6SdKq8UpgSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlDeFl1bY5NYtc94YfnLrlnlvJC+NgwEgrbD5vuTnCgVpnBwCkqROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnnAxOq2rXweOcPHvugvLJrVvG0BpJwwwAraqTZ8/x4sHbxt0MSXMYaQgoyT9M8mySbyb5UpI/luTaJI8neT7Jl5Nc2upe1tan2vbtK/EHSJKWZ9l7AEkmgb8PvKeqziV5CLgTuBW4t6oOJ/l3wF7gvvb4elW9O8mdwKeAvzHyXyBtEN4oRuvNqENAbwG2JPlD4K3AKeD9wN9s2x8Efp5BAOxuywAPA/82SaqqRmyDtCF4oxitN8seAqqqk8C/Al5i8MX/BvAkcLaqzrdq08BkW54EXm7PPd/qv2P26ybZl+REkhNnzpxZbvMkSYtYdgAkuZzBr/prgT8JvA24ZY6qM7/ws8C2HxRUHaqqnVW1c2JiYrnNkyQtYpSDwH8V+E5VnamqPwR+FfjLwNYkM0NL24BX2vI0cA1A2/524LUR3l+SNIJRAuAl4H1J3pokwE3At4CvAbe3OnuAR9rykbZO237c8X9JGp9RjgE8zuBg7m8D32ivdQi4G/iZJFMMxvjvb0+5H3hHK/8ZYP8I7ZYkjWiks4Cq6h7gnlnFLwDXz1H394E7Rnk/SdLKcS4gSeqUASBJnTIAJKlTTgYnjZlTRGhcDACtCKd9Xj6niNC4GABaEU77LG08HgOQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpU94PQFqnvFOYVpsBIK1T3ilMq80hIEnqlAEgSZ0aKQCSbE3ycJLfSfJckr+U5Iokx5I83x4vb3WT5DNJppI8k+S6lfkTJEnLMeoxgF8Efr2qbk9yKfBW4OeAR6vqYJL9wH7gbuAWYEf7dwNwX3vUBrLr4HFOnj13Qfnk1i1jaI2kUSw7AJL8CeDHgZ8CqKo/AP4gyW7gxlbtQeDrDAJgN/D5qirgsbb3cHVVnVp267XmTp49x4sHbxt3MyStgFGGgP40cAb4lSRPJflskrcBV818qbfHd7b6k8DLQ8+fbmU/JMm+JCeSnDhz5swIzZMkLWSUAHgLcB1wX1X9GPB/GQz3zCdzlNUFBVWHqmpnVe2cmJgYoXmSpIWMEgDTwHRVPd7WH2YQCK8muRqgPZ4eqn/N0PO3Aa+M8P6SpBEsOwCq6n8CLyf5s63oJuBbwBFgTyvbAzzSlo8AH21nA70PeMPxf0kan1HPAvp7wBfbGUAvAB9jECoPJdkLvATc0eoeBW4FpoA3W11J0piMFABV9TSwc45NN81Rt4C7Rnk/SdLK8UpgSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqe8J7C0wXizeK0UA0Bz8sYv65c3i9dKMQA0J2/8Im1+HgOQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CmngpA2CSeJ08UyAKRNwknidLEcApKkThkAktSpkQMgySVJnkrya2392iSPJ3k+yZeTXNrKL2vrU2379lHfW5K0fCuxB/AJ4Lmh9U8B91bVDuB1YG8r3wu8XlXvBu5t9SRJYzJSACTZBtwGfLatB3g/8HCr8iDw4ba8u63Ttt/U6kuSxmDUPYBPAz8LfL+tvwM4W1Xn2/o0MNmWJ4GXAdr2N1p9SdIYLDsAknwQOF1VTw4Xz1G1lrBt+HX3JTmR5MSZM2eW2zxJ0iJG2QPYBXwoyYvAYQZDP58GtiaZub5gG/BKW54GrgFo298OvDb7RavqUFXtrKqdExMTIzRPkrSQZV8IVlUHgAMASW4E/nFV/WSS/wTcziAU9gCPtKccaev/vW0/XlUX7AFobe06eJyTZ89dUD65dcsYWiNpLa3GlcB3A4eT/ALwFHB/K78f+EKSKQa//O9chffWRTp59hwvHrxt3M2QNAYrEgBV9XXg6235BeD6Oer8PnDHSryfJGl0XgksSZ1yMjhpk5tvltCZbc4U2i8DQNrkFvqCd6bQvjkEJEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjrlaaCdcM4fSbMZAJ1wzh9JsxkAUsfmu0rYK4T7YABIHZvvS94rhPvgQWBJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnXIqiE3GWT8lLZUBsMk466ekpVr2EFCSa5J8LclzSZ5N8olWfkWSY0meb4+Xt/Ik+UySqSTPJLlupf4ISdLFG2UP4Dzwj6rqt5P8CPBkkmPATwGPVtXBJPuB/cDdwC3AjvbvBuC+9ihpnXGa6D4sOwCq6hRwqi3/7yTPAZPAbuDGVu1B4OsMAmA38PmqKuCxJFuTXN1eR9I64jTRfViRs4CSbAd+DHgcuGrmS709vrNVmwReHnradCuTJI3ByAGQ5I8D/xn4B1X1vxaqOkdZzfF6+5KcSHLizJkzozZPkjSPkQIgyR9l8OX/xar61Vb8apKr2/argdOtfBq4Zujp24BXZr9mVR2qqp1VtXNiYmKU5kmSFjDKWUAB7geeq6p/M7TpCLCnLe8BHhkq/2g7G+h9wBuO/0vS+IxyFtAu4G8D30jydCv7OeAg8FCSvcBLwB1t21HgVmAKeBP42Ajv3bX5LvYCL/iStHSjnAX0X5l7XB/gpjnqF3DXct9PP+DFXhoXTw/dXLwSWNKSeXro5uJkcJLUKfcAJI3MoaGNyQCQNDKHhjYmh4AkqVMGgCR1ygCQpE55DGAd8+5eklaTAbCOecGXpNVkAEhaNZ4eur4ZAJJWjaeHrm8eBJakThkAktQpA0CSOmUASFKnDABJ6pRnAUlac54euj4YAOuAV/yqN54euj4YAOuAV/xKGgcDQNK6Md/Q0EL1HTJaPgNgDTnUIy3sYr/MHTIajQGwhhzqkbSeGACrwF/6kjYCA2AV+EtfWhueTjoaA0DShjXfl/yug8cNhiVY8wBIcjPwi8AlwGer6uBat2GlONQjrU9eZ7A0axoASS4Bfgn4a8A08ESSI1X1rbVsx0pxqEfaWBY6zbTHvYO13gO4HpiqqhcAkhwGdgPrOgD8pS9tDgt9wc83bDSfzRAYax0Ak8DLQ+vTwA1r3IZ5v9DnM7l1i7/0pU3uYr/MLzYwLtZaBMxaB0DmKKsfqpDsA/a11f+T5Nur3qpF/C6QA3NuuhL4vTVtzMZkPy2N/bQ0XfTTAt87S/GnllJprQNgGrhmaH0b8Mpwhao6BBxay0YtV5ITVbVz3O1Y7+ynpbGflsZ+WjlrfT+AJ4AdSa5NcilwJ3BkjdsgSWKN9wCq6nySjwNfZXAa6ANV9exatkGSNLDm1wFU1VHg6Fq/7yrZEENV64D9tDT209LYTyskVbV4LUnSpuM9gSWpUwbAIpJckuSpJL/W1q9N8niS55N8uR3MJsllbX2qbd8+znavtTn66XNJvpPk6fbvva08ST7T+umZJNeNt+VrK8mLSb7R+uREK7siybH2mTqW5PJW3m1fzdNPP5/k5NBn6tah+gdaP307yQfG1/KNxQBY3CeA54bWPwXcW1U7gNeBva18L/B6Vb0buLfV68nsfgL4J1X13vbv6VZ2C7Cj/dsH3LeGbVwv/krrk5lTGfcDj7bP1KNtHeyr2f0Eg/97M5+powBJ3sPgjMIfBW4GfrlNO6NFGAALSLINuA34bFsP8H7g4VblQeDDbXl3W6dtv6nV3/Rm99MidgOfr4HHgK1Jrl7VBq5/w5+d2Z8p+2pxu4HDVfXdqvoOMMVg2hktwgBY2KeBnwW+39bfAZytqvNtfZrB9BYwNM1F2/5Gq9+D2f0045Nt6OLeJJe1srmmA5mkHwX8RpIn21XvAFdV1SmA9vjOVt5zX83VTwAfb5+pB2aGyui7n0ZiAMwjyQeB01X15HDxHFVrCds2rXn6CeAA8OeAvwhcAdw985Q5XmbT99OQXVV1HYPhnbuS/PgCdXvuq7n66T7gzwDvBU4B/7rV7bmfRmIAzG8X8KEkLwKHGQz9fJrBbvjM9RPDU1n8/2ku2va3A6+tZYPH5IJ+SvIfqupUG7r4LvAr/GCXfNHpQDazqnqlPZ4GvsKgX16dGdppj6db9W77aq5+qqpXq+p7VfV94N/jZ2pkBsA8qupAVW2rqu0MDjAdr6qfBL4G3N6q7QEeactH2jpt+/Hq4CKLefrpbw19oYXBmPY321OOAB9tZ7i8D3hjZvhjs0vytiQ/MrMM/ASDfhn+7Mz+THXXV/P106zjH3+dH/5M3dnOxLuWwUHz31rLNm9U3hLy4t0NHE7yC8BTwP2t/H7gC0mmGPzyv3NM7VsvvphkgsHu+dPAT7fyo8CtDA7UvQl8bDzNG4urgK+0cwPeAvzHqvr1JE8ADyXZC7wE3NHq99pX8/XTF9rpxAW8CPxdgKp6NslDDO4rch64q6q+N5aWbzBeCSxJnXIISJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktSp/wfI5Z/J6/xSywAAAABJRU5ErkJggg==\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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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 or emcee, where you don't know, or don't trust, the uncertainties on the parallaxes.*\n",
    "\n",
    "*3) Reproduce the file \"parallaxes.csv\". How would you sample from the prior distribution?*"
   ]
  }
 ],
 "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
}