notebooks/5_posterior_SHMF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "import logging\n",
    "import numpy as np\n",
    "from scipy.stats import uniform, norm\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.legend_handler import HandlerPatch\n",
    "%matplotlib inline\n",
    "\n",
    "sys.path.append('../')\n",
    "import paper_settings\n",
    "from simulation.units import *\n",
    "from simulation.prior import get_grid\n",
    "from inference.bayes import Posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "paper_settings.setup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "n_events_for_expectation = 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"calibrated_alices_full_sgd1e2_grid\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "llr = np.load(\"../data/results/llr_{}.npy\".format(filename))[:,:n_events_for_expectation]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 25\n",
    "\n",
    "f_sub_1d = np.linspace(0.001, 0.200, resolution)\n",
    "beta_1d = np.linspace(-2.5, -1.5, resolution)\n",
    "\n",
    "theta0, theta1 = np.meshgrid(f_sub_1d, beta_1d)\n",
    "theta_grid = np.vstack((theta0.flatten(), theta1.flatten())).T\n",
    "\n",
    "bin_size = f_sub_1d[1] - f_sub_1d[0]\n",
    "alpha_edges = np.linspace(f_sub_1d[0] - bin_size/2, f_sub_1d[-1] + bin_size/2, resolution + 1)\n",
    "bin_size = beta_1d[1] - beta_1d[0]\n",
    "beta_edges = np.linspace(beta_1d[0] - bin_size/2, beta_1d[-1] + bin_size/2, resolution + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate posterior on f_sub, beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def narrow_prior(thetas):\n",
    "    f_sub = thetas[:,0]\n",
    "    beta = thetas[:,1]\n",
    "    \n",
    "    return uniform(0.001, 0.199).pdf(f_sub).flatten() * norm(-1.9, 0.1).pdf(beta).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayes = Posterior(llr, theta_grid)\n",
    "posterior_narrow = bayes.posterior_based_on_expected_llr(n, narrow_prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drawing from posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_fsub_beta(posterior):\n",
    "    posterior /= np.sum(posterior)\n",
    "    \n",
    "    # Draw bin\n",
    "    i_bin = np.random.choice(list(range(625)), p=posterior)\n",
    "    \n",
    "    # Draw f_sub and theta uniformly within bin\n",
    "    f_sub = theta_grid[i_bin][0] + (f_sub_1d[1] - f_sub_1d[0]) * (np.random.random() - 0.5)\n",
    "    beta = theta_grid[i_bin][1] + (beta_edges[1] - beta_edges[0]) * (np.random.random() - 0.5)\n",
    "    \n",
    "    return f_sub, beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SHMF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def M_200_sigma_v(sigma_v, scatter=False):\n",
    "    \"\"\"\n",
    "    Relate central velocity dispersion to halo virial mass\n",
    "    From https://arxiv.org/pdf/1804.04492.pdf\n",
    "    \"\"\"\n",
    "    a = 0.09\n",
    "    b = 3.48\n",
    "    if scatter:\n",
    "        sigma_log10_M_200 = 0.13  # Lognormal scatter\n",
    "        log10_M_200 = np.random.normal(a + b * np.log10(sigma_v / (100 * Kmps)), sigma_log10_M_200)\n",
    "    else:\n",
    "        log10_M_200 = a + b * np.log10(sigma_v / (100 * Kmps))\n",
    "    return (10 ** log10_M_200) * 1e12 * M_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shmf_by_mhost(fsub, beta, msubs):\n",
    "    \"\"\" Evaluates dn / dm for a given set of m (msubs)\"\"\"\n",
    "    \n",
    "    sigma_v = np.random.normal(225, 50)\n",
    "    mhost = M_200_sigma_v(sigma_v * Kmps)\n",
    "    \n",
    "    m_min_calib=1e6 * M_s\n",
    "    m_max_calib=0.01 * mhost\n",
    "    M_0=M_MW\n",
    "    m_0 = 1.e9 * M_s\n",
    "                \n",
    "    alpha = fsub * ((2 + beta) * M_0 * m_0 ** beta) / (m_max_calib ** (beta + 2) - m_min_calib ** (beta + 2))\n",
    "    \n",
    "    return np.where(msubs > m_max_calib, 0., alpha * (mhost / M_0) * (msubs / m_0)**beta / mhost)\n",
    "                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find mean and credible intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_and_ci(samples, levels=[0.68, 0.95]):\n",
    "#     n_samples = samples.shape[0]\n",
    "    mean = np.mean(samples, axis=0)\n",
    "    all_down = [np.quantile(samples, (1. - level)/2, axis=0) for level in levels]\n",
    "    all_up = [np.quantile(samples, 1. - (1. - level)/2, axis=0) for level in levels]\n",
    "    all_cis = np.array([all_down, all_up])\n",
    "    \n",
    "#     # Loop over data points\n",
    "#     all_cis = []\n",
    "    \n",
    "#     for samples_, mean_ in zip(samples.T, mean):\n",
    "        \n",
    "#         # Find position of mean in samples\n",
    "#         sorted_ = sorted(samples_)\n",
    "#         mean_pos = np.searchsorted(sorted_, mean_)\n",
    "    \n",
    "#         # Loop over levels and find lowest and highest position\n",
    "#         cis = []\n",
    "#         for level in levels:\n",
    "#             half_n_ci = int(round(level * n_samples / 2,0))\n",
    "#             low_pos, high_pos = mean_pos - half_n_ci, mean_pos + half_n_ci\n",
    "#             if low_pos < 0:\n",
    "#                 cis.append([sorted_[0], sorted_[2*half_n_ci]])\n",
    "#             elif high_pos >= n_samples:\n",
    "#                 cis.append([sorted_[-2*half_n_ci - 1], sorted_[-1]])\n",
    "#             else:\n",
    "#                 cis.append([sorted_[low_pos], sorted_[high_pos]])\n",
    "                \n",
    "#         all_cis.append(cis)\n",
    "#     all_cis = np.array(all_cis)  # Shape (n_datapoints, n_cis, 2)\n",
    "#     all_cis = all_cis.T  # Shape (2, n_cis, n_datapoints)\n",
    "\n",
    "    return mean, all_cis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate SHMF data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_realizations = 20000\n",
    "n_realizations_show = 0\n",
    "fsub_true, beta_true = 0.05, -1.9\n",
    "msubs = np.geomspace(1.e7, 1e12, 200) * M_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "shmf_posteriors = []\n",
    "shmf_trues = []\n",
    "for _ in range(n_realizations):\n",
    "    fsub_, beta_ = draw_fsub_beta(posterior_narrow)\n",
    "    shmf_posteriors.append(shmf_by_mhost(fsub_, beta_, msubs))\n",
    "    shmf_trues.append(shmf_by_mhost(fsub_true, beta_true, msubs))\n",
    "shmf_posteriors = np.asarray(shmf_posteriors)\n",
    "shmf_trues = np.asarray(shmf_trues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "shmf_true, _ = mean_and_ci(shmf_trues, levels=[])\n",
    "shmf_posterior, shmf_posterior_cis = mean_and_ci(shmf_posteriors, levels=[0.68, 0.95])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 1./M_s * msubs\n",
    "y_mult = msubs**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANYAAADTCAYAAAD9Lu2dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXtwW9d95z/nAiDBB0jwLZEiJYF6y/GDohI/4lhxqG0ax3k4dJU26abNNFKz02l32llzPbN/uH9svdI2s5uZzrSWu51Jt9lENlM3cZqkER3bTRPbtUi/LUuUIInvNwDi/bj37B8AaIriAwTvxYO8nxlJuNDF/R0A94tzzu/8fr8jpJSYmJjoi5LvBpiYbEZMYZmYGIApLBMTAzCFZWJiAKawTEwMwBSWiYkBmMIyMTEAU1gmJgZgCsvExACs+W7AagghHgYedjgc39i3b1++m2Oyxejv75+RUjZk81pRDCFNnZ2d8sKFC/luhskWQwjRL6XszOa15lDQxMQATGGZmBiAKSwTEwMwhWViYgCmsExMDMAUlomJAZjCypIzZ87Q09NDb28vQgjOnj3LmTNnOHXqlG7X7+3tpbe3lzNnzgAwMDBAe3s7p06dwuv10tfXR01NDWfOnGFgYIDjx48vnJumr6+P9vZ2+vr6Fh6n233mzBmOHDmiS3tNliClLPg/R44ckYXGU089tfDY5XIt+3y2nD9/Xp4/f37h+OTJkwuPu7u7ZX9//8JxR0eH9Hg8C6/r6Oi45VpdXV0rvv706dMbbu9mBbggs7xnzR4rSzo7l183XOn59eD1ejl//vzC8Xp6QZfLhdvtXjiura1d1U5HR0d2jTRZFVNYWbLSDel2uzl+/PjCUKuvr4/jx48D0NPTQ09PD5C8qdP/f/bs2Zuu0d3dTV9fH0eOHOHMmTPruvlPnDhBb28vkBw6ulyuFc/t6+ujq6sr42ubZM6mEdYTTzzBE088AcC+ffu4fPky/f39C3OIP/uzP+Nb3/oWAM3NzYyNjfHSSy9x7NgxAE6ePLlwgzscDvx+f1bt6O7uxu12093dzcmTJ2+6cRf3PE8++SQdHR10dXVx9erVW67T39/P6dOnef3112/psc6dO7cw//J6vbfYX9zbLce5c+c4deoUc3Nz2bxFkwwo6CDc9ZAWFcDly5cXHvf39wMsiApgbGwMSArspZdeArip18hWVGnSPYzT6VzxHLfbjdfrZWBggLq6ulv+z+Vy0dXVRVdX10KPl+bEiRMLNp588slbru10OnG73SvaT7++r69vXe/LJHM2TY9VyKRv8MVzn6NHj+Jyuejo6ODkyZM3nd/X13fTuasJdDlOnDhBT0/PqsNAwBwGGsim6bHygdfr5ZlnnsHtdnPmzJmFYeDAwAADAwMLvcrRo0fp6+tbcJG73W4ee+wxzpw5szAcW3yT19bWLlxjbm6Oxx9/HGDhuXPnzuFyubhw4QJut5uzZ8/S0dFBT08Pjz/++EI7INkTu91u+vr6Fq577tw5nE7nmsIzyR4zbcTEZAXMtBETkwLDHAqaFCxSSmYjUa75AgRiCdqqKtAk7HCUUWYt7Fu3sFtnsiUJxOIM+UO8MjZNIJ5ASlCl5N0ZH5B8vK/GwYHaKloc5ZRaLHlu8a2YwjIpGELxBK+Oz/D2tBdFQFy7ef4f07SFx4MeP1e9ATQkD+xo4s7Gmlw3d1VMYZnknaiq8osbE1zyJNcPVSlR1/CpJaSElOPt5eFJvJEYH9/RgFUpDLfBphHWX795mVBC1e165VYL37xz5cpQAwMDPProo3R3d3P06FFef/11jh8/bsjakNvtpqenh2effVb3a+eTqKryxqSHf5+YQUsN97IhISVvTntQBHyitUnnVmbHphGWnqLK5HodHR10dHQsRDF0d3dTU1ODx+PJ6Pq9vb10d3dndK7L5eLpp5/O6NxiIKaqvD4xS//kHBJIaBtf8lGlZGDKg6II7t5en/eea9MI6/U3Zogn9FuTs1kFrHMFo7a2diEcaTXS0euZCisd/rQZItFvzAf5iXuUmKolh3M6okpJ/+QcU8EIj+xr0/Xa62XTCEtPUWVzPa/XuxDNcPbs2YX0jZMnTy4bk3fhwoWFXsvr9S5ET7jdbmpra3nqqac4derUQkBvT0/PQnDt0uv39vbedP5jjz2my2egJ3FV48XhSS7O+nQX1GISmmQ4EGI8EGZ7ZZlhdtZi0wjL/dTf6X/Rb35izVP6+vqYm5vD7XbzwgsvLKR5dHV10dfXx5kzZ5idnV2Yf6V7tNra2oUe68knn1z4/56eHk6ePElPT8+C6JxO50K84HLXf+yxx246v5CQUvLerI9fDE0ipTRUVGkSmuTXY9N8KY+91qYRluvU1/NiN32Tp3n99dcXBONyuXjqqad4+umnefLJJ+np6Vl2rrRcpPtKEfLLXX+18/OJJxLjZ9fGmApHdJlHrYdhf4j5aJyqUltO7aYpDN+kDtisoiCutziD1+12LwTgnj59mv7+/oVhYVoAAwMDq0a6Z3L9QkPKZI/xnffcjAfDORdVqhW8OZW/fLNN02Mdvas+p/YGBgZwu908++yzdHZ2Lgjl9OnTNxV/SQ/T0izubXp7e+nq6rol0j392nSEfPqx2+1e9vp9fX23RNTnC180xg+vjOCJxrJ2n+uBKuGtaS+d2+oot+X+Nt800e3funBRd7t/1nlQ92tuZkb8IZ4bHCauaRTCXaUAOxzlPLp/Z1av30h0+6bpscqtFt0XiE0yI65qvDA0waW5+Zw4JzJFA0YDYXzRGNWlJTm1vWmEtVqUhIlxhOIJnrl0A280nteh38ok17YebNuWU6t5d14IIVxCiA4hxGNCCDOltYiYCkX4+/ev4Ynkdz61GqqEd6a9zIajObWrm7CEEB1LjruFEF1CiLVWKzsAN9AHZBaKYJJ33p328L0PrhOMJ9DWPj2vqFLy3mxu1/d0EZYQogt4dtFxB4CUsg/wLhXdYqSUvVJKL9AF9OrRHhPjiGsaz18d4YXhyTy50dePBK54Azm1qcscS0rZJ4RwL3rqBJAubucmKZoBIcTSHqlPSulNCbNPSunGpGAp/PnUyviicSIJFXuOnFJGOS+cwOLVuTpI9k5LT0yJqgdwCyHOL3eOSf6ZDIb5weAw0YRa8EO/5bAKwbA/yN6aqtzYy4mVVUgNF2+JUhVCnAROArS15TdSeatz1evnx+7Rohn6LUdM03D7AjkTllFeQS+QrsbvBGbXewEp5VkpZaeUsrOhoUHXxplkhiYlvxyZ4sdXi1tUaa75gjmzZVSPdY4Ps5lcLNMjmRQ2kYTKc1eGmQpFCmrRdyOE4gkSmpaTJEi9vILdQGfaOSGlHEg93wV408cmxUFyfcrNZN4CaI3Bqij4ovHc2NLjIimHQ++S586ucLpJAXPNF+BHV0c2laDSCJEMEq4rKzXcVt6dFyaFgZSSgUkP/zY6tWmGfktRNYm3mHosoxBCPAw8vGfPnnw3ZVOjapK+oXE+KLAgWr1RpWQukpvQpoIWlpTyeeD5zs7Ob+S7LZsVfyzOc4PDeKKxvAz/pJTEYhqRqEosrpEsFygBgb1UoaLcitWqn7NhJkcxgwUtLBNjmQiG+cHlIWKqltNF32hMxeuLMTUTIRBMIAChiFT9TZmMQRIghEDTJNVVJbh2VmIv3XjUhCcS2/A1MsEU1hZESsmFyTl+PTqds6FfQtWYnokwPhkmGtMQAtIVoyVwS+lbufAXXl+Mt971sG9PFTXVG8uriiRUYqpKicH13k1hbTHimsaPr44y7A/mpmJSQmNoNMjUdAQWi2mdplVNcumKj1pnKXt2O1CU7GqSWBWFmXCMZoNLo5nC2kIEYnF6Lw/lJIhWVTVGJ8KMTYRYKLO+QZOaBnOeKO/HNQ7vr0aI9YtLk5KZcMQUlok+LATRqipG+ig0TTI1E+bGSAipSd1taRICwTiT0xG2Na5fHAkpmQhGuN3gKLmCFpbpbteH92e8nL8xYfjQzzcfY9DtJ6FqaAZ6QzQNrg8HqKywUlmx/rqBU6GIAa26mbyn5q+GlPJ5KeXJ6urqfDelKNGk5Pz1cc4PGSuqaEzl4qCPi4M+YnFjRZVG0+Di4DzqWvv9LIMnarxnsKB7LJPsmYtE+dGVEXzRuGGi0jTJ2ESIkbGQocPLlVBTjpHdbZXrel1C0wxPejSFtQm5OOvj5zfGUTVpWH0/33yMwWt+EgktL6KC5HxrcipMc1MZpetY47IKhdlwlBZHuWFtM4W1iQgnVH5+fYzr80HDoihiMRX3jQDe+VhOhnxroUkYGg2y15V5AmPSM2gKyyQDrnj8/Oz6GAlNW3Ob0WzQNMnEVJih0WBBCGoxM3NRWlvUjCMzElIyFgxxB8btW2wKq8hJaBp9Nya45Jk3rJea98cZvDZPPEeOifUiJQyNBNnXnnmvNR4w1jNY0MIy3e2rMxOO8NzgSDIz1gAHRSKhcW0owOxcNG/zqEyZ9UQJR1TK7Jn1Wr5YDE1KlCwWmTPBdLcXIcncqTm+e/E68zFjvH6znigDb88xUwSigmSvNTYZyvh8ixDMG5ibVdA9lsmtRBIqz18dYcygtPlYXGPQPY8/EC/IYd9qTM9E2NVaiSWDOEIhBP54HKfdmM0SshKWEEKVUprbceSYofkgP3aPElNV3R0UUkqmZyO4bwTRiqGLWgaBwOuNUVe7duq9lBCIJQxrS7Y9ljEDU5NliasaL49M8t6MMRtjRyIqV677CQSLr5dajKpJJqfDGQkroWn4Y4U3FCzOn7QiZCIY5odXRogkVN1FpWmS4bEQ4xP5iZwwAp8/TiyuUWJb3X0gwdD6F7rNsYQQXyJZQ/CqlPIf9bruVkWTktfGZ/j38VlDeql5f5zL7vm8Rk4YxcRkmLYdFWue5zUwZlBP58WAlPIHQojdOl5zSzIZDPOTa2OGePxicY3rQwHmPMXh7VsvUsL0bCQjYRXiUHA5XKleawC4puN1twxxVeNXY9O8NeUxZNg3NhFiZHxR4uEmJRZPFqdZKxIjFNdva92l6CksJ0mnhlOvC26lBeIb80F+4h4lpmq6i+qmgNkidk5kiiBZJ2OtRMiENC7KvaCHgluh/Nl8LM756+OMBEK6r0vF4hru6/6cBczK1A9C2DPF9MV+Glo/wuCLz2K3O6mqauWdX/4dh+/7XeLxEHGitHQew76tCcWq73KqJsGTgbCsQsETibHdgDT9db0jIcSdUso3V/hvcyi4DlQtGT3xyvg0CZ3TOxYHzErNWBducHaceMhP1XYXL/6Pk+w98kVikzN4Ji5TcpugtnEvJTW1VDQ0c3vjH+Goa8U3PEjwxiU8/W8w+v4vab7zAbbd/wDWyrXnRZniD6w9f5IkC3jmTVhCiP9CcvM4CbzJ8utYHSS37+kAfqFXAzcjY4EwP3GPEkokiOvcS+Vi2KepCd7u/Stu+9wfMtX/Kr7rg1TbtrPnwGdwNLVS3nkP7fV1txR7qSQ5mKnY1Ubz/Z9CSkmgxEtpdQ1Xv/v31Nx2Ow333IPQYTcQVZVEYyqlJSsP8+KaZDZsjGcw0x7r/JKearm7oY/ktqinN9yqTUownuCl4UmuePy6z6Oi0WSelM+v/7BPSokQghuv/ox4KMD23Z0oPpUb//A9yltaqLv/81TsasNStr5ffiEE+3/zqwDMzF1i4t1XSMx4afoPD677WktRFEEolFhVWABTYWOi3DMV1h8KISTQL6X8W5bvsWqllG8IIZ7Sr3mbA01K3p728K8jU2hS6hqOpKqSkXFjFnk1NQFS8sJ//zr3//G3KYmUEL8+j9f/LruPfY6KnTtRSvWJtTv42d9DSsngP/5fJs/+L277vT/CVp397ouaJglF1DUzrsIJY8KaMhXWUynRpN/pLV+hlPKF1BzsDf2aV9xIKbni9fPi8CSRhKrrsE9KydRMhBvDQTSpX5mxeDiIrayCqy/+gIh/DlfHQ+y/q5uJH/6EynYXbZ//AiV1tWtfKAuEEGy75358H3zA+E/+hZYvPozFbs/qWlJCMLi2aKKqMePlVYUlhPgG0A5IIcQJ4C7gN1jSY6XmYO2Ah+QcbMtz3RfgxeFJ/LG47vMo73wM9/UAsbiqy7BPS8SJhQOUVlTxiyf/gAf+9K9wVu0mOKEy+8prOA8fpHJPO0rJ+kuNrZeq5t1UNe9m4Km/IP7cM7hOfDXrOVcovLawYvkQFtAnpXxaCHFXqsf61ArnnZdS/k+9G1eMeKMxzl8fZzwY1l1QoXAC942ALsGyUkrmx69RtX0311/5KdH5OXYffYg7jp1i4rmfUt7WSuMnH6C0qTGrirMbpeXjXUQvjzL32uvU3fOxrK4Ria69ABw3yMOzqrCklGm3eXqoutKd8mUhxCmScYJ/qVfjiompUIRfjU5zYz45NNNTUtGoytBoUJdM3ojfQyzgw9HUxlvnvs3Rr/036ur345/6gKm+F6k6dICGez62YefBRmk6/FHC22d4++m/pKLdhb1x/aVrpUxmQa+2DVBCkwvOGT3JdI7lEUL8DfA3qeOlX+/3V1nfyppCj7yQUjLsD/HK2DQTwQiqzoKKxzWGx5IbCmxEUFJKxt/6N7Z95B7mrr1PYHIYu93JoY/+NpPP/xx7YyO1Rzsoa92Rl95pJUoqHVS07mDmV6/Q8oWH1902RQgiUZXKVYSliGSvpffuI0Jm4fZdmuiYEt1ir6GudHZ2ygsXLuh92ayRUnLZ4+flEf2dEgDxhMbIWIjJ6fCG4vpCc5MEZ8Zp2Hcnb5773+w//hVUb4D59z4g7vHiOLAXx8ED2BzrK3iZS6SUDH7n/7Dt3k9QtW/ful5rUWDP7qpV87OsiuDrt7XjWGb+KITol1J2rrvR6BfStNRruCnRpOT9WR+/Gp0mquovqISqMToWYnxqY4IafeNlmg4eJeKbxTd6hZrmvbS5PsHUP/dhq66m6tABKnbvRBi8R5QeCCHwRIaQv3wRx9696+q1NAmR2OrzLAVBVFVxoK9jZkPCWsVruGnQpOTGfJBLc/Nc8frRpNRfUAmN8ckwoxNhpJRZCSrsm8U3Msi2w3fjG3VTvWMPdms1lTEnI73PUbmnne0PfZqSWuNq6RnFR0/9OWP/9DwB93Uc7ZmHokoJ4Qw8g5GE/g6MjfZYmXoNiwopJaOBMO/Nerk850dijPcoHtcYHQ8xMR0GSVbzqPF3fk1N237UaBjv0CANrjto3nGEub5fo9hsVB06QOODD6DYjHeVG4UQgrnoEOM/6ufIf35iXb1WOLK2ZzCq6p8+siFhrcNrWPDEVI0b80EGPfO4fQE0KXUPjk0TiaiMToSYnokgWf+QLxrwMnP5TVo6jhGYGqGspgm7pYLakjaGv/8DynftpPHBByhtbCgoZ8RG2PuZrzD1sxcID41QvrM149dFo6v/IEqkIWtZehWT8Qgh/hoomnCmmKoyFggz4g/h9gWYi8SwCEHMoHUNKSX+QILR8RDe+VhWw73Ji69TUbcdW1klnqHLbL/tXhqbbmP+lbfxaxpVBw9Qd9/dWOxrF1MpNuxVNSjN1dzoe54DX/9mxj8YiTWGeZqESKH0WFJKZcnxG8A3dWmRzqhaMjVgLhLDF40xGggzFYoQSiSwCoW4pi30SkZsH6qqybJioxOhrEo0x8NBxt76JTvv/jQR7ww2ewU2pYymmkMMffcZynY0U3/v3dibt22a3mklSuprSbwdJjI6TtmO5oxeo8nkd2CxLP/ZaFISLcA5Vk5QU/vGSglCfDh0UqUkqiazQFUpUzuia0RVlclQBG80RjCeFBBwk4gAYtK47L9wOMHEdITJ6TDAugU1M/gWltIyqrbtZH7UjRqNUVu7h/k3LxIIvYPj4H52nHgEa7lxO2YUGo0HjmDHwdxbb9GSobAUBWJxlTLL8re6JFkKTW+KomDnTDjK9y7eWKExgEyOlTW5fK9jpIAWE49rTM9GmJyOEI2qSZf5Ol6vxqPceOWn7L7/80SDPmyaihoMs31bB8PffxZ7YyPOjjspb92hS85SMTJ86V8JvHeZ7V2fyiiyXghBLKZRtkosrxGOqaIp2GnU3GejaJpkzhtjYiqMPxBHsH7v3tz1i6ixKPV7bifkmSIRDlNVvh3/xcuMv3YJx4G9tDzy+YJeyM0V+3/zq0zZXiY0PELlHtea50spicZWv3f0Xj4BHQt2bqW6gpommffHmZqJMOeNAmKhLHOmX5HUNK682Ev7sUdQo2Hi0TCx2Tm21d7G6DPPYW9qpPr2w5S3tSEsW7N3Wg5riZ1AfIr4u8GMhKVpyc3yViOvPVZKOLtZuZ6F7sVk0rGCdTva9Lpk1kRjKl5fjKmZCIFgAkUI1IVfusx/8XwjVwl5Jtn+kXtR41Fifh8lcTvRy0NMvjuK48A+dnR/Qdf6D5sNa7WD0PvDSFXNKHpkrSj3fM+x1hKO7sVk0lWaWg/dlvMqTYmEhs8fx+uL4fHGiCc0hPjQCbFeD+KVXzzLrnsfQtNUEtEw4bFxam1tTDz3U8paW6i7+6PYW7Zves+eHriOfYFR748Jj45R3rb2mlYisfp3ZcSuLesRVlo4K2UI615XMJdomsQfSAppzhsjElVRFIG6KI9+vd54/+Qw3qFLtB7tQgiFiGcOZsPIQQ8z116l6sA+6u77WNZZslsVKSVvv/Yd7ij/g4yEFV/DnZ5vYa0lnKIpMS2lJBrVCITizPsTzPtjhCO3CknNsjiF+5c/YvtH7kUoColYhMAVN+VhBzM/e4kK1y4aP/UApQ31Zu+UJUIIPvEn32bm5y8hNW1ND+la32PCAK9xQQ8F9SAtolA4QSCUwDcfJxhKABIBLI5myVZIkEzRmHjvNVz3fw7FYiU8Po465UNc8+L3XqZy316ajj+IYiuKpcOCx1pZgc8/Qmx2jtKG+lXPXet7VQt8KJjXuoJpt2oonCAcVgkE44QjKpGICiK5KZkRH+DQa/+Cs20/pQ4najiEp/8NlBtBgsODOPbvYcejX8RaUXyOCEUkewYECxH36aGwEMkkQgRILbl+KESy5FgymDi7CP31IDWVqAwSHh1bW1hrfO9G7OayrLBSxWE8gHNRqv1i4SyH4XUFNU0mC95HVCJRNSmilHhi8aRzQQiBpi4JnpULf+lC1O/lxqs/Y9/xL6MoNoLXbhCYeR+bR0Mtj9D4qWMFPdRbLIz0MoHNqlBSolBiU7DbLZSWWCgtUbDZFKxWBatFYLGIhc94MVJK4glJJKKSUJPfjz8QxzsfR0oM2SHSXlXLoYd+n/n3L+G88/ZVz9UKqMdaTiRrCcewuoLhcIJ/H5ghoUoUJdn7LPerKHUW0FJG33gZe3Ud1c3tJHx+JvteJDE8hrV5G1UfOUR5W2tBJA8qAoQiUr1M8vMoKVGwl1oot1spLVUoLbUsiMhqERv6ERBCUGITt2z2lgw8jjM+FWHOE9W9F5saeZfxN15g26e7Vv3cNflh0dHlMCJGdCVhLSeSGimlTwhxFvgGgBDiQSnlLyBZVzD1r+51BVUNEqlfnaS7O3fZKfFIiMHz3+PAQ7+PFowwP3QR/y/fotq5Hfu2JurvuzuvhVcsSvJHRgiSwimz4qi0UWa3UJLucTYonGwRQlDlKKHKUUIonOCDQR/RmKabwJrv+jhiNEhkapqy7dtWaUeyV7KuEIibsx5rBZHMpZ4bWPQl/ZYQ4tHU435gbrNEXUy89ypSQk2Di8SMj6F/+B7W8gpq9uyn4kFXzsOLko4vgdQkthIlJSArleVWysutlNiUgh16ApSXWbnjcA3vX/YRDCZ0KTBaWulEq1AIXHOvLayEZKXdejQDfqjXdF6sVoxTSvmHKS/hHFALnAKKVliamuC9fzrL3gdOEL4+RmR0jLhjmB23PUDlHhclNblJa0/PY6SUlJdZqXLYcFRaKS+zYi+1JJ0ERYjFonB4v5M33/VkVPMvEzwzV9D8YRruvXfFcwSChCpZKUvNiDlgJl7BVYtxSimvpYrISOAvdGtZDpkZfIvg5AjVlW0kxuYY//E/49i3n+1H78+JE0IRLHguqxw2nNUlVFZYqaywFnQvlA2KIjiwt4q33/foUsX3zt/+U4b+37nV17PE6i53LYdzrMWsWIxTCLGLpJewnQ+3+Xlc5zYagpSSt77/bXbs/Tiei+8Q8cxQfkc1Bx75Ovbt2wxPy0jPjSrKrTTU2XFW2bDbLZtOSMtRXmZlV2sF14eDGxZXLOJneOQVmuceorS+bsXzEquk32us7tzIhkyEtVoxzrPAz4FnpZTXdWuVgUxfHGDqzVdx2lvRxnzE6uZoeaCLspZmwz16ipIc4tU4S6irKcVZZcOyRSPXmxrKmJyOEAxtbEhY6nCybf/HiExMriwsuXqPJUh6Bq05FtbSLXwW82jKU3iXEKId6JJSFlyPFZ+f5+1zf0VtxU5igQC26jJqOu+i+XMPGZ6SoSjJuUVDXSl1NaWbcniXDUIIdjRXMOie31CvpVis1B/swH/DTfVth5Y9R7J6xqkQKWFl34xbyORa6WKc1YvbAiCl9KX+fSPZQFEw5Wrjvnmm+l9h9I1/pbH+MBWVjdR0dFC9Z39OxGS1fiiminJTTMtR6yxJfS4bm+NMXu8nePkqOx76bFavV9JROToOWNYUVlo0aRGlHi97Zy4+Rw/S+VgVDS1rniulJDY7R+j6EIMv9+Io345z734aO+5m5wMPGT5nUpRkNENDvZ3GejsV5WZM4FoIIaivLWVyemO7Kh54+Pe48Z3vkgiGsFZkUQNE6B/hXtDffjofq2bngWXzsaSmERmfJHj9Bp4P3mV68j32f/LLbL/7GHUH76C8ttHQ9ikKIKGmppRtDXaqHDazZ1onjfV2pmejG3J5S01lcvY96qc+TuXuXet+vUDoHuFe0MJaDi0eJzw8SvD6DULDI8x6rlCz9zA7PvswZaN7qOv4GCv7hvRBEVBWZmV7Uxl1NaUrltYyWZvKCutC2FW2KBYr9pp6otPTWQpL/+iLohCW1DTmP7hM6NoNwuMTWGoqGBvt58jX/iuO6RuUVddTVtfEjqYmw9qQjnxorC9lW2MZ5WVF8dEVPEIIysusqVSe7Gl/4Iv4B91ZNkL/eMGsyjm8AAAOUUlEQVSiuDvicx7Cw6PMxyYo+0gzO458Estb9VgrKqhz3GaobUUBm02htbmC+trSoo16KGScVbYNC2vkg1/heftNtn+6a92vFWyxOVaaOFEaHryf0vHrWEvsCEWh5a5PGGYvPU2qqylle1OZ6SI3GIfDhmU6sqFE0/YHH2F4XEWNRrGUrq/EthFFO4tCWCUVyW23nDuM3dkxlaJEU4OdluaKW9IgTIyhsty64Xg9a6mdQGKW6Mws5S2ZVcldwIBkiaIQlq2sEmWFEsF6kBZUY4Od1uYKbAYJyioEihDENY0yq4USi4JVUXCUJN9bNJEsjx1VNeKahqpJNCmxKgoS/fflKhRKSiy31BtZN0Jhfn6YyOTk+oVlAEUhLKMQJEVlhKAEYFMUVClpKCul1VGB026jocxOU4U9mcGbAaommY1EmQpF8EVjhOIq87E4gXgCfyxOQpPJuENNGpJinivKyyz4A9nPs4QQ3P7wN4lMTOnYquzZksJK39MNdXZaW8opLdFnyV2QXCQut1rYW+Og3emgpbIcywYcHhZF0Fhup7H81hJpUkoC8QTz0TizkSgj/hCzkRjeSAxNShSRDPQtBsFVVtg2JCyAqeG3CVwapPGT9+vUquzZcsJSRHJBd3drBSU6CcoqBBJod1ZyX0sDtTnan0oIgaPEhqPERoujnNsbkvliUkr8sQSzkSjz0TjeaIzRQIjpUBRJUqxx1Yj0vuwps1tQRHa7WqZpvO0oiSsTGVfINZItIyxFSX557bscVFbos22oRQisiqCzqY47G2uwr5SimmOEEFSV2qgqvfl9SikJJVSmQxFGU5vuzUWiRFUNRRizOUCmlJYoCEXABuZZjm1t2CqriHm8q6aQ5IJNL6x0dLlrZ+WioM+NYxWCfbUOHmzdRmmBCGothBBU2KxUVFeyq/rD0gKheIJrvgBvTXuZCkXyIrLSUsuGa2Go8SiX3/0hDR+9xxSWUSTT26G1uYLtTWW6LewKwGZR+MzuZtqdDl2umW/KbVYO1zs5XO9kPhpn2B/kkmeeoflQzkRWWmLRweVexj1f+XNis3M6tWoDbcl3A4xAEVBXW8qu1kpdPX1WReCw2fjSvlaqM9j0rBipKrVxuDQpskhC5f1ZL78em0GTxrr7LRaBRdl4UVXv7DVi47PU33e3Ti3Ljk0lrOQ8ysqe3Q7d0zasQnBHQw33tzRuyMtXTNitFjqa6rijoZZ3Z7y8PDJpqLhsNoEa3WCvVVFJOLDC7p85ZFMIK13ZdVdbBY31dl3DjwRgVRQ+62rGtUmGfuvFogjuaKyhubKMH7tH8UXjxhS5tCoQ3VhoUcvRY1x/ZxgtkUCx5u/2LvqYHUVAfW0pHXfU0tRQpm9BECGoLrXxu4d2b1lRLaah3M5XD+2mpbJM1/oQaWzWjd+OgekRrrh/TmLer0OLsqdohaUoSRftof1O9rqqdPlSFmNTBC5nJV877KLGvjnnU9lgUxQe2dtGi6MMi87iKinZ+HdYUd/M/o8+StyrazL7uilKYSlK0tvXcXstVQ591qQWYxWC+3c08llXC9Ytujv9algUwefad1Cic+0QPRxNisWKqsQJT0/q0KLsKao5lqJAtaME165K3cKQlmJL3TSL13lMbqXEYuFYaxN9NyZ02xzbZlUQYv07Zy5lauQdBBZW39zHWIrm57jEpnBgbzUH91UbIiqLgEqbld85uMsUVYbsr6nCqqOH1GYVGQcnr8btX/hP2LT8Dt8LWlhCiIeFEGdlLEjH7bU4q4z5sKxCUF9m52uHXdSXmfsBZ4pFEXxmd7Nu4rJalVRhvY3hmRxk4srrG7/QBihoYUkpn5dSnnTUVBuWEm9VBK1V5Xz5wM6CifUrJnZVV3KwtloXL6HVKnRJOqzc1ordXo0aiW78YllS0MIyGquSXPT9wp5W00mxAY61NqGHg9BqTSZ0bpSq7btwNO0k7sufZ3DL3k1WIXhgRyPHWpt0GddvZUosCh2NtWy0CpzNqqBHeb94OMhbv/pb4nlcy9qSwrIKwfFd27izsTbfTdk07HZWYhEbu50URZ/yE7ayCu75rSdIBAI6XC07tpywrIrg83t2cKjOme+mbCoay+y6VzraCH7fCIHxobzZL6p1rI2gkByydO9ro6kif3sGb1ZsFoXKEhvzsXi+mwJANDpPIpi/If6W6LEsqYza/3jYZYrKQNqqyvXwluuC69gjlCoVebO/6YVlFcliLF89tBtHif7hTyYfcnt9ja4LxhvBO3mVoUsvIfM0PN3UwrIpgt3VlZzYv5PSPBcX2Qpsq7BjK5Bli+qWdhqab0MNhvJivzA+BQOwCsGRploebm/ZMomJ+UYIQWVJYUzbSyqqsFRWEJufz4v9TSksqxB07dzGfS2NZs31HFNhKwxhCSG4+u4/E52bzYv9wvgUdMSmKHx+zw52VuVv4rqVcdgKZx579Is9yEh+vJSbpscSgN2i8JWDu0xR5RFHia1gPINzE5fwDl/Oi+1NISyF5BDkq4d2U1eWmyq0JstTbrPonlmcLdbKStRQOC+2i15YCqTWqHZv2pJkxUS5zVowsZctnZ+khPysWxa1sCxC4LSX8DsHd1GWx4o8Jh9SXkDfw+zQe1z74DxSVXNuu3A+hXViFYKG8lK697VRYq5RFQwVNosuqR96UL/ndtr2PUDcHwByGxtalMKyCsH2yjIe2WvmURUadqvFkJqD2WCxlRKT4bykjxTdXWlVBG1VFXxpb5spqgJEQZ8sYF0QgrHrrxHLQym0ouqxrIpgj9PBb+5uLpgJsknhIoTgrof/JC9hTUUjLKsiOFxXzafatpnRFCYZM3blVexKBfBgTu0WzViqo7HWFJXJuqlubUckcn/PFIWwKm1W7t9hxv2ZrJ+m2+7GEs/9bV4UwiqUwE6T4mPi4qu4PziPluO1LPOONdnUtBw5hnZphuh8EJpyF0NaFD2WiUnWSIkvOELYk9u1LFNYJpsbIZj3DRH25jbh0RSWyaZGsVg5+ImvEA/kNsrdFJbJpuf6O+eZuPp+Tm2azguTTc/2Q3djNXBT8uUweyyTTU/1jnbiwUhObZrCMtn0jH3wCtfe+7ec2sz7UFAI4QQ6gQ7grJTSm+cmmWwy2o99kdGZ3Ebt6NZjCSE6lhx3CyG6hBCPrfHSTuAC4AVcerXHxCRNIh5mcuQdtETuoi90EZYQogt4dtFxB4CUsg/wLhXdYlLnADillAN6tMfE5CaEQCoaEX/u0kd0EVZKHO5FT50g2QORer4LFnqxxX+cQoiTqeHfQAa9m4nJurFX1bDr4H2EfcGc2TRqjuUE5hYd1wFIKXuXniiEuJDq8VzALf9vYqIH7154loO/cSd1u7blxF7enRfm8M8kF9x+fzdaLHdrWUYJywuk9yF1AusuoC2EOAmcTB1GhRDvZvCyamBpgYOlzy0+Xu5xderYt+i5emAmC9ur2Vtqc6V/bRnYXsm++d4XHb/8jyu/9z9e/r3vz8D28kgpdfkDnF/0uAM4mXr8GNCxwWtfyPC8s2s9t/h4ucfA2cWPM7W/nO3V7C21udK/5nsv7Pe+0h+9vILdQGfq34XhXWru5JW5G+49n8Fzz6/x+PkljzdiezV7S22u9m+29s33vnZbNvrel0WklFnQCCEuSCk7t6J9870X53vPu/MiQ86mH6R6xRMkvY5XpZRncmzfSdKDWUtyqGB0pMhi2y6Sw+xa4Jkc2F6wn1qLrJVS9i2KlnECfQa24xbbyx0byHLv3UXyfXcBvVJK93IvLLhYweUiOAD3ojWuASnloyQXpM8ufX0O7HeSFHZaYLm0fQroS/05iU6sFjUjpUx/xun3DPB46qbecDuysL3ccS7td5Bcm+0Dule6bkEJK5MIjsW/EHr/UmYYQXKB5K/VaW5eFM+F7askeysXqbXBHNkFOLfoZc7UOV6gPce2lzvOmX0pZW/qfXexyrprQQlLriOCg5sXoHNp/7eklEeAR9Gx18jEduoX1Jn683qu7C7zMm8qasZJUuy5tK0b2dpPCbJvpWEgFJiwlmHZCA6Sv9i69RbrtO9OfbBgbKTILbZTv6AukuN9o2yv9Jl3AUdTYnoqddyFvsPxTGwvd5wz+6nvvgc4lfaCL0exOC9uIkcOi5VsGz1hXs32AJCXSJUln7mX3PywLWc759//EnvpueWqFHqPteEIjiK2ny/bW/E9626/0IV1jg89by4y+KXYRPbzZXsrvmfd7ReUsPIdwZFP+/myvRXfcy7sF0XkhYlJsVFQPZaJyWbBFJaJiQGYwtoiCCEeW632yDqu073a+o1JElNYWwevlHJACNEhhDi/tL5IKj7u6qLF72VJLUzXrnaOiSmsLUfK03WaZPjOUtz5XADfTJjC2rq4UykQaXSPvdzKFGVI02YjNfc5AZwneYN3kVycrAVci9IX9OQcybSHMyn7OQtR2gqYPVZhkI6onksN1Y6TvNEvAEcgmeSYmh89lnrsTM2L0vUZbzpey2BqrnR8tXPSNjf21rYmprAKgFT6gWvRSr83lfPTSbIXg1sT7JYmG2aTfOhNDQdvyWtLp+akHB5dBkSSb2pMYRUYqR4inWt1nGSF4I5lEuyWJhtmk3x4Dji9NK9osT0hhDMlWLOu/jowhVUApHqNdM+0OPhzlmTpuMUxbOkEu6XJhhklH6aucVoI0Z0aDqZFfBJwpf6/dlG7HtfzvW4VzFjBImFRgp2bpAgHSA4P4UNHx8Lx0rIFIlkjPyMnSEqctemeLF1AZ5HAM77WVsX0ChYJi+ZPi1k8hNMt+TA1BKxdtFg8J81S4OvCFNbWwZmaO2UkkFRvdYtQjao3stkwh4ImJgZgOi9MTAzAFJaJiQGYwjIxMQBTWCYmBmAKy8TEAExhmZgYwP8H1rHC2FI+tWoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 204.52x204.52 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = paper_settings.figure()\n",
    "ax = plt.gca()\n",
    "\n",
    "for i, shmf_ in enumerate(shmf_posteriors[:n_realizations_show]):\n",
    "    plt.plot(x, msubs**2 * shmf_,c=\"0.3\", lw=0.25)\n",
    "    \n",
    "fillpost2 = plt.fill_between(x, y_mult * shmf_posterior_cis[0, 1], msubs**2 * shmf_posterior_cis[1, 1], facecolor=paper_settings.COLOR_FIX, edgecolor=paper_settings.COLOR_FIX)\n",
    "fillpost1 = plt.fill_between(x, y_mult * shmf_posterior_cis[0, 0], msubs**2 * shmf_posterior_cis[1, 0], facecolor=paper_settings.COLOR_ALIGN, edgecolor=paper_settings.COLOR_ALIGN)\n",
    "lpost, = plt.plot(x, y_mult * shmf_posterior, alpha=1., c=paper_settings.COLOR_FULL, lw=1.0, label=\"Posterior\")\n",
    "ltrue, = plt.plot(x, y_mult * shmf_true, alpha=1., c=\"black\", lw=1.0, ls=\":\", label=\"True SHMF\")\n",
    "\n",
    "ax.add_patch(matplotlib.patches.Rectangle(xy=(1.93e7,1.147e-2), width=6.4e7, height=0.2e-2, facecolor=paper_settings.COLOR_FIX, edgecolor=paper_settings.COLOR_FIX))\n",
    "ax.add_patch(matplotlib.patches.Rectangle(xy=(1.93e7,1.195e-2), width=6.4e7, height=0.09e-2, facecolor=paper_settings.COLOR_ALIGN, edgecolor=paper_settings.COLOR_ALIGN))\n",
    "\n",
    "ax.legend(\n",
    "    [ltrue, lpost],\n",
    "    [\"True SHMF\", \"Posterior\"],\n",
    "    loc=2,\n",
    "    frameon=False\n",
    ")\n",
    "\n",
    "plt.xlabel(r\"$m_{200}$ [$\\mathrm{M}_{\\odot}$]\")\n",
    "ax.xaxis.set_label_coords(0.5, -0.11)\n",
    "plt.ylabel(r\"$\\frac {m_{200}^2}{M_{200}} \\frac{\\mathrm{d}n}{\\mathrm{d}m_{200}}$\")\n",
    "ax.yaxis.set_label_coords(-0.115, 0.5)\n",
    "ax.yaxis.set_tick_params(pad=1)\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.xlim(1.e7, 1.e12)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylim(0.001,0.02)\n",
    "#plt.ylim(0., None)\n",
    "    \n",
    "plt.savefig(\"../figures/shmf.pdf\", dpi=300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (lensing)",
   "language": "python",
   "name": "lensing"
  },
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}