Commit c845b0c6 authored by Moritz's avatar Moritz
Browse files

Final version

parent 8db1f704
......@@ -68,57 +68,67 @@
"source": [
"## Tutorial Setup\n",
"An example of how to use bilby to perform paramater estimation for\n",
"non-gravitational wave data consisting of a Gaussian with a mean and variance.\n"
"non-gravitational wave data consisting of a Gaussian with a mean and variance.\n",
"\n",
"### Setup on Google Colab\n",
"\n",
"Run the following line below:\n"
]
},
{
"cell_type": "code",
"metadata": {
"id": "4wXC8Zlg3HQA",
"id": "vrm7nHyRlxwl",
"colab_type": "code",
"outputId": "463299b2-4c9e-4457-b524-00bdf5782fcd",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 319
}
"colab": {}
},
"source": [
"!pip install bilby ptemcee nestle"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "s8B1aaCflsnS",
"colab_type": "text"
},
"source": [
"!pip install bilby emcee nestle\n",
"### Setup on CIT via Jupyter Lab\n",
"You need to change to the right environment. Click on the grey circle next to `Python 3` on the top right on your page and the select the `ligo-py37` environment.\n",
"\n",
"![alt text](https://i.imgur.com/ngHqVaB.png)\n",
"\n",
"![alt text](https://i.imgur.com/kwH6d5g.png)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "OKZsbcFcoG-K",
"colab_type": "text"
},
"source": [
"## Imports"
]
},
{
"cell_type": "code",
"metadata": {
"id": "4wXC8Zlg3HQA",
"colab_type": "code",
"colab": {}
},
"source": [
"import bilby\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"label = 'gaussian_example'\n",
"outdir = 'outdir'"
"%pylab inline"
],
"execution_count": 0,
"outputs": [
{
"output_type": "stream",
"text": [
"Requirement already satisfied: bilby in /usr/local/lib/python3.6/dist-packages (0.5.5)\n",
"Requirement already satisfied: emcee in /usr/local/lib/python3.6/dist-packages (2.2.1)\n",
"Requirement already satisfied: nestle in /usr/local/lib/python3.6/dist-packages (0.2.0)\n",
"Requirement already satisfied: matplotlib>=2.0 in /usr/local/lib/python3.6/dist-packages (from bilby) (3.0.3)\n",
"Requirement already satisfied: numpy>=1.9 in /usr/local/lib/python3.6/dist-packages (from bilby) (1.16.4)\n",
"Requirement already satisfied: dynesty>=0.9.7 in /usr/local/lib/python3.6/dist-packages (from bilby) (0.9.7)\n",
"Requirement already satisfied: scipy in /usr/local/lib/python3.6/dist-packages (from bilby) (1.3.1)\n",
"Requirement already satisfied: future in /usr/local/lib/python3.6/dist-packages (from bilby) (0.16.0)\n",
"Requirement already satisfied: dill in /usr/local/lib/python3.6/dist-packages (from bilby) (0.3.0)\n",
"Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (from bilby) (0.24.2)\n",
"Requirement already satisfied: corner in /usr/local/lib/python3.6/dist-packages (from bilby) (2.0.1)\n",
"Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=2.0->bilby) (2.5.3)\n",
"Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=2.0->bilby) (2.4.2)\n",
"Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=2.0->bilby) (0.10.0)\n",
"Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=2.0->bilby) (1.1.0)\n",
"Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from dynesty>=0.9.7->bilby) (1.12.0)\n",
"Requirement already satisfied: pytz>=2011k in /usr/local/lib/python3.6/dist-packages (from pandas->bilby) (2018.9)\n",
"Requirement already satisfied: setuptools in /usr/local/lib/python3.6/dist-packages (from kiwisolver>=1.0.1->matplotlib>=2.0->bilby) (41.2.0)\n"
],
"name": "stdout"
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -127,9 +137,9 @@
"colab_type": "text"
},
"source": [
"# Priors in Bilby\n",
"## Priors in Bilby\n",
"\n",
"Priors in Bilby are implemented as classes in Python. Using the object-oriented approach helps us doing a lot of things dynamically in the background without the user having to worry about them.\n",
"Priors in Bilby are implemented as classes in Python. Using an object-oriented approach helps us doing a lot of things dynamically in the background without the user having to worry about them.\n",
"\n",
"We implemented a range of different classes of priors into bilby. We are going to look at a power law prior of the form:\n",
"\n",
......@@ -153,16 +163,6 @@
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "GeP6w0L5PtPC",
"colab_type": "text"
},
"source": [
"Note that Bilby takes care of the normalization in the background for us."
]
},
{
"cell_type": "markdown",
"metadata": {
......@@ -170,7 +170,9 @@
"colab_type": "text"
},
"source": [
"If you are not familiar with object orientation in Python, take a moment to look at the next cell. *Objects* in Python have *attributes* that you can access and set with the '.' operator."
"Note that Bilby takes care of the normalization in the background for us.\n",
"\n",
"If you are not familiar with object orientation in Python, take a moment to look at the next cell. **Objects** in Python have **attributes** that you can access and set with the '.' operator."
]
},
{
......@@ -178,11 +180,7 @@
"metadata": {
"id": "wuJOd61zY-bM",
"colab_type": "code",
"outputId": "2eacdcbc-1e6e-4fa9-c4ab-eee43b24818c",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 67
}
"colab": {}
},
"source": [
"print('The minimum of the prior is ' + str(power_law_prior.minimum))\n",
......@@ -190,17 +188,7 @@
"print('The spectral index of the prior is ' + str(power_law_prior.alpha))"
],
"execution_count": 0,
"outputs": [
{
"output_type": "stream",
"text": [
"The minimum of the prior is 1\n",
"The maximum of the prior is 5\n",
"The spectral index of the prior is -3\n"
],
"name": "stdout"
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -209,7 +197,9 @@
"colab_type": "text"
},
"source": [
"The prior object also has a number of *methods*. Since a prior is a probability distribution, we implemented taking the probability at any given value, the associated cumulative probality density function, and also a method to sample from the distribution."
"The prior object also has a number of **methods**. Since a prior is a probability distribution, we have implemented the typical properties of such a distribution.\n",
"\n",
"You can evaluate the probability density at any given value, calculate the cumulative density function, draw random samples from the distribution, etc."
]
},
{
......@@ -219,7 +209,7 @@
"colab_type": "text"
},
"source": [
"### Probability distribution\n",
"### Probability density function\n",
"\n",
"We can get the prior probability at any given value by using the `prob` method. Let's evaluate the prior at some value."
]
......@@ -229,30 +219,13 @@
"metadata": {
"id": "K0cswskpQgl7",
"colab_type": "code",
"outputId": "0a5c2de4-b0e7-4a9f-f830-aa3a55ded175",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
}
"colab": {}
},
"source": [
"power_law_prior.prob(val=2.4)"
],
"execution_count": 0,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.15070408950617284"
]
},
"metadata": {
"tags": []
},
"execution_count": 70
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -269,11 +242,7 @@
"metadata": {
"id": "GUlkZP5-NX1m",
"colab_type": "code",
"outputId": "0e76cfe2-65f7-4589-bc24-c9beea6d272f",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 306
}
"colab": {}
},
"source": [
"grid = np.linspace(1.0, 5.0, 100)\n",
......@@ -286,31 +255,7 @@
"plt.clf()"
],
"execution_count": 0,
"outputs": [
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEQCAYAAACnaJNPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VHXe/vH3J4WEGnox9G4IvUkL\nFpooINW26lpABQVhXctjWXzUXXWVoiKKigXsSFdpFkKXIl16EQQBQXon398fifvwYykzyUzOzOR+\nXddcIScnM/fMurnnlPkcc84hIiLiqyivA4iISHhRcYiIiF9UHCIi4hcVh4iI+EXFISIiflFxiIiI\nX1QcIiLiFxWHiIj4RcUhIiJ+UXGIiIhfYrwOEAxFixZ15cuX9zqGiEjYWLx48e/OuWK+rBuRxVG+\nfHkWLVrkdQwRkbBhZlt9XVe7qkRExC8qDhER8YuKQ0RE/KLiEBERv6g4RETELyoOERHxi4rjLLPW\n72H7H0e9jiEiEtJUHBlOnUnj0THLaTs4lQ/nbSEtTddiFxE5HxVHhtjoKD67twn1yhXi6QmruHHE\nPDbuOex1LBGRkBNRxWFmHcxsxIEDBzL1+2UK5+HDuxrxcvfarNt1mGuHzuKNHzZw6kxagJOKiIQv\ncy7ydsk0aNDAZXXkyO5Dxxk4cRVfr/iNGpcV4MWutUhOTAhQQhGR0GJmi51zDXxZN6K2OAKpeP54\n3ri1Pm/+pR67Dp6g07A5vDRlDcdPnfE6moiIp1Qcl9AuuRTfDmhJ57qJvPHDRtq/OotFW/Z5HUtE\nxDMqDh8k5Inl5e61+fCuRpw4lUb3t+bxjwkrOXzitNfRRESynYrDDylVizGtfwp3NCnPh/O30nZw\nKjPX7fE6lohItlJx+ClvXAwDO9bgi3ubEB8bxR0jf2TA50vZf/Sk19FERLKFiiOTGpQvzFd9W/DA\nVZWZsHQHrQbN5OsVO72OJSISdCqOLIiPjebhttWY+EAzSibE0/ujJdw3ajG7Dx73OpqISNCoOAKg\nxmUJjO/djEfbVee7tbtpNWgmny/aRiR+RkZERMURIDHRUdx/ZSWm9GtB9ZIFeGTMcm4f+SPb9mlo\noohEFhVHgFUslo9Pe13Bs51qsGTrH7Qdksp7czZzRkMTRSRCqDiCICrKuK1JeaYNaEnD8oV5ZtJq\nerw1jw27D3kdTUQky1QcQZRYMDfv39mQQT1qs3HPYdoPnc3r363X0EQRCWsqjiAzM7rUK830/i1p\nXaMEL09bR8fX57Bie+Ym+IqIeE3FkU2K5Y9j2C31eOu2+uw9fIIb3pjDC99oaKKIhB8VRzZrW6Mk\n0we0pFu90rw5cyPXDp3Fgk17vY4lIuIzFYcHEnLH8mK3Wnx0T2NOp6Vx44j5PDV+JYeOn/I6mojI\nJak4PNSsclGmPpTCXc0qMHpB+tDE79fs9jqWiMhFqTg8lidXDE93SOLL+5uSNy6GO99fSP/PlrLv\niIYmikhoUnGEiHplCzG5b3P6XlOFSct20HrQTCYt26GxJSISclQcISQuJpoBrasy6cHmJBbKzYOf\n/ESvUYvZpaGJIhJCVBwh6PJSBRh7f1P+p311UtftodWgmXy28BdtfYhISFBxhKiY6Ch6pVRi6kMp\nJJUqwKNfruDWdxbwy14NTRQRb6k4Qlz5onn5pOcVPN85meXbD9B2SCrvztbQRBHxTsgXh5nlNbMP\nzOxtM7vV6zxeiIoybm1cjukDUmhSqQjPTl5N1+FzWbdLQxNFJPt5UhxmNtLMdpvZynOWtzOztWa2\nwcwey1jcBRjjnOsJdMz2sCGkVEJu3r2jAUNvqsPWvUe47tVZDJ2xnpOnNTRRRLKPV1sc7wPtzl5g\nZtHAMOBaIAm42cySgNLAtozVcvxgJzOjU51EZgxoybXJpRg8Yx0dX5/Nsm37vY4mIjmEJ8XhnEsF\n9p2zuBGwwTm3yTl3EvgU6ARsJ7084CJ5zayXmS0ys0V79uwJRuyQUiRfHK/eXJd3bm/A/qOn6PzG\nHJ7/ajXHTub4bhWRIAulYxyJ/N+WBaQXRiIwFuhqZsOBSRf6ZefcCOdcA+dcg2LFigU3aQhplVSC\naQNSuLFhGd6etZl2Q1OZt1FDE0UkeEKpOM7LOXfEOXenc+5+59xHXucJRQXiY/lXl1p83LMxADe/\nPZ/Hx67goIYmikgQhFJx/AqUOev70hnLxEdNKxVlSr8UeraowGcLf6HNoFS+/XmX17FEJMKEUnEs\nBKqYWQUzywXcBEz0OFPYyZ0rmieuS2Js72Yk5I7l7g8W0feTn9h7+ITX0UQkQnh1Ou4nwDygmplt\nN7O7nXOngQeAqcDPwOfOuVV+3m8HMxtx4IAuy1qnTEEmPdich1pV4ZuVO2k9OJUJS3/V2BIRyTKL\nxD8kDRo0cIsWLfI6RshY+9shHvlyOcu27eea6sV5rnMypRJyex1LREKImS12zjXwZd1Q2lUlQVKt\nZH7G3t+UJ6+7nDkbf6f1oFQ+WrCVNI0tEZFMUHHkENFRxj0tKjL1oRRqJibwxLiV3PLOfLb8fsTr\naCISZlQcOUy5Inn5uGdj/tWlJqt+PUjbIamMSN3I6TMaWyIivomo4tDBcd+YGTc3Ksv0AS1pUaUo\n//x6DV2Hz2XNbwe9jiYiYUAHx3M45xyTl+9k4MRVHDh2it5XVabPVZWIi4n2OpqIZCMdHBefmRkd\nal/G9AEtub5WKV79dj0dXpvNT7/84XU0EQlRKg4BoHDeXAy5qS4j/9qAQ8dP02X4XJ6dvJqjJ097\nHU1EQoyKQ/4/V1cvwbT+KdzauCzvzt5MuyGzmLvhd69jiUgIiaji0MHxwMgfH8tzN9Tk015XEB1l\n3PLOAh77cjkHjmlooojo4LhcwvFTZxg8Yx1vp26iaL44nrshmTY1SnodS0QCTAfHJWDiY6N5/NrL\nGd+nGYXz5qLXqMU88PESftfQRJEcS8UhPqlVOn1o4t9aV2Xaql20GjSTcT9t19BEkRxIxSE+i42O\n4sFrqvBV3+ZUKJqX/p8t4873F/Lr/mNeRxORbKTiEL9VKZGfMfc15R8dkliwaR9tBs1k1HwNTRTJ\nKSKqOHRWVfaJjjLubFaBaf1TqFu2EE+NX8lNI+azac9hr6OJSJDprCrJMuccXyzeznOTV3PidBr9\nW1flnuYViImOqPclIhFNZ1VJtjIzejQow4wBLbmyWjFe+GYNN7wxh9U7NDRRJBJdsjjMbLGZ9TGz\nQtkRSMJX8QLxvHVbA4bfWo/fDpyg4+uzeXnqWo6fOuN1NBEJIF+2OG4ELgMWmtmnZtbWzCzIuSSM\nXVuzFDMGpNCpTiKvf7+B616dxeKt+7yOJSIBcsnicM5tcM49AVQFPgZGAlvN7BkzKxzsgBKeCubJ\nxSs9avPBXY04fiqNbm/OY+DEVRw5oaGJIuHOp2McZlYLeAX4N/Al0B04CHwXvGgSCVpWLcbU/inc\nfkU53p+7hTaDU0ldt8frWCKSBT4d4wAGAwuBWs65vs65Bc65V4BNwQ4o4S9fXAzPdErmi/uaEBcb\nxe0jf+ThL5Zx4KiGJoqEo0uejmtmFZ1zm85ZVsE5tzmoyTLBzDoAHSpXrtxz/fr1XseR8zh+6gyv\nfruet1I3UThvLp7tVIN2yaW8jiWS4wX6dNwxPi7znHNuknOuV0JCgtdR5ALiY6N5pF11JvRpRrF8\ncdw3egn3j17M7kPHvY4mIj6KudAPzKw6UANIMLMuZ/2oABAf7GAS2ZITE5jwQDNGpG5i6Lfrmbtx\nL09dn0TXeonopD2R0HaxLY5qwPVAQaDDWbd6QM/gR5NIFxsdRZ+rKvN13xZUKZ6Ph79Yxh3vLWT7\nH0e9jiYiF+HLMY4mzrl52ZQnIDRyJPykpTlGzd/Ki1PWAPBou+rcdkU5oqK09SGSHfw5xnHB4jCz\nR5xzL5nZa8B/reSc65u1mMGj4ghf2/84yv+MW0nquj00KFeIF7vVolKxfF7HEol4gTo4/nPG10XA\n4vPcRAKudKE8fHBnQ17uXpv1uw9z7dBZDPt+A6fOpHkdTUQy+DUd18yigHzOuZCeXqctjsiw+9Bx\nBk5cxdcrfiOpVAFe6laL5ESdMScSDAE9HdfMPjazAmaWF1gJrDazv2c1pMilFM8fzxu31ufNv9Rj\nz+ETdBo2hxenrNHQRBGP+fI5jqSMLYwbgG+ACsBtQU2VSbqQU2Rql1yKGf1b0qVuIsN/2Ej7obNY\nuEVDE0W84ktxxJpZLOnFMdE5d4rzHCwPBfoAYORKyBPLv7vXZtTdjTh5Jo3ub87j6QkrOayhiSLZ\nzpfieAvYAuQFUs2sHOkDDkWyXYsqxZj6UAp/bVqeUfO30nZwKj+s3e11LJEcJVOXjjWzGOdcyL7V\n08HxnGHx1n08MmY5G/ccoUu9RJ6+PomCeXJ5HUskLPlzcPyCI0fOurM4oCtQ/pz1/zdT6UQCpH65\nwnzVtwWvf7eBN2duJHXdHv63UzLta2pookgw+bKragLQCTgNHDnrJuK5+NhoHm5bjQkPNKNkQjy9\nP1rCvaMWsfughiaKBIsvI0dWOueSsylPQGhXVc50+kwab8/azOAZ64iPieLJ65PoXr+0hiaK+CDQ\nY9XnmlnNLGYSCbqY6Cjuv7ISU/q1oHrJAjwyZjm3vfsj2/ZpaKJIIPlSHM2BxWa21syWm9kKM1se\n7GAimVWxWD4+7XUFz3aqwU+//EGbwam8N2czZ9JC8ixykbDjy66qcudb7pzbGpREAaBdVfKnX/cf\n44lxK/hh7R7qlS3Ii11rUaVEfq9jiYScgO6qyiiIMsDVGf8+6svviYSCxIK5ee+vDRl8Y202/36E\n616dzWvfrtfQRJEs8GVW1T+AR4HHMxbFAqODGUokkMyMznVLM31AS9rUKMEr09fR4bXZrNiu0TQi\nmeHLlkNnoCMZp+A653YAIbmtr1lVcjFF88Xx+i31GHFbffYdOUmnYbP51zc/a2iiiJ98KY6TLv1A\niAPImJIbkjSrSnzRpkZJpg9oSY8GZXhr5iauHTqLBZv2eh1LJGz4Uhyfm9lbQEEz6wnMAN4ObiyR\n4ErIHcsLXWvx0T2NOZ2Wxo0j5vPU+JUcOn7K62giIc+nWVVm1hpoAxgw1Tk3PdjBskJnVYk/jp48\nzSvT1jFyzmZKFYjn+c41uap6ca9jiWSrgFxzPJypOCQzlvzyB4+OWc763YfpXDeRp65PonBeDU2U\nnCEgp+Oa2SEzO3ihW+DiioSGemULMblvc/peU4VJy3bQetBMJi3bQSS+uRLJigsWh3Muv3OuADAU\neAxIBEqTfmrukOyJJ5K94mKiGdC6KpMebE5iodw8+MlP9Bq1mF0amijyH758cnyZc672pZaFEu2q\nkkA4fSaN9+Zs4eVpa8kVE8UT7S/nxoZlNDRRIlKghxweMbNbzSzazKLM7FY0Vl1ygJjoKHqmVGTq\nQykklSrAY2NXcOs7C9i6V//5S87mS3HcAvQAdmXcumcsE8kRyhfNyyc9r+D5zsks336AtkNSeWfW\nJg1NlBxLZ1WJ+GHngWM8OW4l367ZTe0yBXmpay2qlQzJQQoifgn0rioRyVAqITfv3NGAoTfVYdu+\no1z/2iyGzFjHydMamig5h4pDxE9mRqc6iUzvn0L7mqUYMmM9HV6bzdJt+72OJpItfJmOG50dQUTC\nTZF8cQy9qS7v3N6AA8dO0eWNOTz/1WqOndTQRIlsvmxxrDezf5tZUtDTiIShVkklmDYghZsaleXt\nWZtpNzSVeRs1NFEily/FURtYB7xjZvPNrJeZFQhyLpGwUiA+ln92rsnHPRsDcPPb83l87AoOamii\nRCC/zqoys5bAx0BBYAzwrHNuQ5Cy+c3MOgAdKleu3HP9+vVex5Ec6tjJMwyavpZ3Z2+meP54nrsh\nmVZJJbyOJXJRAT2rKuODfx3NbBzpo0ZeASoCk4Cvs5Q0wHQ9DgkFuXNF88R1SYzr3YyCeWK558NF\nPPjJT+w9fMLraCIB4dMxDqAT8G/nXF3n3CDn3C7n3BhgSnDjiYSv2mUKMvGB5vRvVZUpK3fSatBM\nJiz9VUMTJez5Uhy3O+fuds7N/XOBmTUDcM71DVoykQiQKyaKfq2q8FXfFpQrkpd+ny7lng8WsWP/\nMa+jiWSaL8Xx6nmWvRboICKRrGqJ/Hx5f1Oeuj6JuRv30mZwKh8t2EqaxpZIGIq50A/MrAnQFChm\nZgPO+lEBQJ/tEPFTdJRxd/MKtL68BI+PW84T41YycekOXuhaiwpF83odT8RnF9viyAXkI71c8p91\nOwh0C340kchUtkgeRt/dmBe71mT1zoO0G5LKiNSNnD6jsSUSHny5Hkc559zWbMoTEBpyKOFi18Hj\nPDl+JdNX76JW6QRe7FqLy0vpY1KS/QJyzXEzG+Kce8jMJgH/tZJzrmPWYgaPikPCiXOOr1f8xj8m\nrmT/0VP0vrISfa6uTFyM9ghL9vGnOC54jAMYlfH15axHEpELMTOuq1WKppWK8Ozk1bz63Qa+Wfkb\nL3arRb2yhbyOJ/JfdD0OkRDz/drdPDF2BTsPHufOphV4uG1V8uS62Hs8kawL1K6qFZxnF9WfnHO1\nMhcv+FQcEu4OHT/FS1PWMmr+VsoUzs0LXWrRrHJRr2NJBAtUcZS72C+G8gFzFYdEigWb9vLY2BVs\n/v0INzUsw+PtLychd6zXsSQCBaQ4wpmKQyLJ8VNnGDJjPW/P2kSRvLl47oZk2tQo6XUsiTABGXJo\nZrMzvh4ys4Pnfg1UWBG5uPjYaB67tjrjezejSL44eo1aTJ+Pl7DnkIYmijcuWBzOueYZX/M75wqc\n+zX7IooIQM3SCUx8oBl/a12V6at20XrwTMb9tF1DEyXb+XTNcTOrZ2Z9zexBM6sb7FAicn6x0VE8\neE0Vvu7XnIpF89L/s2Xc+f5CftXQRMlGvlyP42ngA6AIUBR438yeDHYwEbmwysXz88V9TRnYIYkf\nN++jzaCZjJq3RUMTJVv4MnJkLVDbOXc84/vcwFLnXLVsyJcpOjguOcm2fUf5n3ErmLX+dxqVL8wL\nXWtSsVg+r2NJmAnoFQCBHUD8Wd/HAb9mJpiIBF6Zwnn48K5GvNStFmt+O0i7obMY/oOGJkrwXGys\n+mukfwDwALDKzKZnfN8a+DF74omIL8yMHg3KcGXVYjw1YSUvTlnDVyt28FLX2iRdpnNZJLAu9gHA\nOy72i865D4KSKAC0q0pyum9W7OSpCavYf/Qk97WsxANXVyY+VkMT5cL0AUAVhwj7j57kua9+Zszi\n7VQqlpeXutWifrnCXseSEBXQYxxmVsXMxpjZajPb9Oct6zFFJJgK5snFy91r88FdjTh+Ko1ub85j\n4MRVHDlx2utoEuZ8OTj+HjAcOA1cBXwIjA5mKBEJnJZVizG1fwq3X1GOD+Ztoc3gVGat3+N1LAlj\nvhRHbufct6Tv1trqnBsIXBfcWP/HzCqa2btmNia7HlMk0uSLi+GZTsl8fm8T4mKjuO3dH3n4i2Uc\nOHrK62gShnwpjhNmFgWsN7MHzKwz6dcivyQzG2lmu81s5TnL25nZWjPbYGaPXew+nHObnHN3+/J4\nInJxDcsX5uu+Leh9ZSXG/fQrrQbPZMrKnV7HkjDjS3H0A/IAfYH6wG3ARc+4Osv7QLuzF5hZNDAM\nuBZIAm42syQzq2lmk8+5FffxcUTER/Gx0TzSrjoT+jSjeP447hu9hPtHL2b3oeNeR5MwccnLijnn\nFgJkbHX0dc4d8vXOnXOpZlb+nMWNgA3OuU0Z9/sp0Mk59y/gel/vW0SyJjkxgfF9mvH2rE0MmbGe\nuRv38uR1l9OtfmnMzOt4EsJ8OauqQcbVAJcDK8xsmZnVz8JjJgLbzvp+e8ayCz1+ETN7E6hrZo9f\nZL1eZrbIzBbt2aMDfyK+iI2OoveVlfmmXwuqFM/H38cs5/aRP7Jt31Gvo0kI82VX1Uigt3OuvHOu\nPNCH9DOtsoVzbq9z7j7nXKWMrZILrTfCOdfAOdegWLFi2RVPJCJUKpaPz+9twjMda7B46x+0HZLK\n+3M2a2iinJcvxXHGOTfrz2+cc7NJPzU3s34Fypz1fWk0+0rEc1FRxh1NyzOtfwoNyhdm4KTV9Hhr\nHht2H/Y6moSYi10BsJ6Z1QNmmtlbZnalmbU0szeAH7LwmAuBKmZWwcxyATcBE7NwfyISQKUL5eGD\nOxvySvfarN99mPZDZzHs+w2c0tBEyXCxWVXfX+T3nHPu6kveudknwJWkX8djF/AP59y7ZtYeGAJE\nAyOdc8/7G/wCj9cB6FC5cuWe69evD8RdiuRoew6dYODEVXy1YidJpQrwUrdaJCcmeB1LgkCzqjSr\nSiSgpqz8jacmrGTfkZP0SqlIv2uqaGhihAn0rKoEMxv05xlLZvaKmekth0gO0i65JDP6t6RrvUSG\n/7CR9kNnsXDLPq9jiUd8PavqENAj43aQbDyrSkRCQ0KeWF7qVpvRdzfm5Jk0ur85j6cnrOSwhibm\nOL5cOnapc67OpZaFEu2qEgmuIydO88q0dbw3dzOXJeTm+c7JXFlNgx7CWaAvHXvMzJqfdefNgGOZ\nDRdMZtbBzEYcOHDA6ygiES1vXAxPd0hizH1NyZ0rmr++t5ABny3ljyMnvY4m2cCXLY7apI9S//O4\nxh/AHc655UHOlmna4hDJPidOn2HYdxt444eNFMwTyzMdk2lfs6TGloSZgG1xZMynquacqw3UAmo5\n5+qGcmmISPaKi4lmQJtqTHygOaUSctPn4yXcO2oxuw9qaGKkumhxOOfSgEcy/n3QOXcwW1KJSNhJ\nuqwA43o35fFrqzNz3R6uGTSTzxduIxJP+c/pfDnGMcPMHjazMmZW+M9b0JOJSNiJiY7i3paVmPJQ\nCpeXKsAjXy7ntnc1NDHS+HKMY/N5FjvnXMXgRMo8fXJcJHSkpTk++vEXXvj6Z9Ic/L1tNe5oWp7o\nKB37CEX65LgOjouEjB37j/E/41bww9o91C1bkJe61qJKifxex5JzBPqT4/FmNsDMxprZl2b2kJnF\nZz2miOQElxXMzXt/bcjgG2uz5fcjXPfqbF77dj0nT2toYrjy5RjHh0AN4DXg9Yx/jwpmKBGJLGZG\n57qlmT6gJW2TS/LK9HV0fH02y7fv9zqaZIIvxzhWO+eSLrUslGhXlUhom756F0+OX8GeQyfo2aIi\n/VtX1dBEjwX6k+NLzOyKs+68MaC/yiKSaa2TSjCtf0t6NCjDW6mbaDcklfmb9nodS3zkS3HUB+aa\n2RYz2wLMAxqa2Qoz0wcBRSRTEnLH8kLXWnx8T2PSHNw0Yj5PjFvBoeOnvI4ml+DLrqpyF/u5c25r\nQBNlgU7HFQlPR0+eZtC0dYycs5kSBeJ5vnMyV1cv4XWsHEWn4+oYh0hY+umXP3j0y+Ws23WYG+pc\nxtMdalA4by6vY+UIgT7GISKSLeqWLcTkB1vQ75oqTF6+k1aDZjJx2Q6NLQkxKg4RCSm5YqLo37oq\nk/s2p0yh3PT95Cd6friY3w5oaGKoUHGISEiqXrIAY3s344n2lzN7wx5aD5rJJz/+oq2PEKDiEJGQ\nFR1l9EypyNSHUkhOTODxsSu49Z0FbN17xOtoOZqKQ0RCXrkiefm4Z2P+1aUmK7YfoO2QVN6ZtYkz\nadr68EJEFYcuHSsSucyMmxuVZdqAFJpVKspzX/1Ml+FzWfvbIa+j5Tg6HVdEwo5zjknLdzJw4ioO\nHT9F7ysr0+eqyuSKiaj3wtlKp+OKSEQzMzrWvozp/VNoX7MUQ79dT4fXZrN0m4YmZgcVh4iErSL5\n4hh6U13evaMBB46dossbc3hu8mqOnTzjdbSIpuIQkbB3zeUlmDYghZsaleWd2ZtpOySVuRt/9zpW\nxFJxiEhEKBAfyz871+STnlcQZXDL2wt4fOxyDmpoYsCpOEQkojSpVIRv+qVwb0pFPlu4jdaDZjJj\n9S6vY0UUFYeIRJzcuaJ5vP3ljO/TjEJ5cnHPh4vo+8lP7D18wutoEUHFISIRq1bpgkx8oDn9W1Xl\nm5XpQxMnLP1VY0uyKKKKQx8AFJFz5YqJol+rKnzVtwXli+al36dLufuDRezYf8zraGFLHwAUkRzj\nTJrj/blbeHnqWqKjjMeurc4tjcoSFWVeR/OcPgAoInIe0VHG3c0rMK1/CrXLJPDk+JXc/PZ8Nv+u\noYn+UHGISI5TpnAeRt/dmBe71mT1zoO0G5LKWzM3cvpMmtfRwoKKQ0RyJDPjxoZlmTGgJSlVi/Gv\nb9bQZfhcft550OtoIU/FISI5WokC8Yy4rT7DbqnHjv3H6PDabAZNW8uJ0xpbciEqDhHJ8cyM62qV\nYnr/lnSscxmvfreB61+dzZJf/vA6WkhScYiIZCiUNxeDetThvTsbcuTEaboOn8szk1Zx9ORpr6OF\nFBWHiMg5rqpWnKn9U/hL43K8N2cLbQanMnu9hib+ScUhInIe+eNjefaGZD6/twm5oqP4y7sLeHTM\ncg4c09BEFYeIyEU0qlCYr/u14P4rKzFmyXZaD5rJ1FW/eR3LUyoOEZFLiI+N5tF21RnfuxlF8sVx\n76jF9PloCXsO5cyhiRFVHJpVJSLBVLN0AhMfaMbDbaoyffUuWg+eydgl23Pc0ETNqhIRyYQNuw/x\nyJjlLPllPy2rFuOfXWqSWDC317EyTbOqRESCrHLx/HxxX1MGdkhi4ZZ9tBk0k1HztpCWFnlvxs+l\n4hARyaToKOOvzSow9aEU6pUrxFMTVnHjiHls2nPY62hBpeIQEcmiMoXz8OFdjfh3t1qs/e0Q7YbO\nYvgPkTs0UcUhIhIAZkb3BmWY8beWXF2tOC9OWUOnYXNY+Wvknayj4hARCaDi+eN587b6vHFrPXYd\nPEGnYXP499Q1HD8VOUMTVRwiIkHQvmYpZgxI4YY6iQz7fiPXvTqLxVv3eR0rIFQcIiJBUjBPLl7p\nUZsP7mrE8VNpdHtzHgMnruLIifAemqjiEBEJspZVizGtfwp3NCnPB/PShyamrtvjdaxMU3GIiGSD\nvHExDOxYgy/ubUJcbBS3j/yRv32+jP1HT3odzW8qDhGRbNSgfGG+7tuCPldVYvzSX2k9OJUpK3d6\nHcsvKg4RkWwWHxvN39tWZ0KPpToFAAAGl0lEQVSfZhTPH8d9o5dw/+jF7D503OtoPlFxiIh4JDkx\ngfF9mvFou+p8u2Y3rQel8sWibSE/NFHFISLiodjoKO6/shLf9GtB1RL5+PuY5dw+8ke27TvqdbQL\nUnGIiISASsXy8VmvJjzbqQZLtv5B2yGpvD9nc0gOTVRxiIiEiKgo47Ym5ZnaP4WG5QszcNJqur81\njw27D3kd7f8TUcWhCzmJSCQoXSgP79/ZkEE9arNxz2HaD53NsO83cCpEhibqQk4iIiFsz6ETDJy4\niq9W7OTyUgX4d7daJCcmBPxxdCEnEZEIUSx/HMNurcebf6nP74fThya+8I23QxNVHCIiYaBdcklm\n9G9J13qJvDlzI+2HzuLHzd4MTVRxiIiEiYQ8sbzUrTaj727MyTNp9HhrHk+NX8nhbB6aqOIQEQkz\nzasUZVr/FO5qVoHRC7bSZtBMfli7O9seX8UhIhKG8uSK4ekOSYy5ryl542L463sLGfDZUg4cOxX0\nx1ZxiIiEsfrlCjG5b3P6Xl2Zn7btJybKgv6YOh1XRCRCnDh9hriY6Ez9rk7HFRHJgTJbGv5ScYiI\niF9UHCIi4hcVh4iI+EXFISIiflFxiIiIX1QcIiLiFxWHiIj4JSI/AGhme4Ctmfz1osDvAYwj4SsB\n0FXB/BOpr1m4PK+s5CznnCvmy4oRWRxZYWaLfP30pEQ2MxvhnOvldY5wEqmvWbg8r+zKqV1VIhc2\nyesAYShSX7NweV7ZklNbHOfQFoeIyMVpi+O/jfA6gIhIKNMWh4iI+EVbHCIi4pcYrwOIRBIzywu8\nAZwEfnDOfeRxpJAXqa9ZpD4v0BaHRCAzK2Nm35vZajNbZWb9snBfI81st5mtPM/P2pnZWjPbYGaP\nZSzuAoxxzvUEOmb2cbObmcWb2Y9mtizjNXsmC/cVcq+ZmUWb2U9mNjkL9xFyz8srKo6LMLO8ZvaB\nmb1tZrd6nUd8dhr4m3MuCbgC6GNmSWevYGbFzSz/Ocsqn+e+3gfanbvQzKKBYcC1QBJwc8ZjlAa2\nZax2JovPIzudAK52ztUG6gDtzOyKs1cI89esH/Dz+X4Q5s/LEzmuOC70riGnvWOIZM65nc65JRn/\nPkT6H4zEc1ZrCYw3szgAM+sJvHae+0oF9p3nYRoBG5xzm5xzJ4FPgU7AdtL/YEAY/f/LpTuc8W1s\nxu3cM2fC8jUzs9LAdcA7F1glLJ+XlyLqyfjofc5515AT3zHkFGZWHqgLLDh7uXPuC2Aq8FnG1uRd\nQHc/7jqR//tvA9L/SCQCY4GuZjac8PnQGPCf3TlLgd3AdOdcpLxmQ4BHgLTz/TCMn5dnctzBcedc\nasYfk7P95x0DgJmd+45hKTmzZMOameUDvgQecs4dPPfnzrmXMv63Hg5UOusdd6Y5544Ad2b1frzg\nnDsD1DGzgsA4M0t2zq08Z52wes3M7Hpgt3NusZldeZEMYfW8vKY/huly3DuGSGdmsaSXxkfOubEX\nWKcFkAyMA/7h50P8CpQ56/vSGcvCnnNuP/A959+fH26vWTOgo5ltIX0X0tVmNvrclcLweXlKxXER\nzrkjzrk7nXP3R9KpdJHOzAx4F/jZOTfoAuvUJX1KQCfS3xUWMbPn/HiYhUAVM6tgZrmAm4CJWUvu\nHTMrlrGlgZnlBloDa85ZJ+xeM+fc48650s658hmP951z7i9nrxOOz8trKo50Oe4dQ4RrBtxG+rvL\npRm39ueskwfo4Zzb6JxLA27nPKP4zewTYB5Qzcy2m9ndAM6508ADpO8b/xn43Dm3KnhPKehKAd+b\n2XLS/xBOd86de+pqpL5mkfq8giZHjhzJOMYx2TmXnPF9DLAOuIb0wlgI3BLp/+OLiGRGjtviON+7\nhpz4jkFEJLNy5BaHiIhkXo7b4hARkaxRcYiIiF9UHCIi4hcVh4iI+EXFISIiflFxiIiIX1QcIiLi\nFxWHiIj4RcUhkg0yrnUxNOOyrCvMrKLXmUQyS8Uhkj0eBzY552oArwK9Pc4jkmk57kJOItnNzPIC\nnZ1z9TMWbSb9UqYiYUnFIRJ8rYAyGZdlBSgMzPAwj0iWaFeVSPDVAZ52ztVxztUBppF+OWKRsKTi\nEAm+QsBR+M+1X9qgSxFLGFNxiATfOuCKjH/3B75yzm32MI9Iluh6HCJBZmaFgG+AoqRfRKyXc+6Y\nt6lEMk/FISIiftGuKhER8YuKQ0RE/KLiEBERv6g4RETELyoOERHxi4pDRET8ouIQERG/qDhERMQv\n/w+gfmeBh5lSKgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"tags": []
}
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 0 Axes>"
]
},
"metadata": {
"tags": []
}
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -331,7 +276,9 @@
"source": [
"### Cumulative Distribution Function\n",
"\n",
"Bilby also implements the cumulative distribution function (`cdf`). We go through the same two steps to convince ourselves that this looks sensible"
"Let's continue with the cumulative distribution function (`cdf`). We go through the same two steps to convince ourselves that this looks sensible. The `cdf` of the prior $\\pi(\\theta)$ is defined as.\n",
"\n",
"$\\mathrm{CDF}(\\theta) = \\int_{\\theta_{\\mathrm{min}}}^{\\theta}\\pi(\\theta ') d\\theta '$"
]
},
{
......@@ -339,30 +286,13 @@
"metadata": {
"id": "RbcZxWCAZLWb",
"colab_type": "code",
"outputId": "07c25e37-fe67-405b-e475-90cd6df84c1a",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
}
"colab": {}
},
"source": [
"power_law_prior.cdf(val=2.4)"
],
"execution_count": 0,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([0.86082176])"
]
},
"metadata": {
"tags": []
},
"execution_count": 39
}
]
"outputs": []
},
{
"cell_type": "code",
......@@ -373,9 +303,9 @@
},
"source": [
"grid = np.linspace(1.0, 5.0, 100)\n",
"plt.plot(grid, p.cdf(grid))\n",
"plt.plot(grid, power_law_prior.cdf(grid))\n",
"plt.xlabel('$\\\\theta$')\n",
"plt.ylabel('cumulative distribution')\n",
"plt.ylabel('Cumulative Distribution')\n",
"plt.show()\n",
"plt.savefig('prior_cdf')\n",
"plt.clf()"
......@@ -400,30 +330,13 @@
"metadata": {
"id": "OVQFl765OWtA",
"colab_type": "code",
"outputId": "9ee3d0f9-fa5d-44b7-d807-c3960ef4ca21",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
}
"colab": {}
},
"source": [
"power_law_prior.sample(size=1)\n"
"power_law_prior.sample(size=1)"
],
"execution_count": 0,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([2.15582499])"
]
},
"metadata": {
"tags": []
},
"execution_count": 37
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -438,41 +351,19 @@
{
"cell_type": "code",
"metadata": {
"id": "NXzjo0V6cKkp",
"id": "lD4uDgirPzxH",
"colab_type": "code",
"colab": {}
},
"source": [
"# Binning\n",
"samples = power_law_prior.sample(10000)\n",
"plt.hist(samples, bins='fd', density=True, label='Binned samples')"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "2rgf1K21cLBY",
"colab_type": "code",
"colab": {}
},
"source": [
"samples = power_law_prior.sample(size=10000)\n",
"plt.hist(samples, bins='fd', density=True, label='Binned samples');\n",
"\n",
"# Analytic pdf\n",
"grid = np.linspace(1.0, 5.0, 100)\n",
"plt.plot(grid, p.prob(grid), label='Analytic pdf')"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "lD4uDgirPzxH",
"colab_type": "code",
"colab": {}
},
"source": [
"plt.plot(grid, power_law_prior.prob(grid), label='Analytic pdf')\n",
"\n",
"plt.loglog()\n",
"plt.legend()\n",
"plt.xlabel('$\\\\theta$')\n",
......@@ -491,9 +382,9 @@
"colab_type": "text"
},
"source": [
"#### Exercise:\n",
"#### Exercise (in your own time)\n",
"\n",
"Try the above steps with some of the other bilby priors. I pasted some of the class declarations below, just fill in the . \n",
"Try the above steps with some of the other bilby priors. You can use the templates below, just fill in the misssing parameter values.\n",
"\n",
"`bilby.core.prior.DeltaFunction(peak=)`\n",
"\n",
......@@ -517,7 +408,7 @@
"\n",
"We now want to look at how we can implement a likelihood in Bilby. We already provide a number of analytic likelihoods and an implementation for gravitational-wave transients, but it is instructive to see how we can implement one for ourselves.\n",
"\n",
"In order to implement a custom likelihood into bilby, you need to define a class that *inherits* from the base likelihood class in Bilby.\n",
"In order to implement a custom likelihood into bilby, you need to define a class that **inherits** from the base likelihood class in Bilby.\n",
"\n",
"Let's have a look at the implementation below."
]
......@@ -549,9 +440,7 @@
" sigma = self.parameters['sigma']\n",
" res = self.data - mu\n",
" return -0.5 * (np.sum((res / sigma)**2) +\n",
" self.N * np.log(2 * np.pi * sigma**2))\n",
"\n",
"\n"
" self.N * np.log(2 * np.pi * sigma**2))"
],
"execution_count": 0,
"outputs": []
......@@ -565,9 +454,9 @@
"source": [
"We need to implement an `__init__` method which is called upon instantiation. This takes in our data and defines the two parameters that are relevant for our problem, $\\mu$ and $\\sigma$.\n",
"\n",
"Additionally, we implemented the log-likelihood function for normal distributed data, which is defined as.\n",
"Additionally, we implemented the `log_likelihood` method for normal distributed data, which is defined as:\n",
"\n",
"$\\ln \\mathcal{L} = -\\frac{1}{2} \\sum_{i=1}^N \\left(\\frac{x_i-\\mu}{\\sigma}\\right)^2 - \\frac{N}{2} \\times \\ln (2\\pi\\sigma^2)$\n",
"$\\ln \\mathcal{L} = -\\frac{1}{2} \\left(\\sum_{i=1}^N \\left(\\frac{x_i-\\mu}{\\sigma}\\right)^2 - N \\ln (2\\pi\\sigma^2)\\right)$\n",
"\n",
"Here the $x_i$ are our data points sampled from a normal distribution, $\\mu$ and $\\sigma$ are the parameters that we want to estimate."
]
......@@ -619,6 +508,8 @@
},
"source": [
"plt.hist(data, bins='fd')\n",
"plt.xlabel('Value')\n",
"plt.ylabel('Count')\n",
"plt.show()\n",
"plt.savefig('normal_distributed_data')\n",
"plt.clf()"
......@@ -657,9 +548,11 @@
"colab_type": "text"
},
"source": [
"Let's continue by setting up the priors. Since we have multiple parameters, we need to bundle them together. In Bilby, we have defined a `PriorDict` as a container for this. The `PriorDict` can be used the same way as any normal python `dict`, but it also implement the `prob`, `cdf`, `sample` (and a lot more) methods of the `Prior` class. \n",
"Let's continue by setting up the priors. Since we have multiple parameters, we need to bundle them together. In Bilby, we have implemented a `PriorDict` as a container to do this. The `PriorDict` can be used the same way as any normal python `dict`, but it also implements the `prob`, `cdf`, `sample` (and a lot more) methods of the `Prior` class. \n",
"\n",
"We also add `latex_labels` to our priors. These will be used in our plots later.\n",
"\n",
"The syntax below shows how you set this up.\n"
"Here is how you can set this up:\n"
]
},
{
......@@ -671,8 +564,8 @@
},
"source": [
"priors = bilby.core.prior.PriorDict()\n",
"priors['mu'] = bilby.core.prior.Uniform(minimum=0, maximum=10)\n",
"priors['sigma'] = bilby.core.prior.Uniform(minimum=0, maximum=10)"
"priors['mu'] = bilby.core.prior.Uniform(minimum=0, maximum=10, latex_label='$\\mu$')\n",
"priors['sigma'] = bilby.core.prior.Uniform(minimum=0, maximum=10, latex_label='$\\sigma$')"
],
"execution_count": 0,
"outputs": []
......@@ -684,7 +577,7 @@
"colab_type": "text"
},
"source": [
"#### Excercise\n",
"#### Excercise (in your own time)\n",
"\n",
"Try out the `prob`, `sample`, `cdf` methods of the `PriorDict`.\n",
"\n",
......@@ -700,11 +593,10 @@
"source": [
"### Running inference\n",
"\n",
"We wanted to be able to use different samplers in Bilby without having to implement any of them ourselves. We managed to implement multiple samplers, which can all be accessed via the `run_sampler` function.\n",
"We can use different samplers in Bilby thanks to the interfaces we built. All the samplers are accessed via the `run_sampler` function.\n",
"\n",
"The `run_sampler` function takes the defined likelihood and prior objects, outdir and label arguments to define where we want our output to go.\n",
"Additionally, it takes all the keyword arguments that are present in the sampler. For example, for dynesty we define the number of live points `n_points`.\n",
"Dynesty also implements checkpointing. This allows us to stop and restart runs as well."
"The `run_sampler` function takes the defined `likelihood` and `priors` objects, `outdir` and `label` arguments to define where we want our output to go.\n",
"Additionally, it takes all the keyword arguments that are present in the sampler. For example, for dynesty we define the number of live points `npoints`."
]
},
{
......@@ -712,48 +604,15 @@
"metadata": {
"id": "k7GZPPTS40UR",
"colab_type": "code",
"outputId": "48d2c6e6-38ea-4c67-a1a3-bf1ba4c9fa56",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 423
}
"colab": {}
},
"source": [
"result = bilby.run_sampler(\n",
" likelihood=likelihood, priors=priors, sampler='emcee', nwalkers=1000, iterations=500, outdir=outdir, label=label, resume=False)"
" likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500, \n",
" walks=10, outdir='outdir', label='dynesty_run', resume=False)"
],
"execution_count": 0,
"outputs": [
{
"output_type": "stream",
"text": [
"10:33 bilby INFO : Running for label 'gaussian_example', output will be saved to 'outdir'\n",
"10:33 bilby INFO : Search parameters:\n",
"10:33 bilby INFO : mu = Uniform(minimum=0, maximum=10, name=None, latex_label=None, unit=None, boundary=None)\n",
"10:33 bilby INFO : sigma = Uniform(minimum=0, maximum=10, name=None, latex_label=None, unit=None, boundary=None)\n",
"10:33 bilby INFO : Single likelihood evaluation took 5.532e-04 s\n",
"10:33 bilby INFO : Using sampler Emcee with kwargs {'nwalkers': 1000, 'a': 2, 'args': [], 'kwargs': {}, 'postargs': None, 'pool': None, 'live_dangerously': False, 'runtime_sortingfn': None, 'lnprob0': None, 'rstate0': None, 'blobs0': None, 'iterations': 500, 'thin': 1, 'storechain': True, 'mh_proposal': None}\n",
"100%|██████████| 500/500 [09:19<00:00, 1.27s/it]\n",
"10:42 bilby INFO : Checkpointing sampler to file outdir/emcee_gaussian_example/sampler.pickle\n",
"/usr/local/lib/python3.6/dist-packages/emcee/autocorr.py:41: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" acf = np.fft.ifft(f * np.conjugate(f), axis=axis)[m].real\n",
"/usr/local/lib/python3.6/dist-packages/emcee/autocorr.py:43: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" return acf / acf[m]\n",
"/usr/local/lib/python3.6/dist-packages/emcee/autocorr.py:105: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" tau = 1 + 2 * np.sum(f[m], axis=axis)\n",
"10:42 bilby INFO : Max autocorr time = 20\n",
"10:42 bilby INFO : Discarding 60 steps for burn-in\n",
"10:42 bilby INFO : Sampling time: 0:09:37.814559\n",
"10:43 bilby INFO : Summary of results:\n",
"nsamples: 440000\n",
"log_noise_evidence: nan\n",
"log_evidence: nan +/- nan\n",
"log_bayes_factor: nan +/- nan\n",
"\n"
],
"name": "stderr"
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -772,31 +631,13 @@
"metadata": {
"id": "C5nsTcXEiGYt",
"colab_type": "code",
"outputId": "4686d4a4-9bcc-4b91-cf6c-86a94d6b0285",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 379
}
"colab": {}
},
"source": [
"result.plot_corner(truths=dict(mu=mu, sigma=sigma))\n"
],
"execution_count": 0,
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAFqCAYAAABBF3JVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl83Fd56P/PM/OdVaN9sWRbtmI7\nduI4iZOYrGRjDwkQCBAoW6AQCty2ub+0vYUbfqVtCv210PKCW0LCFpbcXLZLSBqgAUIWshk7G3Zi\nO/Eiy5Zl7ctoNMt3vuf3x8woY2lkz4xnpNHoeb9efkmaOV/pSNY8Ot9znvMcMcaglFKqdFwL3QGl\nlKo2GliVUqrENLAqpVSJaWBVSqkS08CqlFIlpoFVKaVKTAOrUkqVmAbWKicia0Xk9yKyR0SeEZEt\nhbYTkc+mH3dE5Nr5671Si5MG1iogIg+JSNccT38d+K4xZj3wN8BdIiIFtvs18CbgkZJ2XKkqpYG1\nACLyThH5qYh0i8iUiOwWkS+ISG2e11+ZHhVOiciwiHxfRJYV2y6Pr9cKXAjcCWCM+TUgwHmFtDPG\nPGmM2Vfo11dqqdLAWpi/ApLAZ0iN4G4DPgH8WkSO+7MUkUuBB4BR4DrgL4HLgN+KiK/QdnlaBRwx\nxiSyHjuQfryYdkqpPFgL3YFF5i3GmIGsjx8WkWHgu8AVwIPHufbvgG7gWmOMDSAiLwJ/AP4U+Foh\n7UTkh8CG9DXrgF+ISDz98fXGmN0n8X0qpU6CjlgLMCOoZvwh/XbFCS6/EPh1JlimP982YAh4e6Ht\njDHXG2M2G2M2A9uAN2c+zgqqB4EOEfFkff6u9OPZ8m2nlMqDBtaTd3n67YsnaJcE4jkejwGbimh3\nQuk/BFuBGwBE5PWk5k63F9NOKZUf0bKBxRORFcAzwHPGmNefoO1WwBhjLsh6bDWwH0gYY3yFtJvx\nuR8CbjDGHMjx3KmkpipagAhwozFma/q5bwL3GmPuPUG7zwEfBVqBCSAKXGiMOXSin5FSS5EG1iKJ\nSAh4CFgOnH+iICMi7wN+APwT8BWgCbgDuJhUwAwU0k4pVbk0sBZBRALAL4CzgcuNMX/M87p/JJVZ\n4AcM8EOgBthkjFlTaDulVGXSwFqg9ALPPaRSoF5vjHmywOtrgDVAvzHmaGbF3xjzwWLaKaUqjy5e\nFSCdq3oX8BpS6VAFBVUAY8ykMeaP6WD5JuA0UrueimqXR591S6tS80zzWAvzH8C7SM1/TorIhVnP\nHcrMs4rI5cBvgY8YY76Xfuwc4Crg6XT7VwN/DfyLMebxzCfJt10BMltVv5Fe7b9LRE4zs29Vjtfu\n16T+oHy7iK+v1NJjjNF/ef4jtRvJzPHvc1ntrkg/dkPWY2cAvye1o2qKVOD8cI6vkVe7PPubWcX3\nZD22B9hSZLuHSI3UF/z/Qv/pv0r+pyPWAhhjuvJs9xCpPNDsx3aSGn2e6Nq82uXpeFtVtxXRTimV\nBw2si5SI/AbYPMfTbzPGPDaf/VFKvUID6yJljHldHs2mt6pmjUa7OM6W1hO0U0rlQbMCqpjRLa1K\nLQjNY61yuqVVqfmngVUppUpMpwKUUqrENLAqpVSJVWRWQEtLi+nq6lrobuS0dygKwNpmf+4Ggy+l\n3raceuK2qiy2b98+aIxpXeh+qKWrIgNrV1cX27Yt0rz071ydevvh+xe2H0uYiHQvdB/U0qZTAUop\nVWIaWAv06fv38+n795e8rVKqelTkVEAle6J7vCxtlVLVQ0esSilVYhpYlVKqxDSwKqVUiekca4FW\n1s86fbokbZVS1UMDa4F+8L7TytJWKVU9dCpAKaVKTANrgW66Zy833bO35G2VUtVDpwIK9GxvuCxt\nlVLVQwPrSeq69Sm6R2IArG70cWDFAndIKbXgNLCepO6RGOZLlwEgNz8CGliVWvIqZo5VRG4UkW0i\nsm1gYGChu6OUUkWrmBGrMeYO4A6ALVu2VOx5MetbA2Vpq5SqHhUTWBeLO961vixtlVLVo2KmApRS\nqlpoYC3QjT/ew40/3lPytkqp6qFTAQXaMzA153OrG308tHcUgBtufYquJj3rSqmlSANrCR245QL4\nTgMA3TtiGliVWqJ0KkAppUpMA6tSSpWYTgUUaPPyUFnaKqWqhwbWAn352rVlaauUqh46FaCUUiWm\ngbVA779rF++/a1fJ2yqlqodOBRTo0FisLG2VUtVDR6xKKVViGliVUqrENLAqpVSJ6RxrgS5aXVeW\ntkqp6qGBtUBfuPqUsrRVSlUPnQpQSqkS08BaoOvufIHr7nyh5G2VUtWjYqYCRORG4EaAVatWLXBv\n5jYUSfBk93jqRFZSNViP11YptfRUTGBdLIcJAsRsM33ktVJKzaRTAUopVWIaWJVSqsQqZipgsXjt\nqQ08vHcs77ZKqaVHR6wF+uzrVxfUtpD2SqnqoIFVKaVKTANrga76xh8LaltIe6VUddA51gJNJZyy\ntFVKVQ8dsSqlVIlpYFVKqRLTwKqUUiWmc6wFuub0przzWK85vanMvVFKVSIdsRbor67sLKhtIe2V\nUtVBR6xlsrrRd0wFrAO3XLDAPVJKzRcdsRboiq89l1e7A7dcwOVr67l8bT3dI3oMtlJLiQZWpZQq\nMQ2sSilVYhpYlVKqxDSwKqVUiWlgLdC7z24tqG0h7ZVS1UEDa4E+ecnygtoW0l4pVR0qJrCKyI0i\nsk1Etg0MDCx0d+YUiScLaltIe6VUdaiYDQKL5ZTWN39zR1naKqWqR8WMWJVSqlpoYFVKqRLTwKqU\nUiWmgVUppUpMA2uBbtiyrKC2hbRXSlUHDawFuuH89oLaFtJeKVUdKibdarEYDCfK0lYpVT00sBbo\nnd97oSxtlVLVQ6cClFKqxDSwKqVUiWlgVUqpEtPAqpRSJaaBtUCfuKijoLaFtFdKVQcNrAW6/py2\ngtoW0l4pVR003apAPSPRsrRVSlUPDawF+sDdu8vSVilVPXQqQCmlSkxHrHnquvUpukdiAPgsWeDe\nKKUqmY5Y89Q9EsN86TIuX1vPhavrFro7SqkKVjGBdbEcJqiUUidSMVMBi+UwwZsvX1lw24f3jpWr\nO0qpClQxgXWxeMsZzWVpq5SqHhpYC7S7PwLAhrZg3m2VUkuLBtYCffwnLwHw0CfPzrutUmppqZjF\nK6WUqhYaWJVSqsQ0sCqlVIlpYFVKqRLTxasC3fK6VQW3fXjvH8vVHaVUBdLAWqDXrW8suO3qRh9y\n8yPT7x+45YKy9E0pVRk0sBbo2cNhADavCOXdNjuQZgKsUqp6aWAt0E0/3wvkl8daSFulVPXQxSul\nlCoxDaxKKVViGliVUqrENLAqpVSJ6eJVgT5/VVdZ2iqlqocG1gJdfEp9WdoqpaqHTgUU6PH9Yzy+\nP78TAQppq5SqHjpiLdBnfnkAyC83tZC2SqnqoSNWpZQqsYoJrHpKq1KqWlRMYDXG3GGM2WKM2dLa\n2rrQ3VFKqaJVTGBVSqlqoYtXBfry29aWpa1SqnpoYC1QPuUCi2mrlKoeOhVQoN/sGeE3e0ZK3lYp\nVT10xFqgW39zEMjvJIFC2iqlqoeOWJVSqsQ0sCqlVIlpYFVKqRLTOdbj6Lr1KbpHYkDqdFWllMqH\nBtbj6B6JYb502TGP3f7OU/O+vpC2SqnqoYG1QBvagmVpq5SqHjrHWqD7dg5x386hkrdVSlUPHbEW\n6EsPHwLgLWc0l7StUqp66IhVKaVKTAOrUkqVmAbWeba60Yfc/Ahdtz610F1RSpWJzrHOswO3XACA\n3PzIAvdEKVUuGlgL9P33bihLW6VU9dDAWqDORn9Z2iqlqkfFzLEulsMEf/hMPz98pr/kbZVS1aNi\nAutiOUzwtieOcNsTR0reVilVPSomsCqlVLXQwKqUUiWmgVUppUpMA6tSSpWYplsV6Ccf3FiWtkqp\n6qGBtUAtIU9Z2iqlqodOBRTozq193Lm1r+RtlVpsRGStiPxeRPaIyDMisqXQdiLy2fTjjohcO3+9\nLy8NrAW6c9tR7tx2tORtlapEIvKQiHTN8fTXge8aY9YDfwPcJSJSYLtfA28Cqqp4hgbWHLpufQq5\n+RE9QFAtGBH5lYgYEbk1z/adIvITERkTkXER+b8ismpGm5Ui8lUReUJEIunP31Vk/1qBC4E7AYwx\nvwYEOK+QdsaYJ40x+4rpQyXTwJpD5hDBTCUqpeaTiLwXOLuA9kHgQeA04EPAB4BTgd+JSE1W03XA\nu4ER4NGT7OYq4IgxJpH12IH048W0qyq6eLVAMnVZM+9rEFcAItII/Dvw34H/nedlHwPWABuMMS+n\nP8/zwEvAx4F/S7d7xBizLP38R4E3zNGHHwKZ0mzrgF+ISDz98fXGmN0FfVNLkAbWBZIdSLU2q8ry\n/wE7jDF3i0i+gfWtwJOZoApgjNkvIo8BbyMdWI0xTj6fzBhzfeZ9EXkIuMEYc2BGs4NAh4h4skaj\nXenHi2lXVTSwFugXH91UlrZKicirgQ9SwDRA2hnAz3M8vhN418n2KxdjzICIbAVuAL4hIq8nNXe6\nvZh21UbnWAsU9LoJet0lb6uWNhHxArcDXyziVruJ1LzpTMNA48n27Tj+DPiwiOwB/hV4nzHGAIjI\nN0XkrXm0+5yIHAIuAr4pIodEZGUZ+zwvdMRaoK891gvAJy9ZXtK2asn7GyAA/NNCdySbMeaK4zz3\nEnDxHM99NM92nwM+dzJ9rEQ6Yi3Qj54b4EfP5VeIu5C2aulKp0X9T+CzgE9EGkSkIf105uPj3fqM\nkHtkOtdIVpWZjljTum59iu6RGIDmr6r5tgbwAz/I8dxfpf+dAzw7x/U7Sc2zzrQReKEUHVSF0cCa\nlsldXQiaerXkPQtcmePx35EKtt8CXs7xfMa9wBdFZE0m2T6d+H8J8Lcl7WkWEVkLfBdoAyaBjxlj\nthXSTkQ+Syrvdh3wDmPMPeXq73zSwFoBNPVqaTPGjAIPzXw8veuz2xjzUNZjlwO/BT5ijPle+uFv\nAP8N+LmI3AIY4B+BHlILYtmf853pdzM7pK4SkQFgwBjzcIFdz2xVzaz23yUip2UWpvJs92vgLuDb\nBX7tirakA2sl3v7r6FWdgABustZHjDGTIvIaUhsLvp9u81vgJmNMeMb1P57x8dfSbx8Grsi7E69s\nVX1zug+/Tu//Pw/Ylm87Y8yT6Xb5fulFoWICq4jcCNyY/jAsIvO6u6MbkM/m314+ddynW/iIDObZ\ntmR9qlAtwOAJW5XW6nn+emVhjJkVbdKj11yPHwSuK+ZzFul4W1W3FdGuqlRMYDXG3AHcsdD9KAUR\n2WaMyVlCbanRn8XiIyK/ATbP8fTbjDGPzWd/FqOKCaxKqcpgjHldHs10S+txaB6rUqpgxpgBILNV\nleNtac2nXbXRwFoeVTGlUSL6s6heuqV1DjI7M0IppdTJ0BGrUkqVWEUuXrW0tJiurq6F7kZxBl9K\nvW05Ne9LHANJx+B2Ca7qSudbENu3bx80xrTmem4+frf2DkUhFmatd6Sg34O9Q1EA1jb7y9U1dZKO\n97uVrSIDa1dXF9u2LdIUt+9cnXr74fvzvsRxDJGEQ9DjwqWR9aSJSPdcz83b71YRvweq8h3vdytb\nRQbWpcblEkI+rduqVLXQOValSuzT9+/n032FF/T59P37+fT9+8vQIzXfdMSqVIk90T0OkRXFXaeq\ngo5Y1YJxHEM4lsRxNOVPVRcNrGrBRBIO41GbSCKvw0OVWjR0KkAtmKDHBVjpt0pVDw2sasFUazbE\nynofDE4Ud52qChpYlSqxH7zvNPjOzcVdp6qC3oMppVSJaWBVqsRuumcvNx15bXHX3bO3DD1S802n\nApQqsWd7wzDVVtx1qipoYFVqHs08wFIPi6xOGliVmkfdIzHMl1LbXfWo8+qlc6xKKVViOmJVqsTW\ntwZgeLi461RV0MCqVInd8a718J3/Xtx1qiroVIBSSpWYBlalSuzGH+/hxsNvLO66H+8pQ4/UfNOp\ngCqmR74sjD0DUxBrKu46VRV0xFrFtCyfUgtDR6xVbCHK8ukoWSkNrFUlHI6ztTfM+ctDhELeBSnL\nlxklg1WVJQGLMXO3lap+GliryNbeME8cSNUBfc36wuf4SkGLV8Pm5SEY7Z/+OHu31QmvU1VBA2sV\nOT/9wjx/AV+g1Vq8uhBfvnYtjPz2hO1WN/qQmx+Zrhnw5WvXzkPv1HzQwFpFQiHvgo1UK5mI3Ajc\nCLBq1aoF7s0rMgVYtGZA9Vm692tqyTDG3GGM2WKM2dLa2lr2r/f+u3bx/p5rirvurl1l6JGabzpi\nVdN0Rb80Do3FIFFb3HWqKuiIVU3TvFelSkMDqwJSo1XHMYS87iW9oq9UKSz5V5DjGMKxJI5jFror\nCyqScAjHk7hcotMASp2kJT/HqgntKZp/WjoXra6DscPFXaeqwpIPrBpQUjT/tHS+cPUp0F94CtUX\nrj6lDL1RC2FpRxNeCSiluP11HEPSGAzzN62gUxlKVZ4lH1hLKZJwsB3DfMY4XcmvPNfd+QLXHby2\nuOvufKEMPVLzbclPBZRS0OPCuIT5WvuJx5McGYvRHCh8KsO2HQYjNi1BC8vSv6+lNBRJgF34+VVD\nkUQZeqMWgr6iSsjlEtwiCOWJrLbt0Dcex7ZTo9Oe8Tj7hqOMxJIFT2UMRmx6x2MMRuxydFWpJU1H\nrItIJhgCtNd56azzAky/LURL0DrmbS46qlWqOBpYF5GZwdDrdbO2pbgjky3LRfsJAvLMQK6Uyo8G\n1kWi2NHjyez/z2dUq2Z77akNMNFd3HWqKugrZpEodvR4Mhsg8hnVqtk++/rVcOjx4q5TVUED6yJR\n7OhRN0AoNf/01bZIZEaPhS4ilXIDhMrPVd/4I1cdeFdx133jj2XokZpvOmJdBLRO6uIylXDAKfyl\nNaWbPKqGjlgXgfnaXTU+HuP+nQOMj2vBZaVOhgbWLMXuu49EEjx5YIxImXbOBD0u6vzlnyd9tGec\nR/aO82jPeFm/jlLVTqcCshS7gv58f4TtPWEALixDv+ar8tSlnXXHvF1oOgWiFisNrFmKXUE/qy14\nzNvFqq7Ox9VnlP+wveMZH4/xaM84l3bW4fJZi7JW7jWnN8Hk3uKuU1VBA2uWYkeGwaCHC7vqcz4X\njyfpGY/TWefF662c4FCufmUHxro6X8HXZ6YjAK46vYXFmCr2V1d2woGtebdf3eibPgJ7daMvdb1a\n1MoWWEVEjDFLvkhoz3iclwenAIrefnqyct1Sl6tf2YGxmNFv9nTEYi6+/WT3OFHb4cqbH2F14/H/\nwBy45YLp9zMBVi1uJQ+sIlJrjJnQoJpyMoVSSjXHmGvuuLPOSzyZJB5P/SvVqPVk52krYTriZF3x\ntecYnfxTnt34Y8yHLyvoOlUdSnqPJSLXAN8SkR+KyMr0Y0t61SFTKKWYwFWqNKtcWQVerxuv283B\n8Tg94/GT+vzZMoGxmGkApapFyQKriFwK/AvwVWAS+GeAfEeuInKjiGwTkW0DAwOl6taiVqo0q7l2\nX3XWeVnXEihqNK2UmlspR6xXAvcYYx4FPg+4ROQWEdksIie8LzTG3GGM2WKM2dLaurhvBUul0O2o\nhebhnsxoWik1t1IG1u3AeSLyGeBRoBtoBf4SOAd0WqDc9Pyr3PRuSM23kwqs6dHo6SKyzhhzP/At\noAZ40BjzaWPMXwI9wIcg/2mBpSQatXm+N0w0aud8LByO8+CeYcLhE8+DeiW139yrf76OMd93Q+8+\nu5Ur3IUXU3n32XqnVi2KDqwichVwH/Ap4B4ReZ8x5kek5ljDIrIx3XQPYIuITuTlsGc4yo4jk+wZ\njuZ8bGtvmCcOTLC1N3zCzzUcTTIWtRmOJsvZZXUCn7xkOddaTxZ1naoOBadbpW/na4A/Bz5ljLlX\nRC4Cvi8i9cDtgAP8dxFJAq8G3muMKd3ScwU52YT49U3+Y97OfGxNemHp/OWhWdfOTMcqZ8X/fL9P\n3YYKkXiSqPEUdZ2qDgW/AtO382ER2QbUiYjHGPOEiLwX+BFwGLgFuAI4HXi7MealEva5ovyue5QH\ndo9iOw5vO3NZXtdkH7Pi91ucNSNoznzsNetzb3WcmZ9azor/+Sb+n8yJBdXizd/cwWj8Bt7Ejwu+\nTlWHkxna9AGvBe4FEsaYP4jIB4EvA28zxvy0FB2sdBubgxxoi7GxOf86AXMdsxKN2rzQH6EpZLGy\nzn/CotZBjwvHceM4BscxZR0h5pv4rycWKFXEHKuIuAGMMV8DgsBtIlKfHrk+CjzPIq5BUGjK0ilt\nNVy/uZ1T2mry/hotQYv2Wi9+t2DbzvTX2zMc5cmDE/z+5VG2HZogfoJbQ5dLcLmEcDxJJOEUXfYw\nH/km/uuJBUrlGQBF5NXAKcaY7xtjkiLiNcbEjTHXi8jdpEapT4qIBVwG2Mf9hBUscyvrOKngcKK5\nwmJuv6NRm2cOjXNqq59o0odjDGCxvsmPbTuMxxMMRRL0jMdZ2xI47rxl9ghRb8OVqgzHDawi4iI1\nKr099aHUGGO+boyJi4jfGBM1xrxXRD4CLAfOBt5qjDlU/q6XRyZQOY4pS5CKRBJ85+mjDERSa3lv\n3FBDNGmmg+a5q+qOqTwFx5+3zC5UorfhSlWG4wZWY4xDaqHqu0ASuFhEAsaYfzfGRLPafRtARHzG\nmEV9rkcmUGXmLOcKUtkLUIUc8Pd8f4RIwqY16OXy1fVYlotQ+n8hO6B21nmn388OmMcbvS7malDV\n5IYty9h1cHtR1z28d6wMPVLzLd+IYAOdwHeB80Xk30TkCwAicrGInJtuVzUpVZmgmpm7nCmzADUY\nKWzW46y2IJevbeK9ZzazY3DqmONcMqX8etKFUTLvZ/clHEsu6O6qXBsa1LFuOL+dN1lPF3Wdqg75\nBtafA33GmN8C24BPAA3p514F9EL17ayaa4toPJ5kLJKgpcYqOGc0UxT75bEE23vCPN8fAVKLZo0+\nN2ua/dMj1uwCKa9MBzAv51/NJdeGBnWswXCCMVP4aRKD4fKcmabmX75RYQrYICIfA/4M+AKpkeuf\nAF+ptoCaMdecZc94nAOjMda1BFKr8rFk3gnxmVv5jU1+IvEkm9JFpiMJh6hj6KjzTRdFyS5A7XcL\nYUmNXLOnHuY6CaBcifq5NjRUEscxIK4FnWR+5/deYDT+Pt5WYB7rO7/3Qpl6pOZbXoHVGNMrIj3A\nZ0nttrpPRK4EXq7WoApzz1lmF68udCV+dDLBrsEIKxt8nLasBjzuVCbAlE3Q45pzJBpNGhxjiCbN\n9JwszH0SQLkyBHJtaKgkkYQD4tKJZrWgCvnL/g1Su6juS3/8sDGmpwx9qnjZ5faya6Y6jiFpDIa5\n/9Z0j0X49e4Rjo5N4RLB7xYGIzZ94TjR5Owk/0xuqivpMDiZwJrxd6yzzktHwM1LA+FjCrXM15HZ\nlSbocYFxdG+oWlB5v+qMMT3GmO2Z0n/pjIElLzshPpJwsB1D0jH0jcex7dk/ov0DEZ45PMGe/sj0\nCLQlaLG8zpdzvjaScBiNJHi+P8Lh0RiHwwls25n+/F6vm8FYku09kWMKtSzVRH2XS8A4+rupFlSx\ntQJUDkGPC+MS7KTJuWUVIGEMPsuFyy3TI0qXS+bcZBD0uAi7XbTVWgS9bjrrvLO2xGYKtOQq1DKX\nSDpYn9UWJBgsvGCIUmpuS+s+sYRmph1FozY7+iYRA163HDMCzd5qeumqOs7qqOHy1XV5jyiDHhct\nQR+b2mvwet2zRrihkJfXrG8iFDo2OGePbGd6vj9yTFZCKcXjSfYOTp1wS+6JjI/HuH/nAOPjiys1\n+hMXdfBWd+FlAz9xUUcZeqMWggbWImWnHTmO4dm+SZ7vDROxHURSI9DM6n122taOwRgT0dTbfIRj\nSY6m504zQTizjfZEGxOOl2t7VluQ8zpDnNVWeFrQiWTn4J6MTEWtR3vGS9Sz+XH9OW28xiq80PX1\n57SVoTdqISzaYikLLTvtKJJwqLcEYwx+9+wRaHba1tp6P/uaY6yt9xeUEpUZ9WYWo8Kx1GjQMoaX\nR2Osb/Lj9x/733m8+qyZfNpyOJkjv7Od7FHaC6VnJEq/U/jPtmdEc4OrhQbWImWnHTmOYcw2iAjR\npKHGbZjMym3N3jnV2Rzg2hovLUGL0ckEO/rDtIa8nNoczDkCTU0X+I6pXQBMj2KHJhPsG4riGMO6\nluD018zOby1ky20pZLImTlamotZi84G7dzOaeDfvLjCP9QN37y5Tj9R808BaoFw1AlwuYXN7DUGv\nm+BzLhzDrBzSzHRAnd+aXqgamrLZOxhjaNKmMeDNuYCVCcrhWJJQOr0LoDXoIZJw6GhJLWotr/Ec\n8zXnym9VSpWfBtYCzVWk2rJctAQ9OBhcyHSOakauXVyrG3w4q+sI+FzTt+u5Anck4RCOJ6nzW8fM\ns2I72CK0hbzU+d3EDdOfvyNoMeB30zFjGqDY4jG5lPJzKVVNKvLVkEianCvZlSCzIt/kdx9TVHow\nYrNveIqphEPMdrCTDtFk6rnMXKrfnc51TRe3tiwXG9prWNWYGlH2jcfpC8enF5wy86r+dGqW3y3T\nXzOzASBiO/SOxxiOJo/JMhi3DV7Lxbh9bHZc70SMp7pH6Z3If6V9ruyCYgvRKFXtKnLEmkgaBiN2\nzlvj+Tqsbq49+JkV+UyVqcytd0vQwmkK4HW7MBgs9yvbUzPTAGERHGMIi2AnnVR+asiDyyXTQarJ\n58ZyCQ3e1JzscCROIpka3UaTx9aIDfnc+N2p0fHMBarMx5k/AJmf11TMIZIwTMXy+8Nl2w4v9EeI\nJZJAYPr/xHFSC3Xttd6yHF6o1GJWkSNWj3t2oMiYq+JUqR0vZShzxlTI68Yyhr2DUziOYXmDD68l\nWC4XXoHesRijkwlcSYf+cByX4xBK56FabhdR22b/cBTbdqZHwpbHxfB4lO8/fZT4VJyYbegZmWL/\nSHR65Jo9nTBX6lXm8Xh6vjfz81oR8rAs5GFFaO5NAdkj1MGITdRO4nYLtuNMj1ojCYeI7VDn12mA\nmW6+fCXvth4t6jpVHSpyqOFxy5wv1sxcpSvp8HzvVM40o1LIlTKUOQL6jBY/vZEkp7UEGYglj1kk\nEgQRw8sjU7xwZJKNHUHiySTNVy7/AAAgAElEQVS/2DnK5etqueSUZiwrNVLdP5zk6ESUgXCcc5eH\naPC6eKE/yh/7J3nu0CRR2+H6M1sZsNzTUwuZxbDjHUedGdVbxnBkPE5z4JVgfDSaJOGk3obm2KiV\nPY+c+gMXwLYd+sMJLFcqYOtpBXN7yxnN1Lp38UplzfyvW93oQ25+BIDVjT4O3HJBGXqoyq0iA+vx\nZPbAP987xY4jkwBlqbaUK2Uok7DeMxKlPuil3uem1mvR1eSbDsCOMcRsByfpYLlcYBue6QmzdyDC\nipCHS1YmeXF0CpKp0eOuiE0yCT2hOGORBI93j7Om3stIxKbO42IgkqC1xnNMcMzuC6SOo86eusiM\nUgcnExwZj+NqCdBUe+wfiuPlmGbnv2ZGvrbtpBfoUs/paQWzdd36FN0jqT9I/6emjULrvu/uj/Bf\nN57JhvSmjUyAVYvPogusGfNRFzR7PjccjjM6EaMuYDhzmZuf7RilxkrSUV9DZ2MAy0qlRJmYzWDE\n5tDwFOtaAkSSDo4rSa3fTUOd8KuXRxicTOAgbO6ooTnkASN01nlJOAmGJ+Msbwhw3qpa4rZhKJog\nkHSozcoIgGOT5x3HsGtwisNjMRxjWBbyEvK6afK58QDxZJJ4PInX684rxzTXAYnFHJq41HSPxDBf\nuowrvvYctx25juu5u6DrP/6TlwB46JNnl6N7ah4t2vu4TIJ+qacBsvf1Ty86xZI8sG+UB/aNMzJp\n+NWeCDv6ozxxYAKfJzVnOjwR55F9wwxPJZiI2vzX7mEOj0cxSWFVQw2XnFLLS0fjTMXjGLE5MBSm\nKWCo8wi7jk4yPhnj0ECUkak4koD1jQFsJ0n3UIShyTgN3mP/q1wuIeB5papWyCM0BiwavW7C8SQu\nl+D3W3i9bg6OpI56KdUefqXU8S3aEWu5ZPbmLwul5hFHjeHIWIzTWv3E1jXQ0WQxFYlzoH+c09p8\ndIU89I9G+fbTh3h6f5i/m7JxCbg8Dn88HGbLaoexiGE4bHN4LIpbYCqR2hiw9eAUo1Mx7n5uiD19\no7TWB4nHHRIu6JlIsH8wzngsSY3fw2jcoT1rcL61N8yj+8YYjiS45rRmjLhYns5g8Gcdgph9659Z\nkLMdh3q/Z/q2vtpzUUXkRuBGgFWrVi1wb9RSoIE1i+MYwvFXclOjScPgZIKRyQQul2F1vZddgxF+\nt2uUnf1R9o4N8F97whweT7BvKMLQZJKbfA5eN+ztj7MrEeHBl0dZXeehqy2Ay4nz+N4xxqKGM5fX\ncFannwd2RYjGk/RFbNoahEvWNfKqjjqSboPlrqPBctNQl0ppchzD+JRNxHY4d1mQ4UiC9qCHw+EE\nPrdguV2zKmZlbv2zz9TyiUwvTgE5NzxUU/K/MeYO4A6ALVu2aNlLVXZLPrDOXPRxjKE+4CHocdE3\nMsULfRM0euFXu8fY3jvOjiMxVjYFeXHAxjE2T/RM4bcEr+XC77UAIWobfrsvQq3PTTQW50BwCkdi\nPL4vxmAYwoDIFLv6YjR4PVy9sQG3JOkbn2RzW4COeh+2CF0NQaJJg2UM3aMxGn1uDo7HCcds1rUE\nuXZTK4MRm6asXVeZoDozDzf7TC2/W/BG3MektM1Mb5trh5lS6sSWfGDN3B5nFn2CbhfRpKF3OMIX\nftfN3qEpmgIWhyfhse4YnnTges/mVs7sqGFNs5/WGg/pgxU484EgMdvw6bWdPH04zLOHw/SEE/x4\nR5ykAx7ADTR4k8SScfYOxagLuNnZF6d/IoHtjNBaH6A55CNqDIfCcRxjGAjbrGn209XgI2J7jlmx\nB5gZ+mbWCshOj5pZWDvzfjRqs2c4yvom/3ErY6nju+V1q+C/bi/uOlUVluSrJnuLaeb2uNlvEY4n\nmUo49I/F+MXuAZ7uHuGlEYekuJmyHT5+YTvvOLMF73FujQXBbwlv2NDIGzY0YoyhdzzO7U8c4dH9\n46xdFsSJRehsr+X5g5McDidoSyQ5OmpzdqdwenOQiJ0kNh7DJak535W1Xur9ntSoOp5k19FJ6paH\nCMfi/OC5Id5/djMNNX72DEdZ15Aa7dZ7XKxp8k/PseaTHpWpMQupFDYdqRbndesb4bHu4q5TVWFJ\nBtaZW0w76nzTZf1aAvDC0XF2HBljz5DDqA2ntXn529d0srqx8NQuEWFFvY+/f+NqHnx5lK882kvS\nuHh41zgrQxar2mp58sAke4dtav1Bzl/toc7jAreLRq8bO+Slzi3sG4/jOIatvWGeODBB1Lb51QtD\nPNUTJmE7XL2plR1HJonEkzTXpHZVddT7jtmOeyLrm/zYjkOdxzWdt6oK9+zhMEy1sTnQX/h1wOYV\nlXsKrsrPkgysmdtivztVPzVzexzypY6iHhyJ8/iBKKM2vGZdPX/7mk487pMLMiLCa09tZFN7DX97\n/372jsWxkwnq6iJMRmz8LqgPWQxNxjgw5qHRb9g7GKEt6KXXDXv6o9hJh86Qh4kOPw2Wm9oALAt6\naak1NHmFje1B1jf6sdPTEnPtipqr3oLfb7G8zk/veAz/HLUa1Ind9PO9cOS1PLSmsDzWm36+F9A8\n1mqwJANr9m1xKP0TsG2Hg8NTPPhyP7c+eJj+iMNbz2jmL169HHcJi70sq/Xylbev5ZZfdvP8kUkO\n7IziB85eBp2NASKJVG5rwrHZenCCjnofr15Vx6aOGuo8LrpHY7gc6JmM8ca1LZy70mZg0mbXUJQL\nuxoIZt3u2+nKVzjkPComu15shs6tKnXy9F4vbTBi88j+Cf75oT66xxzes7mVmy4tbVDNqPVZ/Os1\np3BJVx02EAMa6kJ4XLB/eJKn+8boG03QVmuxMmjhsoSukIfBqQROMsmLQ1M8fTCC7XJz9YZWLlzV\nwOb20KwR6mDE5oW+SV4YiBxT2i9TcjDXiDbf87SUUnPTYQmpW2PLGCZjMfYPJ7hsTT03Xtg+vdJ/\nPFOJJC8PRrFcQtDrYpNjsPIIxl7LxWdfv4qP/+QlJmJJfDLJUy9HGU5YRGNJOhttLu9qoKM5wIHh\nGM9NhTk0FmN9q4+LO0McGI2xvjE1h7p5ZW3Or9EStNjYXgPOsSNQ3eevVHktycCaSX6vs4TDkwkc\n2+aXu4f5wu+OsCzk5a+vWHncoBqOJbl35xD37hxiYDKBk5Vy/rvAJF5L+MWLw1zSVUd9YO4fsc9y\nccvrVvGJn75Mf6KWQ/3jeH1JjngcvJabyUSSgXCcejdMGQeMoTXoxee1aLGFg2GbulCSukDu0wcs\ny8XKhvLVUlBK5bYkA+uhsShPdo/THHTx0mCcI2OTfPOpAYanDP/r7WvnHM1FEw4/2zHInX84Sjxp\naApaXNRVx6qGVNm+uO3QcNhiIpbkXx86xL+54COvaufdZ7di5Ti9FWBdS4CPXtDO1584QpsHmgOG\nmG2o8bgxlsNj+yZY3+oliov2eh8d9X4st4uA2zVrIieT1O8YQ9ByEbEd2mo8els/zz5/VRf84j+K\nu05VhSUZWHvHojx/KExXq4dIXHjq0CR9kw4fOK9tumTbTAeGo/y3n73MZNxhRZ2Xqzc20ZFj1bxh\n2KLe7+Zjm5bxq90jfOOpPn783AD/dFVX6rY8h+vObOH+F4exk4akK0F7vZfaoJuRcYfmkJu1TQFs\nhBqfG0tcYGBFg3/6oMGMlqCFYwx2wuHZgUlGIjavWlWno9Z5dvEp9RA8XNx1qiosyaGMR1x4fW48\nboglkzx1cIr1rQE+eN6ynO33Dk3xZz99CcfADa9q46MXtucMqhkiwvJ6Hx85v533bG4h4Rj+/J69\n/KFnImd7yy3cfPkKjkzEeXnIMJmw6R+PMxqz2bQsxMb2Whp88NNn+ukbDTMRtzkyFsMrqfnhTLV/\ny0otSk0lDVE7iTEwGX2l9oGaH4/vH+PxyIrirts/VoYeqfm2JAPrmR0hrj2jhZX1Ae7dOcxkPMn/\nfG1nztv10Smbv/jZXtwi/On5ywreJLChLcinLumg3m9xyy8PMDiZyNnu7OUh3nxaI1MGnt4fY99g\njCSG01qDeL1u7tk5yoMvj/KTHSMMTMR5cSBCz3h81oF+QY+Ljjov56+o58yOEAGfq+zH2KhjfeaX\nB/hM32XFXffLA6XvkJp3SyqwZmqtZlbSV9X5OTSW4LwVtazKETCTjuGTP32JmO3w/vPapnc0FSrg\ncfO+c1tJGsNN9+zFmNwjyLee0UzSAWNBjc/FhkY/kaTD6GSCd29q4PxVtVx3RiNdDQHOWFZDZ7rq\n1fI63zGV/esCFk21Xta2BGgKevX4FKXm2ZJ6xY1OJni+N0zf8BTP94bZdmiC3vE4l62py9n+R88N\ncGQiwZtPzz2fWojmGg9vWN/I4fE49784nLPN+tYA7bUeljeFaK/3MJJweKZ3lG882cvgpMM7zm7D\nF/CzvMHHqW2pkezx8k4zaVXlPM1WKTXbkgqsQ1Gb/nCcp3rH+c2eIb63PbWX+6Ku2YF1Imbz7a1H\n6Wzwcs6K3ItOGZOxJInkiW+3X9UZYlnIw3881stUjttzEeHirjqe650kILC+OcCuvime7pngtwcG\neOLAGHE7yo6jk/SNTPFSf0RPA1CqAi2prIDV9T5cIsSjcbb1TDI0Zeio9dIUnH2Lf+/OYWzH8ObT\nm+bMaT00GuNbW48e89gVgSjttblHtyLC1Rub+PbWo9y3c4h3b26d3cdGP/GkIZI0xI1wwfIQI1MJ\ndvXGSJgEiaRhTYvDnn6D10rt9c8+wypXgepqKlqt1GKwpAKr4xgGJqPsPTqJSSbYfijMDa+anQmQ\nSDrc9XQ/y+u8OYNk3Hb4woOH5vw6fRNxHt47xuVrZ6fPdDb46Kj18J0/9PG2Tc34ZgS6TE7sM4fG\n8To269rrWFnvhaTB7XbznjMbORh2WFXnYSopr5wOmzmja8qmL5w6HTRTREWLVs+vL79tLfz8y8Vd\np6rCkgqse4aj/OrFUQ4NT7JrMBVorswR/HYPTDGVcLhmY+6519/tPXFKzEN7xwh4XJy/avZ20yvW\nNXD3MwPs6JvkvBnbUTOB7+nDNs8dHmdF7Th/dsVyLlrTyOntQWK4sdxxhmOGjW2BWYVV/G6h3m/R\n5H9lk4MWVplfm1eEoMCSgdPXqaqwpO4L1zf5uWxNDe0NFpOxOD4rVSt1pt39qcr7K+tnj+76xuM8\n2Z07H3WmX+4aYSSr+ElGZlT64tGpWc+11njwuIU4qeIsURtGJpJc1FlPe22AzjovfstNLJHMWVjF\nslwEPC7iWYkHWlhlfv1mzwi/Ca8u7ro9I2XokZpvS+qV5vdbBL0+bMeN5fGxqsGfs3rVnoEpAh4X\ntTm2tt7+ZF9BX/ORHAnffo+LOr+bF/sjs55zu2Q68HY1wDtOD3HuinqixhDyufF63WxsC9LZGMhZ\nWCXkc89ZuUrNj1t/c5Bb+y8u7rrfHCxDj9R8W3KvvrPagpy1PMTwpE1bKHde6nO9YWp97lmLVlOJ\nwlfgnz08mfPx5XVeXjia+7m2kIdTW/xsaLJoaKqhMeBmR98E/SPRVCUuTbFaElY3+pCbH0FufoSu\nW59a6O6oAiyZSbfMyrjfODQFLKK2Q2COUV0s6bAsNHsa4OhE7l1TJxK3nVnnZAW9bqKJ3BsFAh43\nkUSMvrBwZDTKg/vGCHjdJGzhyqAHv1t0lb9CdN36FN0jqfn61Y2zp5VOxoFbLph+X25+pKSfW5XX\nkgmsmZXxfYMRXjwaYSrhzFqRz4jEHWr9s6cBjoaLC6yDEZvlM1bjazwuYraDYwyuGSPjgMdFJO5w\n+goPIb+ftY0+ggEPF6xIFbPuDyd0lb9CdI/EMF8qfPuqqm5LJrBm5iNX1bgZnkzgdQuJZO4Roytd\n3GQmT5G31945SgYCzPWMS6B7aIpz1tbT2eznktWN+P3Hru7rKr9SlWnJvDIz85J947CmKUCd381E\nbPaKPUDQ4yYcmz2f2l5bXK2AphwBMJZ0sNySc/NBPOngdbt4cRBCL4ZZ29gwfUBg9veiKtPt7zwV\n/u8Xi7tOVYUlE1gzWoIWQZ8bnyWMTOUOrKsbfRwai816vHWOxa4TmXmrDxBLmDlHwNGEg9cSakOw\not7F+mYtpLKYbGgLgi93PYgTXqeqwpJ5tWYqW7lcwmlNPjpqvRwYjpLMccvf2eBjNJqcVYXK43Zx\nVkfuX/4jNR/A1L6FK6wdXGHtwNS+BVP7FobqP5Sz/YGRKI1z3MrvG47SOxJjfauwrM5PTzhB30Sc\naNTmwECEpw+NE43m/qOgZhORG0Vkm4hsGxgYKPvXu2/nEPeNF76L6r6dQ9y3c6gMPVLzraSBVUSu\nEZHPi8hXRaRFRIob4pVBZmdS71CELz/ei9vYRG3D3qHZSfob2oIkkobhHMn9b9jQmPPzt7tGcz7e\n5MweuYRjScaiSa45vWnWc2NRm8NjcRLAH3sNQY8Xv0t4YNcQz/eHefLQBE8emGDPcPQE37HKMMbc\nYYzZYozZ0to6uz5DqX3p4UN8afD84q57eO6t0mrxKNlUgIicB3wd+CRwLfBV4Ici8jtjTGFl0Qdf\ngu9cXaquYTD4HYPHgbFogvdFbK5PGv5HwGbF7320zqizusFO8ubAFM27LGr9s39EH2xJMjRHwepc\nbtj5p8d8PBG1uSZgs743QHDg2OyDSCLJ7wKpYF/rgTV7A0ReNKwxhjqfxZmWC9sxBHa5mLQdEsZQ\n63Xjdi2Zmw+lKl4pX43rgQeMMfcaYz4CPARcDVwuIpac4Czp7Nu1RKK4tKa5OAbiSUPSGOp9Fk1B\nD13NPnyWi7Ect9R+y4XllpzTAQAhryvnglR+fTGMRJN43JIzjzZqp8oJCpA0cHg0zpSdxHK58Fku\nvJaLWp9FwjEMR21GIzaTcT0hQKlKUsrFq6eAD4nIxcaYx40xt4uIC3g/8PCJRq3GmDuAOwC2bNli\n+PD9JeuYOAY7vcrv97lpdQkv9Uf42LeeY8+wzS+v3zRra+tDe0f5+wcO8q61LWxszz2vGrMdnuye\n4KG9Y5jat8z59e8841vT7z+6b4wHe8f48tvWkFg+u+jGbY/38rOjQySSBvcUrG2AG87u5EPnLWPE\n5U7tFrNcuONJxo5OMhy3aeuohRylD5esj+iuM7WwTmrEKiKbReR0EdlojNkHbAcuFZHTAIwxtwFx\n4K9PvqvFyxxXUhewprd6dtZ5WbcsQMIxHM6RAXDpKfXU+dz8avcIMTv3iNBnubh8bT1/94ZVefVj\nOJLg0X3jrKz3cnaOoGqM4ZnDYdY0+akBOmrg/ee18KHzljGZhIlYqvCK4xhcLqGl3s/FqxsIFhFU\nbduZPoRQKVVaRQdWEbkKuA/4FPAjEbkO+BawBnibiFyebroVmF1tZAHF40le6I9w9frU4tG2Q+FZ\nbdwu4R/etJpwLMndzwzMeU5VRtjTfNzH7aThrqcHEIEvvmVNzrZbD07w0mCU3QNTOMDKILTU1dIS\n8tHst5hKJInaNpGEM+sQwUKd7PVqbt9/7wa+v/I/i7vuvRvK0CM13wqeCkjPldYAfw58yhhzr4hc\nBPwA+AzwBeAG4HMichS4mNRca8XoHo2x/fAEEduwot7L1oMTvOPMllntzl4e4sYL27n9yT4e2Tee\ns3B1xhe3PAi8slCVffsP8MCeEYYjNv90VRfLchTPNsbwve39LKv1kIwlqPHBuvZGNrYEGY6mpjFs\nx5B0hKDHhT+9myt791U8nqRnPE5nnRevd/aW3Gy6e6t8Ohv94M2vtOSs61RVKPhVZVJDt7CIbAPq\nRMRjjHlCRN4L/Aj4f4wxnxORlcA5wN8YYyqqFlqtR3AL1PgFt0my7dAEu/ojnJYjQfv6za3sH47x\nwJ4R3C549SlzB9dcko7hob1j/KEnzLvObuHiHOdrAfx85xAvHI3QZsF7zq/Fwsvbz2xhTUtNVvDz\n0+R3E0k4BD2zd1/1jMd5eTCVUZB9XEsuunurfH74TD+Mncb19bsKvw64/py2cnRLzaOTmWPtA14L\nBACMMVuBDwCfEZG1xphDxpj7Ki2oAkwkDEGfi4DL4qq1bpqDHm79zcGcZQFFhL++YiWnNPn47Utj\nfHtrH/3po09OpD8c57bHj/D7/eNcdVojH7+wI2e7w2MxvvbYES5cXYt4QYyPM5Y30DdpODQeO6ZU\nYNzAeDQ1HTBTZ52XdS2B6eNa1MK47Ykj3DZ0TnHXPXGkDD1S863gwCoibgBjzNeAIHCbiNSnR66P\nAs8DFX106OoGH201Pmr9HqaSAaLxBL1jcb7y+96c7S238M13r+ejF7TTN5Hgtsf7+NkfB9k3FMWe\nsXPLGEMknuTnO4e4/Yk+wvEkn3vDav7mys6cRbUBbn/iCG6X8Fz3BOvr4KxOP6ctC9Je52U4kqBn\nPD69c8zvljkLWXu9bta2BE44DaCUKq+8pgJE5NXAKcaY7xtjkiLiNcbEjTHXi8jdwJeBJ0XEAi4H\nKnpFxOt18+quenrG40SpZ3vvJANhm1/tGmFlvY/3nTv7VswlwvvObeOajU3c/cwAP3l+gOePRLBc\nQq3PRWPQw3jU5jUmhmPghXiEN6xv5MYLO+bcugpw784hHt0/jh/oqocrz+zg0lOWsbYpgOOY6TnT\nzM4x/BahHCcbKKUqx3EDazoPNQjcnvpQaowxXzfGxEXEb4yJGmPeKyIfAZYDZwNvNcZU/L48y3LR\nWuOhxhXk3JUhwHBkUvjmU334LRfXnTV7MQug3m/xZxd18KEty3j2cJgnD06wrWeCiajNKU1+GiIW\n9X6L+6/ehHWccoEAz/aG+crvD3Ph6lqseJhXrWnm7Wc0c2pLcDotLDNXajkG0CNXlqrMaQKZ97OL\nYKvKc9zAaoxxSC1UfZfU7f3FIhIwxvy7MSaa1e7bACLiM8bMTgqtQJGEw9FwnM46P+/Z3I7LA/Fo\nnJjt8L8e68VrCW/ZmDuFClLFqC/qquOiGYtRqx5IreyeKKgOhBP8/X91s7zOB9EIF57aSF1NADuZ\n+1iVzJErauGV89SAuehpAotLvlkBNrAK+C7wURH5NyBmjPm0iFwMRI0xT5PaDLAoBD0uloW8xL1u\nHLcwNhln7wRc2BkknjT828OHGYnYfOC8tpw1U0/G2JTNX//nPmJJQ6vEuXR9G2euDDEwYdNUq+lP\nle5Epwb85IMb4e5/LPjz/uSDG0+mW6qC5Psq/jnwLmPMb0VkM3ArcGf6uVcBP4TpVKxFIbMb66Vw\nnH1DUyyv93D+qhCblvlYWevi4PAk3/nDUQYmE9x06Yo5F54KNRlP8j/u38+R8TgdluGSNSFet6Ge\nZaEAiWVCk19X9Be7lpAHrNlV0/K6TlWFfCfspoANIvIx4M9IbQLoFJE/Ab5ijCnsTOgK0uy3aK/1\nsbzGz4qGAG6X4YmDU3z8omZW17v5zxeG+fOf7eXgyMmX6dszEOEv7tnLy0NTxJMGxw2NwQB7BhIk\nHGgMeI47hxqN2jzfG9ZarBXuzq193Dmyqbjrti7al5LKkldgNcb0Aj3AZ0ltAPgH4EvAo4tplJpL\nQ42Hs5aHaAt5GJqI8b1t/eweCHNwxOGjFzQRFOgdj/Hxn7zE/S8On3Bray6JpMO3t/bxiZ++zNiU\nTcCBejd0tXg4fVkNm1eGaKn1MBa16Q8n5tzDv2c4yo4jk1qLtcLdue0od46cWdx1246WoUdqvhUy\nofcN4OfGmO3pjx9OL24taplFob7xONsPR6jzW1zUVc9fXNzBoYk4b9wwyXNHIjQ0h/jiQ4f4wfaj\nvGZdA5ecUk9Xo49gjpzRpDHsGYiwZ2CKP/SE2X5oYrq031DE5qxa2Li6jnee18alq5rw+yz8bsFv\nJbGTzpwnsK5v8uM4huU1nulCLEqpypN3YDXG9AA9IiImZdEH1WwtQYvLTqllXauPc9tD2CLsHBpk\nTUsN7Y0BRgbH6PHCynof/+fZAf73M6kjPoIe13Sequ0YfuxMknQMH//JywC01njwikN7LfRPgGVB\nR7PQUevhzLYgjsvFeNTGFfTQXufFth0styvnHn6/32Jda3B655VmCShVmYqtFVCVgn4PF6YT8/f0\nR2gNebhoTSMNPjd3PB7DbYVJxCa5dn2QrYciRBLgdjkcGYvjkC5OHXjlx2MB3liCllpYWW9x6akh\nukeiDMcd/nBokq5d47xjsx/L7ZqeWz3RHv5UO81nVaqSaW5PWqaMHqRW7o+Mx+mo83L2+nr8bqF3\nJEJfOMFFq0KsaQ7R1T7Mf/4xzKn1hjEbDg3BKTUQiIIxcFYThAIeLlwdpHvEodZy8eZNzYxNGHb3\nj2N5LN5zVjPB9IJVvrf1ms+qVOXTwJqWufUOueBQOM6qJi+nNgWn991vWFHP5aM2HfUe1rfW8Itd\nwwxPGsaDPv7ysuX8eu8Yzx0cRgCP5eJd57WxprWeV3fW8p+7Bnhg9wTGNvzJuct4ebSB9U1+/H4L\nxzHT1ap0zrQ6/OKjm+D7/29x1+VBd2FVPg2saZlb8Od7w+wZirKpo+aYYiZntYU4OhFj30CU+oCb\nv7y8DTHCuzY3snMwzoHBGEGPF6/lpsbr4vSORl67rolo0hBJgG0cJuIOjtvFWVmnB0zXAEBrAFSy\nQnZbBb1ucBWeEpdrITQX3YVV+TSwzrC+yX/M2wy/36LR52UyEaVvJE5tTZDPX7MOn+Ui4Qwz1h5k\ndauP9qc9xJKGM5f5iCQcmvxuNneGiNlwelto1tyozpkuDifabZXta4/1wtA5fLL5mYK+xtceS1VX\n++Qlywvun6osGlhn8PutY0aU2c5dHsKyXPgsw67+GM3LPbTV+7huUxsb2mvZ1BIg/IwQsx0e3jfJ\neavcuFw+Lutq4pSmEJ113lm3+zpnWn1+9NwAjJ1WcGD90XOpTBMNrIufBtYCBIMeLuyqZ3gizlRC\nWFbjJeRzEwZWNfrB46Yl6MEtwttPb8R2W7QELSzLdcKK/kqp6qGBtQgNNR7O8rqnb9+zb+ddbhet\nIS805j4yW6lS0oWsysRA1GsAAAYvSURBVKSBtQgzb9/1dr56ZRat5qs8YKF0IasyaWAto+HRKe5/\naZSrT22gqUGnAhaLmRkA+S5aLTQdvVYOqcSNVCIyAHQvdD9OQgswiL++EV+wjlhknOjYyEJ3aoGk\nfhbza7UxpjXzgYjcCNyY/nADsHue+rEQ33slqcbv/5jfrblUZGBd7ERkmzFmy0L3oxIs5Z/FUv7e\nYWl//5o8qZRSJaaBVSmlSkwDa3ncsdAdqCBL+WexlL93WMLfv86xKqVUiemIVSmlSkwDq5pXUuqz\nxNWisNT+3zWwqnkhImdAdZ9AoWYTkVpYev/vGljLRESuEZHPi8hXRaRFRJbsofEi8kbgLhE5daH7\nMt+W8u+BiFwDfEtEfigiK9OPLYmRqwbWMhCR84CvA08CNcBXgatFpH5BO7YAROStwN8BnzLGvLRU\nXliwtH8PRORS4F9Ifc+TwD/D0hm5amAtj/XAA8aYe40xHwEeAq4GLhcRa6kEl/T3+Q9A0hjzmIgs\nA24WkX8WkXMyt4lVbCn/HlwJ3GOMeRT4/P/fztmEWllFYfh5NTW6ZgZlIhQZDVIokzDKNG0gIUEa\nDUQSs5KiGghCgmJgRFCTbFD+DCLCSLgE/Ug1qSQkzd+oDMnUQaZlw7rF9aa9Dfa+cTng5Vz9zvng\n2+sZne98e7DWPuu8e+29197AKEnrJd0uaULNtnWcENbOsBeYImk2gO2twCFgGdBTyqid/bwHmCjp\nPeAd0sU/VwOrgNtqNK8blBwHB4E7JK0DdpHu/riW9LvPhGYvC4SwVkQeiadJmm77BCmw5kq6BcD2\nZmAAeK5OO7vB0EzU9l/AXSQR3W37ZdtPAb8Cy2sysWOUHAdDfL/Z9sfAm6QlkC9sr7W9CjgJPArN\nXhYIYa0ASQuBHcCzQK+kh0lBdROwSNK83HQf8Hc9VnYHSYuBtyXNHsxIsrhOB14YkqUcA/okNeYi\n25LjoMX3DyQ9YruXtMbaJ2l6bnoUOCdpbE2mdoU4eXUJZJHoAXqBLbY/knQ3acq7jjQVXAHMA84A\ns4EHbH9fj8WdJe/6fwKcAr4i/dH2tmYmklYCTwPLbf/QdUMrpuQ4GMb3bcCrwFbgddIS0HlgDrC0\nCb4PR1x0fQlkweiTdACYIGmM7T2SlpICbbXtDbnUZCawxvbPddrcYc4CjwGHgfXAEgBJ+22flzSO\nlL0tBlY0QVSh7Dhow/dTpFiYD0wDHrL9U20Gd4nIWCtA0jPALGCV7T/yd3OBjcAS28frtK+bSBpv\nu0/SeGADMBrozX+262yfGWxTr6XVU3IcDOP7a8Ai27/UaV+3iTXWCrC9CbgC2Czpqjxq7wK+I01/\niiGLqrJwvkjy/35JG4EvJV3ZNFEdXCcuMQ7a9L24mXFkrCMkH828Bjhi+/csIs7vtgP9pILwy4DV\nwLymjtatfdHybpTtf/PnPcAU4EHb33bf0uqRNAeYantbfh5reyB/bnQclOx7uxQ3klwKeefzFeAE\nMEbSk7ZP5dH5H9tLJT1OEpEZJCFpZEBdqC8G3w8R1VuBqcCCJmxYSBpFysy2pkf12N5ie0DS5bb7\nmxoHJfs+UiJjbRNJ80kX9y6zvU/S+8Abtj+TNNr2+Zb242yfrcPWTjOSvlA6vtlj+3RN5nYESWtI\n0/sZwDe2N16gXePioGTf2yWEtU0kTQMm294paTLpBM0+UvnM17bfUjobbtuHhi4RNI0R9MW5pkz9\nW5G0GriBVFK2knTg4azttUonrfqbGgcl+94usXnVJraP2N6ZH58ANtleDOwBFkq6EbgXOJ3bNzag\nRtAXZ+qxsCt8CPxm+3PgAKkud2J+N4tmx0HJvrdFZKwVIOlTUpnJ0bptqZtS+kLSFOAlYDewhlQQ\nfyfwLrC9yaJSsu/tEptXI6R1epOPLU4C/qzPqnoouS9sn5Z0EniedCXiDkn3AceaLiwl+94ukbFe\nJPkU0TJSOckS24drNqk2Su0LSdcDk2wfzM//l5g1nZJ9b4cQ1otE6Sb4BcBx2z/WbU+dlN4XJW/S\nlOz7cISwBkEQVExUBQRBEFRMCGsQBEHFhLAGQRBUTAhrEARBxYSwBkEQVEwIaxAEQcWEsAZBEFRM\nCGsQBEHF/AcvDh9ogc8/kQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 396x396 with 4 Axes>"
]
},
"metadata": {
"tags": []
},
"execution_count": 65
}
]
"outputs": []
},
{
"cell_type": "markdown",
......@@ -805,59 +646,13 @@
"colab_type": "text"
},
"source": [
"#### Excercise\n",
"#### Excercise (in your own time)\n",
"\n",
"Re-run the sampling and plotting step with `emcee` and `nestle` instead of `dynesty` as your sampler. By saving the results in different variables, you can plot multiple posteriors like this:\n",
"Re-run the sampling and plotting step with `ptemcee` and `nestle` instead of `dynesty` as your sampler. By saving the results in different variables, you can plot multiple posteriors like this:\n",
"\n",
"`bilby.core.result.plot_multiple([result_1, result_2, result_3])`\n",
"\n",
"Hint: You need to replace `npoints` with `nwalkers` for `emcee`. `emcee` also does not use the `walks` keyword argument. "
]
},
{
"cell_type": "code",
"metadata": {
"id": "qZMGhAr9izbw",
"colab_type": "code",
"outputId": "40fe3193-e301-4843-eec0-1bb6061c4fcb",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 742
}
},
"source": [
"#result\n",
"#result1\n",
"#result2\n",
"bilby.core.result.plot_multiple([result, result1, result2])"
],
"execution_count": 0,
"outputs": [
{
"output_type": "execute_result",
"data": {