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": [
    "import pymc3 as pm\n",
    "\n",
    "with pm.Model() as simple_model:\n",
    "        \n",
    "    # Prior over distance as a pymc3 potential object\n",
    "    \n",
    "    # The model\n",
    "\n",
    "    # Likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the resulting measured distances against the true distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*4. Using PyMC3, infer the distance to each star with the hierarchical model.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the resulting measured distances against the true distances."
   ]
  },
  {
   "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
}