diff --git a/.coverage b/.coverage index af77eba..a9c0853 100644 --- a/.coverage +++ b/.coverage @@ -1 +1 @@ -!coverage.py: This is a private format, don't read it directly!{"lines": {"/home/daniel/repositories/grbeams/grbeams/distributions.py": [1, 2, 4, 8, 9, 11, 14, 15, 16, 17, 26, 27, 28, 30, 39, 40, 42, 44, 47, 48, 49, 50, 59, 61, 62, 63, 65, 68, 69, 70, 71, 82, 83, 85, 86, 87, 89], "/home/daniel/repositories/grbeams/grbeams/scenarios.py": [1, 2, 3, 5, 6, 17, 18, 20, 34, 35, 36, 37, 39, 42, 44, 62, 67, 68, 70, 77, 81, 82], "/home/daniel/repositories/grbeams/grbeams/__init__.py": [1]}} \ No newline at end of file +!coverage.py: This is a private format, don't read it directly!{"lines": {"/scratch/aries/daniel/repositories/grbeams/grbeams/scenarios.py": [1, 2, 3, 5, 6, 17, 18, 20, 34, 35, 36, 37, 39, 42, 44, 62, 67, 68, 70, 77, 81, 82, 83], "/scratch/aries/daniel/repositories/grbeams/grbeams/distributions.py": [1, 2, 4, 8, 9, 11, 14, 15, 16, 17, 26, 27, 28, 30, 39, 40, 42, 44, 47, 48, 49, 50, 59, 61, 62, 63, 65, 68, 69, 70, 71, 82, 83, 85, 86, 87, 89], "/scratch/aries/daniel/repositories/grbeams/grbeams/beamingangle.py": [128, 1, 2, 3, 4, 5, 6, 8, 137, 138, 139, 12, 13, 142, 146, 132, 28, 32, 33, 36, 37, 39, 49, 99, 103, 115, 116, 118, 121, 124], "/scratch/aries/daniel/repositories/grbeams/grbeams/__init__.py": [1]}} \ No newline at end of file diff --git a/figures/overlay_ADE_jets.py b/figures/overlay_ADE_jets.py index 8eb573f..0a17168 100755 --- a/figures/overlay_ADE_jets.py +++ b/figures/overlay_ADE_jets.py @@ -20,6 +20,8 @@ from matplotlib import pyplot as pl import grbeams_utils +pl.style.use('/home/daniel/repositories/burst-style/burst.mplstyle') + resultsfiles = sys.argv[1].split('|') epoch=sys.argv[2] rate=sys.argv[3] diff --git a/final_paper/grb_beams_paper.tex b/final_paper/grb_beams_paper.tex index 5d3640b..e29e36f 100644 --- a/final_paper/grb_beams_paper.tex +++ b/final_paper/grb_beams_paper.tex @@ -201,7 +201,8 @@ \section{Short gamma-ray bursts and compact binary coalescences} \begin{figure} \centering -\includegraphics[width=\linewidth]{theta_dist_grbfrac.eps} +%\includegraphics[width=\linewidth]{theta_dist_grbfrac.eps} +\includegraphics[width=\linewidth]{color_relativenumber.pdf} \caption{\label{fig:thetapopulation} Expected relative numbers of observed GRBs and binary coalescences for different distributions on the GRB beaming angle. Lines in the figure correspond to jet angle diff --git a/grbeams/__init__.pyc b/grbeams/__init__.pyc index ce7275e..2d2bf88 100644 Binary files a/grbeams/__init__.pyc and b/grbeams/__init__.pyc differ diff --git a/grbeams/beamingangle.py b/grbeams/beamingangle.py new file mode 100644 index 0000000..14fcc0c --- /dev/null +++ b/grbeams/beamingangle.py @@ -0,0 +1,198 @@ +import numpy as np +import matplotlib.pyplot as plt +import emcee +import astropy.units as u +from grbeams.distributions import * +from grbeams.scenarios import * + + +from sklearn.neighbors.kde import KernelDensity +def kde_sklearn(x, x_grid, bandwidth=0.2, **kwargs): + """Kernel Density Estimation with Scikit-learn""" + kde_skl = KernelDensity(bandwidth=bandwidth, **kwargs) + kde_skl.fit(x[:, np.newaxis]) + # score_samples() returns the log-likelihood of the samples + log_pdf = kde_skl.score_samples(x_grid[:, np.newaxis]) + + N = np.trapz(np.exp(log_pdf), x_grid) + + return np.exp(log_pdf)/N + +class BeamingAnglePosterior(Distribution): + + def __init__(self, + observing_scenario, + efficiency_prior=DeltaDistribution(1.0),#'delta,1.0', + grb_rate=3 / u.gigaparsec**3 / u.year): + """ + The posterior probability distribution on the beaming angle, theta. + + Parameters + ---------- + observing_scenario : Scenario object + The observing scenario, e.g. O1. + efficiency_prior : str + The shape of the prior on epsilon, the efficiency at which + BNS mergers produce sGRBs. + grb_rate : + The rate of GRBs in the galaxy. + """ + + self.efficiency_prior = efficiency_prior + + # --- Prior Setup + + self.theta_range = np.arange(0.1, 90.0, 0.01) + self.efficiency_range = efficiency_prior.range + + # --- Astro configuration + self.grb_rate = grb_rate + self.scenario = observing_scenario + + def characterise_dist(self, x,y,alpha): + """ + Return the alpha-confidence limits about the median, and the median of the + PDF whose y-values are contained in y and x-values in x + """ + + # Peak + posmax = x[np.argmax(y)] + + # CDF: + cdf = np.cumsum(y) + cdf /= max(cdf) + + # Fine-grained values (100x finer than the input): + x_interp = np.arange(min(x), max(x), np.diff(x)[0]/100.0) + cdf_interp = np.interp(x_interp, x, cdf) + + # median + median_val = np.interp(0.5,cdf_interp,x_interp) + + # alpha-ocnfidence width about the median + q1 = (1-alpha)/2.0 + q2 = (1+alpha)/2.0 + low_bound = np.interp(q1,cdf_interp,x_interp) + upp_bound = np.interp(q2,cdf_interp,x_interp) + + # alpha-confidence *upper* limit + low_limit = np.interp(alpha,(1-cdf_interp)[::-1],x_interp[::-1]) + upp_limit = np.interp(alpha,cdf_interp,x_interp) + + #return [low_bound,upp_bound],median_val,[low_limit,upp_limit] + return [low_limit,upp_limit], median_val, posmax + + def get_theta_pdf_kde(self, bandwidth=1.0): + + self.theta_grid = self.theta_range #np.arange(self.theta_range[0], self.theta_range[1], 0.01) + self.theta_pdf_kde = kde_sklearn(x=self.theta_samples, + x_grid=self.theta_grid, + bandwidth=bandwidth, + algorithm='kd_tree') + self.theta_bounds, self.theta_median, self.theta_posmax = \ + self.characterise_dist(self.theta_grid, self.theta_pdf_kde, 0.95) + + def sample_theta_posterior(self, nburnin=100, nsamp=500, nwalkers=100): + """ + Use emcee ensemble sampler to draw samples from the ndim parameter space + comprised of (theta, efficiency, delta_theta, ...) etc + + The dimensionality is determined in __init__ based on the priors used + + The probability function used is the comp_theta_prob() method + """ + + theta0 = (max(self.theta_range)-min(self.theta_range)) * np.random.rand(nwalkers) + p0 = theta0.reshape((nwalkers, 1)) + + if self.efficiency_prior.ndim==2: + + efficiency0 = (max(self.efficiency_prior.range)-min(self.efficiency_prior.range)) * np.random.rand(nwalkers) + + p0 = np.transpose(np.array([theta0, efficiency0])) + + # Inititalize sampler + if self.efficiency_prior.name=='delta': + print "Delta Distro", self.efficiency_prior.ndim, + self.sampler = emcee.EnsembleSampler(nwalkers, self.efficiency_prior.ndim, + self.logpdf, args=[self.efficiency_prior.x]) + else: + self.sampler = emcee.EnsembleSampler(nwalkers, self.efficiency_prior.ndim, + self.logpdf_nparam) + + # Burn-in + + pos, prob, state = self.sampler.run_mcmc(p0, nburnin) + self.sampler.reset() + + # Draw samples + self.sampler.run_mcmc(pos, nsamp) + + # 1D array with samples for convenience + if self.efficiency_prior.ndim==1: + self.theta_samples = np.concatenate(self.sampler.flatchain) + else: + self.theta_samples = self.sampler.flatchain[:,0] + self.efficiency_samples = self.sampler.flatchain[:,1] + + + # Create bppu posterior instance for easy conf intervals and + # characterisation + #self.theta_pos = bppu.PosteriorOneDPDF('theta', self.theta_samples, + # injected_value=self.sim_theta) + + return self.sampler.flatchain + + def logpdf_nparam(self, x, fixed_args=None): + #print x + #print self.logpdf(theta=x[0], efficiency=x[1]) + return self.logpdf(theta=x[0], efficiency=x[1]) + + + def logpdf(self,theta,efficiency): + """ + Perform the rate->theta posterior transformation. + + Here's the procedure: + 1) Given an efficiency and theta angle, find the corresponding cbc rate + according to Rcbc = Rgrb / (1-cos(theta)) + 2) evaluate rate posterior at this value of the cbc rate + 3) The theta angle posterior is then just jacobian * rate + posterior[rate=rate(theta)] + """ + #print efficiency, self.comp_efficiency_prob(efficiency) + + if theta <= min(self.theta_range): return -np.inf + if theta > max(self.theta_range): return -np.inf + # Get BNS rate from theta, efficiency + bns_rate = self.scenario.cbc_rate_from_theta(self.grb_rate, theta, efficiency) + # Get value of rate posterior at this rate + bns_rate_pdf = self.scenario.comp_bns_rate_pdf(bns_rate) + # Compute jacobian + jacobian = self.compute_jacobian(efficiency,theta).value + theta_prob = bns_rate_pdf \ + + np.log(jacobian) \ + + np.log(self.comp_efficiency_prob(efficiency)) + if np.isnan(theta_prob ): + print bns_rate, bns_rate_pdf, jacobian, self.comp_efficiency_prob(efficiency) + return theta_prob + + def compute_jacobian(self,efficiency,theta): + """ + Compute the Jacboian for the transformation from rate to angle + """ + + denom=efficiency*(np.cos(theta * np.pi/180)-1) + output = abs(2.0*self.grb_rate * np.sin(theta * np.pi / 180.0) / + (denom*denom) ) + if np.isinf(output): + print efficiency, theta, denom, output + return 0 + else: + return output + + def comp_efficiency_prob(self,efficiency): + """ + Prior on the BNS->GRB efficiency + """ + return self.efficiency_prior.pdf(efficiency) diff --git a/grbeams/distributions.py b/grbeams/distributions.py index c20d659..a55e80e 100644 --- a/grbeams/distributions.py +++ b/grbeams/distributions.py @@ -47,6 +47,7 @@ class JeffreyDistribution(Distribution): """ name = "jeffrey" ndim = 2 + range = (0,1.0) def __init__(self): """ A Jeffrey distribution. @@ -83,7 +84,4 @@ def __init__(self, range=(0.0,1.0)): pass def pdf(self, e): - if (e>=min(self.range)) and (emin(self.range)) & (e<=max(self.range))) * 1./(max(self.range)-min(self.range)) diff --git a/grbeams/distributions.pyc b/grbeams/distributions.pyc index 49c9b0c..8d95866 100644 Binary files a/grbeams/distributions.pyc and b/grbeams/distributions.pyc differ diff --git a/grbeams/scenarios.py b/grbeams/scenarios.py index b9870ba..f2a3248 100644 --- a/grbeams/scenarios.py +++ b/grbeams/scenarios.py @@ -15,7 +15,7 @@ def __init__(self, rate_data, rate_pdf): The pdf at each computed rate value. """ self.rates = rate_data.to(u.megaparsec**(-3)*u.year**(-1)) - self.pdf_data = rate_pdf + self.pdf_data = rate_pdf#.to(u.megaparsec**(3)*u.year**(1)) def pdf(self, rate): """ @@ -31,10 +31,16 @@ def pdf(self, rate): probability : float The probability density of that rate. """ - rate = rate.to(u.megaparsec**(-3)*u.year**(-1)) + try: + rate = rate.to(u.megaparsec**(-3)*u.year**(-1)) + except: + pass if rate < 0 : return 0 if rate > max(self.rates): return 0 - return np.interp(rate.value, self.rates.value, self.pdf_data) + try: + return np.interp(rate.value, self.rates.value, self.pdf_data) + except: + return np.interp(rate, self.rates, self.pdf_data) class Scenario(): """ diff --git a/grbeams/scenarios.pyc b/grbeams/scenarios.pyc index 2e59fd2..49a9f51 100644 Binary files a/grbeams/scenarios.pyc and b/grbeams/scenarios.pyc differ diff --git a/grbeams_IDE.py b/grbeams_IDE.py index 5ed9aa9..94809ae 100755 --- a/grbeams_IDE.py +++ b/grbeams_IDE.py @@ -29,7 +29,7 @@ import matplotlib from matplotlib import pyplot as pl -import grbeams_utils +from .. import grbeams_utils from pylal import bayespputils as bppu diff --git a/notebooks/Scenario.ipynb b/notebooks/Scenario.ipynb index ecf9ef6..22c9a87 100644 --- a/notebooks/Scenario.ipynb +++ b/notebooks/Scenario.ipynb @@ -6,16 +6,7 @@ "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/usr/local/lib/python2.7/dist-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n", - " warnings.warn(self.msg_depr % (key, alt_key))\n" - ] - } - ], + "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", @@ -38,12 +29,13 @@ "cell_type": "code", "execution_count": 3, "metadata": { - "collapsed": true + "collapsed": false }, "outputs": [], "source": [ "from grbeams.distributions import *\n", - "from grbeams.scenarios import *" + "from grbeams.scenarios import *\n", + "from grbeams.beamingangle import *" ] }, { @@ -56,7 +48,7 @@ "source": [ "rateunits = u.megaparsec**-3*u.year**-1\n", "data = np.loadtxt('../O1/bns_rate_posterior.txt')\n", - "rate, pdf = data[:,0]*rateunits, data[:,1]*rateunits\n", + "rate, pdf = data[:,0]*rateunits, data[:,1]*rateunits**-1\n", "bns_prior = BNSDistribution(rate, pdf)" ] }, @@ -82,160 +74,6 @@ "grb_rate = 3.0 * rateunits" ] }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "class BeamingAnglePosterior(Distribution):\n", - "\n", - " def __init__(self, \n", - " observing_scenario, \n", - " efficiency_prior=DeltaDistribution(1.0),#'delta,1.0',\n", - " grb_rate=3 / u.gigaparsec**3 / u.year):\n", - " \"\"\"\n", - " The posterior probability distribution on the beaming angle, theta.\n", - " \n", - " Parameters\n", - " ----------\n", - " observing_scenario : Scenario object\n", - " The observing scenario, e.g. O1.\n", - " efficiency_prior : str\n", - " The shape of the prior on epsilon, the efficiency at which\n", - " BNS mergers produce sGRBs.\n", - " grb_rate : \n", - " The rate of GRBs in the galaxy.\n", - " \"\"\"\n", - "\n", - " self.efficiency_prior = efficiency_prior\n", - "\n", - " # --- Prior Setup\n", - " \n", - " self.theta_range = np.arange(0.0, 90.0, 0.01)\n", - " self.efficiency_range = efficiency_prior.range\n", - " \n", - " # --- Astro configuration\n", - " self.grb_rate = grb_rate\n", - " self.scenario = observing_scenario\n", - "\n", - " def get_theta_pdf_kde(self, bandwidth=1.0):\n", - "\n", - " self.theta_grid = self.theta_range #np.arange(self.theta_range[0], self.theta_range[1], 0.01)\n", - " self.theta_pdf_kde = kde_sklearn(x=self.theta_samples,\n", - " x_grid=self.theta_grid, \n", - " bandwidth=bandwidth, \n", - " algorithm='kd_tree') \n", - " self.theta_bounds, self.theta_median, self.theta_posmax = \\\n", - " characterise_dist(self.theta_grid, self.theta_pdf_kde, 0.95)\n", - "\n", - " def sample_theta_posterior(self, nburnin=100, nsamp=500, nwalkers=100):\n", - " \"\"\"\n", - " Use emcee ensemble sampler to draw samples from the ndim parameter space\n", - " comprised of (theta, efficiency, delta_theta, ...) etc\n", - "\n", - " The dimensionality is determined in __init__ based on the priors used\n", - "\n", - " The probability function used is the comp_theta_prob() method\n", - " \"\"\"\n", - "\n", - " theta0 = (max(self.theta_range)-min(self.theta_range)) * np.random.rand(nwalkers)\n", - " p0 = theta0.reshape((nwalkers, self.efficiency_prior.ndim))\n", - " \n", - " if self.efficiency_prior.ndim==2:\n", - " \n", - " efficiency0 = (max(self.efficiency_prior.range)-min(self.efficiency_prior.range)) *\\\n", - " np.random.rand(nwalkers)\n", - " \n", - " p0 = np.transpose(np.array([theta0, efficiency0]))\n", - " \n", - " # Inititalize sampler\n", - " if self.efficiency_prior.name=='delta':\n", - " self.sampler = emcee.EnsembleSampler(nwalkers, self.efficiency_prior.ndim,\n", - " self.comp_theta_prob, args=[self.efficiency_prior.x])\n", - " else:\n", - " self.sampler = emcee.EnsembleSampler(nwalkers, self.efficiency_prior.ndim, \n", - " self.comp_theta_prob_nparam)\n", - "\n", - " # Burn-in\n", - " pos, prob, state = self.sampler.run_mcmc(p0, nburnin)\n", - " self.sampler.reset()\n", - "\n", - " # Draw samples\n", - " self.sampler.run_mcmc(pos, nsamp)\n", - "\n", - " # 1D array with samples for convenience\n", - " if self.efficiency_prior.ndim==1:\n", - " self.theta_samples = np.concatenate(self.sampler.flatchain)\n", - " else:\n", - " self.theta_samples = self.sampler.flatchain[:,0]\n", - " self.efficiency_samples = self.sampler.flatchain[:,1]\n", - "\n", - "\n", - " # Create bppu posterior instance for easy conf intervals and\n", - " # characterisation\n", - " #self.theta_pos = bppu.PosteriorOneDPDF('theta', self.theta_samples,\n", - " # injected_value=self.sim_theta)\n", - " \n", - " return self.sampler.flatchain\n", - "\n", - " def comp_theta_prob_nparam(self, x, fixed_args=None):\n", - " return self.logpdf(theta=x[0], efficiency=x[1])\n", - "\n", - " \n", - " def logpdf(self,theta,efficiency):\n", - " \"\"\"\n", - " Perform the rate->theta posterior transformation.\n", - "\n", - " Here's the procedure:\n", - " 1) Given an efficiency and theta angle, find the corresponding cbc rate\n", - " according to Rcbc = Rgrb / (1-cos(theta))\n", - " 2) evaluate rate posterior at this value of the cbc rate\n", - " 3) The theta angle posterior is then just jacobian * rate\n", - " posterior[rate=rate(theta)]\n", - " \"\"\"\n", - " #print efficiency, self.comp_efficiency_prob(efficiency)\n", - " if (theta>=min(self.theta_range)) and \\\n", - " (thetaGRB efficiency\n", - " \"\"\"\n", - " return self.efficiency_prior.pdf(efficiency)" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -245,32 +83,29 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "text/latex": [ - "$-74.641016 \\; \\mathrm{\\frac{1}{yr\\,Mpc^{3}}}$" - ], "text/plain": [ - "" + "0" ] }, - "execution_count": 16, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "o1scenario.cbc_rate_from_theta(theta = 30, grb_rate=grb_rate, efficiency=0.3)" + "o1scenario.cbc_rate_from_theta(theta = 30, grb_rate=grb_rate, efficiency=-0.3)" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -281,114 +116,66 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { - "data": { - "text/plain": [ - "1857688383.3334398" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "np.exp(thetapos.logpdf(0.3, 0.3))" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": { - "collapsed": false - }, - "outputs": [ + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/astropy/units/quantity.py:821: RuntimeWarning: divide by zero encountered in true_divide\n", + " return super(Quantity, self).__truediv__(other)\n", + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/emcee/ensemble.py:335: RuntimeWarning: invalid value encountered in subtract\n", + " lnpdiff = (self.dim - 1.) * np.log(zz) + newlnprob - lnprob0\n", + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/emcee/ensemble.py:336: RuntimeWarning: invalid value encountered in greater\n", + " accept = (lnpdiff > np.log(self._random.rand(len(lnpdiff))))\n" + ] + }, { "data": { "text/plain": [ - "7.5267680522895066" + "array([[ 1.17961586e-01],\n", + " [ 1.92689986e-01],\n", + " [ 1.72285736e-01],\n", + " ..., \n", + " [ 4.74927905e-07],\n", + " [ 4.74927905e-07],\n", + " [ 4.74927905e-07]])" ] }, - "execution_count": 14, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "thetapos.logpdf(theta = 30, efficiency = 0.3)" + "thetapos.sample_theta_posterior()" ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "emcee: Exception while calling your likelihood function:\n", - " params: [ 11.91949427]\n", - " args: [0.3]\n", - " kwargs: {}\n", - " exception:\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Traceback (most recent call last):\n", - " File \"/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/emcee/ensemble.py\", line 505, in __call__\n", - " return self.f(x, *self.args, **self.kwargs)\n", - " File \"\", line 116, in comp_theta_prob\n", - " bns_rate_pdf = self.scenario.comp_bns_rate_pdf(bns_rate)\n", - " File \"/home/daniel/repositories/grbeams/grbeams/scenarios.py\", line 68, in comp_bns_rate_pdf\n", - " return self.bns_prior.pdf(rate)\n", - " File \"/home/daniel/repositories/grbeams/grbeams/scenarios.py\", line 36, in pdf\n", - " if rate > max(self.rates): return 0\n", - " File \"/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/astropy/units/quantity.py\", line 409, in __array_prepare__\n", - " result = self._new_view(obj, result_unit)\n", - " File \"/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/astropy/units/quantity.py\", line 589, in _new_view\n", - " view.__array_finalize__(self)\n", - " File \"/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/astropy/units/quantity.py\", line 283, in __array_finalize__\n", - " self._unit = getattr(obj, '_unit', None)\n", - "KeyboardInterrupt\n" - ] - }, - { - "ename": "KeyboardInterrupt", - "evalue": "", + "ename": "NameError", + "evalue": "global name 'characterise_dist' is not defined", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mthetapos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msample_theta_posterior\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[1;32m\u001b[0m in \u001b[0;36msample_theta_posterior\u001b[1;34m(self, nburnin, nsamp, nwalkers)\u001b[0m\n\u001b[0;32m 71\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 72\u001b[0m \u001b[1;31m# Burn-in\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 73\u001b[1;33m \u001b[0mpos\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstate\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun_mcmc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnburnin\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 74\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreset\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 75\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/emcee/sampler.pyc\u001b[0m in \u001b[0;36mrun_mcmc\u001b[1;34m(self, pos0, N, rstate0, lnprob0, **kwargs)\u001b[0m\n\u001b[0;32m 155\u001b[0m \"\"\"\n\u001b[0;32m 156\u001b[0m for results in self.sample(pos0, lnprob0, rstate0, iterations=N,\n\u001b[1;32m--> 157\u001b[1;33m **kwargs):\n\u001b[0m\u001b[0;32m 158\u001b[0m \u001b[1;32mpass\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 159\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mresults\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/emcee/ensemble.pyc\u001b[0m in \u001b[0;36msample\u001b[1;34m(self, p0, lnprob0, rstate0, blobs0, iterations, thin, storechain, mh_proposal)\u001b[0m\n\u001b[0;32m 257\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mS0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mS1\u001b[0m \u001b[1;32min\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfirst\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msecond\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0msecond\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfirst\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 258\u001b[0m q, newlnp, acc, blob = self._propose_stretch(p[S0], p[S1],\n\u001b[1;32m--> 259\u001b[1;33m lnprob[S0])\n\u001b[0m\u001b[0;32m 260\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0many\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0macc\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 261\u001b[0m \u001b[1;31m# Update the positions, log probabilities and\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/emcee/ensemble.pyc\u001b[0m in \u001b[0;36m_propose_stretch\u001b[1;34m(self, p0, p1, lnprob0)\u001b[0m\n\u001b[0;32m 330\u001b[0m \u001b[1;31m# Calculate the proposed positions and the log-probability there.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 331\u001b[0m \u001b[0mq\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mrint\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mzz\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnewaxis\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mrint\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0ms\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 332\u001b[1;33m \u001b[0mnewlnprob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mblob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_lnprob\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 333\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 334\u001b[0m \u001b[1;31m# Decide whether or not the proposals should be accepted.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/emcee/ensemble.pyc\u001b[0m in \u001b[0;36m_get_lnprob\u001b[1;34m(self, pos)\u001b[0m\n\u001b[0;32m 380\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 381\u001b[0m \u001b[1;31m# Run the log-probability calculations (optionally in parallel).\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 382\u001b[1;33m \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlnprobfn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 383\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 384\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/emcee/ensemble.pyc\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, x)\u001b[0m\n\u001b[0;32m 503\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 504\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 505\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 506\u001b[0m \u001b[1;32mexcept\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 507\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mtraceback\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m\u001b[0m in \u001b[0;36mcomp_theta_prob\u001b[1;34m(self, theta, efficiency)\u001b[0m\n\u001b[0;32m 114\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 115\u001b[0m \u001b[1;31m# Get value of rate posterior at this rate\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 116\u001b[1;33m \u001b[0mbns_rate_pdf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mscenario\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcomp_bns_rate_pdf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbns_rate\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 117\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 118\u001b[0m \u001b[1;31m# Compute jacobian\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/repositories/grbeams/grbeams/scenarios.pyc\u001b[0m in \u001b[0;36mcomp_bns_rate_pdf\u001b[1;34m(self, rate)\u001b[0m\n\u001b[0;32m 66\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 67\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mcomp_bns_rate_pdf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrate\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 68\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbns_prior\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpdf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrate\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 69\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 70\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mcomp_grb_rate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mefficiency\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtheta\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbns_rate\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/repositories/grbeams/grbeams/scenarios.pyc\u001b[0m in \u001b[0;36mpdf\u001b[1;34m(self, rate)\u001b[0m\n\u001b[0;32m 34\u001b[0m \u001b[0mrate\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrate\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mto\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmegaparsec\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mu\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0myear\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 35\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mrate\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;33m:\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 36\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mrate\u001b[0m \u001b[1;33m>\u001b[0m \u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrates\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 37\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minterp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrate\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrates\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpdf_data\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 38\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/astropy/units/quantity.pyc\u001b[0m in \u001b[0;36m__array_prepare__\u001b[1;34m(self, obj, context)\u001b[0m\n\u001b[0;32m 407\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# normal case: set up output as a Quantity\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 408\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 409\u001b[1;33m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_new_view\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresult_unit\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 410\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 411\u001b[0m \u001b[1;31m# We now need to treat the case where the inputs have to be scaled -\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/astropy/units/quantity.pyc\u001b[0m in \u001b[0;36m_new_view\u001b[1;34m(self, obj, unit)\u001b[0m\n\u001b[0;32m 587\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 588\u001b[0m \u001b[0mview\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mobj\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mview\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msubclass\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 589\u001b[1;33m \u001b[0mview\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__array_finalize__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 590\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0munit\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 591\u001b[0m \u001b[0mview\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_unit\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0munit\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/daniel/.virtualenvs/grbeams/lib/python2.7/site-packages/astropy/units/quantity.pyc\u001b[0m in \u001b[0;36m__array_finalize__\u001b[1;34m(self, obj)\u001b[0m\n\u001b[0;32m 281\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 282\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__array_finalize__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 283\u001b[1;33m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_unit\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'_unit'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 284\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 285\u001b[0m \u001b[1;31m# Copy info if the original had `info` defined. Because of the way the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;31mKeyboardInterrupt\u001b[0m: " + "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mkde\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mthetapos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_theta_pdf_kde\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;32m/scratch/aries/daniel/repositories/grbeams/grbeams/beamingangle.py\u001b[0m in \u001b[0;36mget_theta_pdf_kde\u001b[1;34m(self, bandwidth)\u001b[0m\n\u001b[0;32m 58\u001b[0m algorithm='kd_tree') \n\u001b[0;32m 59\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtheta_bounds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtheta_median\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtheta_posmax\u001b[0m \u001b[1;33m=\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 60\u001b[1;33m \u001b[0mcharacterise_dist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtheta_grid\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtheta_pdf_kde\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0.95\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 61\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 62\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0msample_theta_posterior\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnburnin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m100\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnsamp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m500\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnwalkers\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m100\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;31mNameError\u001b[0m: global name 'characterise_dist' is not defined" ] } ], "source": [ - "thetapos.sample_theta_posterior()" + "kde = thetapos.get_theta_pdf_kde()" ] }, { @@ -417,7 +204,14 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.10" + "version": "2.7.6" + }, + "latex_envs": { + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 0 } }, "nbformat": 4, diff --git a/notebooks/sGRB Beaming.ipynb b/notebooks/sGRB Beaming.ipynb new file mode 100644 index 0000000..e71c53c --- /dev/null +++ b/notebooks/sGRB Beaming.ipynb @@ -0,0 +1,1137 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import astropy.units as u\n", + "import astropy.constants as c" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "plt.style.use('/home/daniel/repositories/burst-style/burst.mplstyle')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is widely believed that short Gamma-ray bursts (sGRB) are produced by compact binary coalescence (CBCs), and observations suggest that the coalescing systems likely to produce them are binary neutron star systems (BNS systems), or neutron star-black hole binaries (NSBH systems). At the full design sensitivity of advanced LIGO these may be visible to a distance of 400 Mpc for NSNS and 1 Gpc for NSBH." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Given this link between sGRBs and CBC events, it is interesting to ask if the beaming angle of the radiation from the GRB can be inferred from gravitational wave (GW) observations." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If the sGRB population has a distribution of beaming angles then the observed rate is related to the BNS coalescence rate, $\\mathcal{R}$, by \n", + "\\\\[ \\mathcal{R}_{\\text{GRB}} = \\epsilon \\mathcal{R} \\langle 1 - \\cos(\\theta) \\rangle \\\\]\n", + "\n", + "\n", + "| Symbol | Name | Note |\n", + "| ------ | ---- | ---- |\n", + "| $\\epsilon$ | efficiency | the probability that the BNS produces a GRB |\n", + "| $\\mathcal{R}_{\\text{GRB}}$ | GRB rate | the observed rate of GRBs |\n", + "| $\\mathcal{R}$ | BNS merger rate | the observed rate of BNS CBCs |\n", + "| $\\theta$ | beaming angle | the angle which radiation is beamed from the GRB |\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "EM observations suggest that beaming angles have a distribution, so the raltive rates of sGRBs and CBCs should provide the mean of that population, $\\langle \\theta \\rangle$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "@u.quantity_input(bns_rate=1/u.year, theta=u.degree)\n", + "def rate_grb(efficiency, theta, bns_rate):\n", + " \"\"\"\n", + " Calculate the rate of gamma ray bursts from the rate of BNS coalescence.\n", + " \n", + " efficiency : float\n", + " The efficiency of GRB production.\n", + " theta : float\n", + " The beaming angle\n", + " bns_rate : The rate of BNS mergers\n", + " \n", + " \"\"\"\n", + " return efficiency * bns_rate * (1- np.cos(theta))\n", + "\n", + "@u.quantity_input(grb_rate=1/u.year, theta=u.degree)\n", + "def rate_cbc(efficiency, theta, grb_rate):\n", + " \"\"\"\n", + " Calculate the rate of CBCs from the GRB rate.\n", + " \"\"\"\n", + " return grb_rate / (efficiency*(1.0 - np.cos(theta)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The sGRB population mean\n", + "\n", + "Here we show the impact of BNS mergers on sGRBs, with the following proceedure:\n", + "\n", + "* We set the number of observed GRBs to zero: $N_{\\text{GRB}} = 0$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "Ngrb = 0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* We draw $N_{\\text{BNS}}$ values for orbital inclination, $\\iota$ from a distribution uniform in $\\cos(\\iota)$ over the range $[0,1]$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "Nbns = 1e5 \n", + "iotas = np.arccos( 0 + 1*np.random.rand(Nbns) )\n", + "# to degrees\n", + "iotas *= 180/np.pi" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* For each $\\iota$ draw a $\\theta$ from a distribution over $(0,90]^\\circ$\n", + "* Iff $\\iota < \\theta$ the event is observable: increment $N_{\\text{GRB}}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "theta_mu=np.arange(5,35,5)\n", + "theta_sigma=np.arange(1,16)\n", + "theta_low=0.01\n", + "theta_upp=90" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def truncparms(low,upp,mu,sigma):\n", + " a = (low - mu) / sigma\n", + " b = (upp - mu) / sigma\n", + " return a, b" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to calculate the efficiency we take \n", + "\n", + "\\\\[ \\epsilon = \\frac{ N_{\\text{GRB}} + 1 }{N_{\\text{BNS}} + 2 } \\\\]\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def compute_efficiency(k,N,b=True):\n", + "\n", + " if b:\n", + " # Bayesian treatment\n", + " epsilon=(k+1)/(N+2)\n", + " stdev_epsilon=1.64*np.sqrt(epsilon*(1-epsilon)/(N+3))\n", + " else:\n", + " # Binomial treatment\n", + " if N==0:\n", + " epsilon=0.0\n", + " stdev_epsilon=0.0\n", + " else:\n", + " epsilon=k/N\n", + " stdev_epsilon=1.64*np.sqrt(epsilon*(1-epsilon)/N)\n", + " return (epsilon,stdev_epsilon)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll draw all of the beaming angles from a truncated normal distribution, and we'll run the simulation for a variety of differently shaped distribution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import scipy.stats as stats\n", + "\n", + "FracGRB=np.zeros((len(theta_mu),len(theta_sigma)))\n", + "deltaFracGRB=np.zeros((len(theta_mu),len(theta_sigma)))\n", + "\n", + "for m,mu, in enumerate(theta_mu):\n", + " for s,sigma in enumerate(theta_sigma):\n", + " a, b = truncparms(theta_low, theta_upp, mu, sigma)\n", + " thetas = stats.truncnorm.rvs(a, b, loc=mu, scale=sigma, size=Nbns)\n", + " # Count GRBs\n", + " Ngrbs=sum(iotas <= thetas)\n", + " \n", + " FracGRB[m,s], deltaFracGRB[m,s] = compute_efficiency(Ngrbs,Nbns)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A population of sGRB beaming angles with a large mean but a narrow width is indistinguishable from the converse, as can be seen in the plot below, where $\\mathcal{N}(15,8)$ and $\\mathcal{N}(10,12)$ both give the same number ratio (see the light grey line)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "linestyles=['-', '-', '-']\n", + "markers=['s', '^', 'o']\n", + "\n", + "f, ax = plt.subplots(1,1)\n", + "\n", + "ax.hlines(FracGRB[1,11]*100,0,16, lw=1, color='gray', alpha=0.5)\n", + "for i in range(len(theta_mu)):\n", + " ax.plot(theta_sigma, FracGRB[i]*100, #deltaFracGRB[i]*100, \n", + " linestyle=linestyles[i%3], lw=2,# c='k',\n", + " #marker=markers[i%3], markeredgecolor='white', \n", + " #label=r\"$\\langle \\theta \\rangle = {}$\".format(theta_mu[i])\n", + " )\n", + " ax.annotate(\"{}\".format(theta_mu[i]), (15.1, .1+FracGRB[i][-1]*100))\n", + "ax.set_ylim(0,18.5)\n", + "ax.annotate(r\"$\\langle \\theta \\rangle$\", (15, 17.6))\n", + "ax.set_ylabel(r\"$N_{\\rm{GRB}} / N_{\\rm{BNS}}$ [%]\")\n", + "ax.set_xlabel(r\"$\\sigma_{\\theta}$ [degree]\")\n", + "plt.tight_layout()\n", + "plt.savefig('color_relativenumber.pdf')\n", + "\n", + "txt = \"\"\"Expected relative numbers of observed GRBs and binary coalescences for different \n", + "distributions on the GRB beaming angle. Lines in the figure correspond to jet angle\n", + "population means, while the x-axis shows the width of the distribution. All distributions are \n", + "Gaussian, truncated at (0, 90] degrees.\"\"\"\n", + "f.text(.1,-.15,txt)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This has the effect of a wide distribution of beaming angles producing an under-estimate of the beaming angle, so population-based constraints must be regarded as upper-bounds." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Rates to Beaming Angle" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The general method to move from a rate to a beaming angle is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate the posterior distribution on the BNS merger rate from observed GW signals" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Gravitational wave detectors identify discrete 'events' which are characterised by their *network signal-to-noise ratios*, $\\rho_{\\mathrm{c}}$\n", + "\n", + "We begin by constructing the posterior on the signal rate; the measured rate, $r$ of events has two components, the signal rate, $s$, and the background rate, $b$:\n", + "\n", + "\\\\[ r = s+b\\\\]\n", + "\n", + "Typically a threshold of $\\rho \\ge 12$ is applied to analysis, placing $b = 10^{-2}$ per year, so only the signal rate needs to be inferred.\n", + "\n", + "Taking a uniform prior on $s$, and assuming the events are generated by a Poisson process, it can be shown that the rate posterior depends only on the number of detected signals, $n$, the observation time, $T$, the background rate, $b$, and the horizon distance, $\\mathcal{D}$ (the maximum distance at from which a detection can be made).\n", + "\n", + "The posterior will then be\n", + "\n", + "\\\\[\n", + " p(s|n,b,I) = C \\frac{T[(s+b) T]^n e^{-(s+b)T}}{n!} \n", + "\\\\]\n", + "\n", + "for\n", + "\\begin{align*}\n", + " C^{-1} &= \\frac{e^{-bT}}{n!} \\int_0^\\infty {\\rm d}(sT)(s+b)^n T^n e^{-sT} \\\\\n", + " &= \\sum_{i=0}^n \\frac{ (bT)^i e^{-bT} }{i!} \n", + "\\end{align*}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "@u.quantity_input(s=1/u.year, \n", + " b=1/u.year,\n", + " T=1*u.year)\n", + "def signal_posterior(s, b, n, T):\n", + " \"\"\"\n", + " The event rate posterior assuming that a Poisson process underlies \n", + " the production of gravitational wave events.\n", + " \n", + " s : float\n", + " The event rate\n", + " n : float\n", + " The number of events observed over the timespan of the observing run\n", + " b : float\n", + " The background rate of events\n", + " T : float\n", + " The total time of the observing run\n", + " \"\"\"\n", + " # Calculate the inverse of C\n", + " C_inv = np.sum([ (b*T)**i * np.exp(-b*T) / np.math.factorial(i) \n", + " for i in xrange(n) ])\n", + " # Then invert it\n", + " C = 1. / C_inv\n", + " # Calculate the posterior piece-wise\n", + " p1 = T * ((s+b)*T)**n * np.exp( - (s+b) * T)\n", + " p2 = np.math.factorial(n)\n", + " return C * p1 / p2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's have a quick look at a plot of this." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "s = np.linspace(0,10, 11)/u.year\n", + "b = 10 / u.year #np.linspace(0,10, 11)/u.year\n", + "T = 1 * u.year\n", + "plt.plot(s, signal_posterior(s, b, 10, T))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll want to express the rate in terms of the number of mergers per *milky-way equivalent galaxy*, or MWEG. We can approximate the number of MWEGS inside the sensitivity range of the detectors as\n", + "\\\\[\n", + " N_{\\mathrm{G}} = \\frac{4}{3} \\pi \\left( \\frac{\\mathcal{D}}{\\mathrm{Mpc}} \\right)^3 \n", + " (2.26)^{-3}\n", + " (0.0116)\n", + "\\\\]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "@u.quantity_input(horizon=u.megaparsec)\n", + "def n_mweg(horizon):\n", + " mweg1 = (4. / 3. ) * np.pi \n", + " mweg2 = (horizon / 1/u.megaparsec)**3\n", + " mweg3 = (2.26)**(-3) * 0.0116\n", + " return mweg1 * mweg2 * mweg3" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This let's us reach the posterior on the binary coalescence rate:\n", + "\\begin{align*}\n", + " p(\\mathcal{R} | s,n,T,b,\\mathcal{D}) &= p(s|n,T,b) \\left| \\frac{\\mathrm{d}s}{\\mathrm{d}\\mathcal{R}} \\right| \\\\\n", + " &= N_{\\mathrm{G}} \\mathcal{D} \\ p(s |n,T,b)\n", + "\\end{align*}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "@u.quantity_input(s=1/u.year, \n", + " b=1/u.year,\n", + " T=1*u.year, \n", + " horizon=u.megaparsec)\n", + "def rate_posterior(s, b, n, T, horizon):\n", + " \"\"\"\n", + " \n", + " \"\"\"\n", + " nmweg = n_mweg(horizon)\n", + " post = signal_posterior(s,b,n, T)\n", + " return nmweg * post" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "s = np.linspace(0,10, 11)/u.year\n", + "b = 10 / u.year #np.linspace(0,10, 11)/u.year\n", + "T = 1 * u.year\n", + "horizon = 100*u.megaparsec\n", + "plt.plot(s, rate_posterior(s, b, 10, T, horizon))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These quantities define our *detection scenarios*.\n", + "\n", + "\n", + "| Scenario name | $T$ | $\\mathcal{D}$ | $V$ | $\\mathcal{R}$ | $n$ |\n", + "| ------------- | --- | ------------- | ------ | ------------- | --- |\n", + "| 2016 | 0.5 | 80 - 120 | 1.05e6 | 1e-6 | 1.3 |\n", + "| 2022+ | 1.0 | 200 | 4.00e7 | 1e-6 | 40 |" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Transform from the rate posterior to an angle posterior" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Inferences of the GRB beaming angle are made from the posterior on the beaming angle, so we need to convert a posterior on the rate prior into a posterior on the angle. To do this we use a Jacobian.\n", + "\n", + "\\\\[ p(\\theta, \\epsilon) = p(\\mathcal{R}, \\epsilon) \\left|\\left| \\frac{ \\partial(\\mathcal{R},\\epsilon) }{ \\partial (\\theta, \\epsilon) } \\right|\\right| \\\\]\n", + "\n", + "We assume that the beaming angle and the rate are probabilistically independent, so\n", + "\n", + "\\\\[ p(\\theta) = \\frac{2 \\mathcal{R}_{\\rm GRB} \\sin \\theta \\ p(\\mathcal{R}) }{( \\cos \\theta - 1)^2 } \\int \\frac{p(\\epsilon) {\\rm d} \\epsilon }{\\epsilon} \\\\]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def jacobian(efficiency, theta, grb_rate):\n", + " denom=efficiency*(np.cos(theta * np.pi/180)-1)\n", + " return abs(2.0*grb_rate * np.sin(theta * np.pi / 180.0) /\n", + " (denom*denom) )" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "So for a given efficiency, $\\epsilon$ and beaming angle $\\theta$ the CBC rate needs to be found, then the rate posterior for this calculated, and finally rotated using the Jacobian." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from grbeams.distributions import UniformDistribution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def p_theta(theta, efficiency, grb_rate, theta_range=(0*u.degree,90*u.degree)):\n", + " if (theta>=min(theta_range)) and (theta" + ] + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAESCAYAAADjS5I+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VPXdLvBn79lzSSCZ3IiUoihY610xBJNgwkgMjamA\nvKWseml939WzcL3W1UpYq6c9R1qsfVf7HpWmtVbK0bqW2kWXWJRKoxUSkigkElgcg4CgKCIQICSZ\nhLllMpfzRzJjZjKZzEz2zJ695/n8lcxl7+9mFvPkd92C3+/3g4iIaAKi0gUQEVF6Y1AQEVFUDAoi\nIoqKQUFERFExKIiIKCoGBRERRcWgICKiqBgUREQUlaR0ARPZv38/3n77bfh8Pvzyl79UuhwiooyV\n8hZFd3c3SkpKkJ2dDZ/PF3y8vr4eVVVVWLt2LQBg+/btWL9+PYqKinD+/PlUl0lERKNSHhSFhYVo\nbm5GWVlZ8LGDBw/Cbrejra0Nbrcb+/fvDz7HHUaIiJSV8q4ng8EAg8EQEgAdHR2oqakBAFRXV6Oj\nowPLly/Hk08+ieHhYVx22WWpLpOIiEalxRiF1WrFvHnzAABmsxlHjhxBaWkpSktLo76vqakpFeUR\nEWlOdXV1zK9VLCgEQQj+bDabMTg4CAAYHBxEXl5ezMe5btdWnDl7Fit3HcDpM2ewYQOwYYPMxSpk\n3bp1eOaZZ5QuI2m0fH1avjaA16d28f6Rrdj0WL/fH+x+Ki8vDxa+a9eukPGLWJh0OrhcLgDAE0/I\nWycRUaZLeVB4PB7U1NSgq6sLtbW16OzsxPz582E0GlFVVQVJkrBgwYKYj+fz+2GURLiGhpJYNRGR\n+rW0tGBDAl0uKe96kiQJO3fuHPd4Q0NDQscTBQHGMS0KLSkvL1e6hKTS8vVp+doAXp9aWSwWWCwW\n9XQ9yUkSRy7D4/EoXIm8KioqlC4hqbR8fVq+NoDXl2k0ERQAYDKZNNmqICJSmuqDwjc6IG4yGuFy\nucDdPoiIIkt0jEL1QSGOTrMNtCi0MjWWiEhuFoslM4MiwDjaoiAiInlpJig4RkFElBxpsYXHVPj8\nfggATCYjHA6H0uUQEaWtlpYWtLS0oLKyMq73qb5FERijyM6exqAgIooi48copmVnw263czCbiEhm\n2gmKaSNBwb2eiIjkpZmgyM6eBrvdrnQZRESao5nB7EDXExERRZbxg9mBriciIoos4wez2fVERJQc\nmgmKadOyYbPZuNcTEZHMtBMUnB5LRJQUmgkKdj0RESWHdmY9cTCbiCgqznpii4KIKCrOemKLgogo\nKTQTFBzMJiJKDk0Ehd9ug66/B5f6+rjXExGRzLQRFD3n4G94ArbeHqVLISLSHE0EBQBkSSKcwx6l\nyyAi0hxNTI8FgGk6HRwer8LVEBGlr4yfHqsXBUiioHA1RETpK9HpsapvUYyVa9Djvv+4AKBY6VKI\niDRD9S2KsXINevzH908oXQYRkaZoKihyDHoMDA4qXQYRkaZoLygGBpQug4hIUzQVFLkGiS0KIiKZ\naSoocvR6DDIoiIhkpa2gMOixbfttSpdBRKQpmgqKXIMeza0WpcsgItIU1a+jCKzMBoAcg+ovh4go\naTJ+ZTYw0qIgIqLIMv7GRQCQbzIoXQIRkeZoKigKTEalSyAi0hyNBYUBubm/U7oMIiJN0VRQFJqM\n8HrXK10GEZGmaCoopuslDA8Pw+l0Kl0KEZFmaCooBEFAUWEhenp4S1QiIrloKigAYEZREYOCiEhG\nmgsKtiiIiOSluaA4d/4/GRRERDLSXFAcProaFy9eVLoMIiLN0FxQAGCLgohIRqrfRW/spoABDAoi\novG4KeAYDAoiovG4KeAYDAoiIvloLige/uGX6O3tVboMIiLN0FxQ/K+f2hkUREQy0lxQFOTno7+/\nHz6fT+lSiIg0QXNBIUkScnJyYLValS6FiEgTtBUUtkE4jx9B/vRp6D75udLVEBFpgqaCwj9ohe1X\na2F22XDx/HmlyyEi0gRNBQUAbDy+BvlGAy72cUCbiEgOmguKhk/XIM+gh9U6oHQpRESaoLmgAACz\nQY9+DmYTEclCm0FhZFAQEclFk0GRZ9DDOsCuJyIiOWgyKMxGPddREBHJRHNB8djVm0fGKNiiICKS\nheaCov6akaDgrCciInloLiiAkVlP1gF2PRERySFtg+LQoUNYuXIljh8/Hvd7R2Y9sUVBRCSHlAVF\nd3c3SkpKkJ2dHbKza319PaqqqrB27dqQ1990001YuXJlQufirCciIvmkLCgKCwvR3NyMsrKy4GMH\nDx6E3W5HW1sb3G43Dhw4gB07dmD79u0AAL/fD3+Ee2JPxqQT4fP54HK5ZKufiChTpSwoDAYDzGZz\nyBd/R0cHampqAADV1dVob2/HPffcgxUrVuDUqVN499138corr2B4eDjm82w8vgaCICDPbEZ/f7/s\n10FElGkkJU9utVoxb948AIDZbMaRI0eCz11xxRX461//OukxfL7QFkfDp2uwrmI7cnNycOzYsYRa\nJOni0qVLOHv2rNJlJI2Wr0/L1wbw+tRm7969aG9vD/5eV1cX1/tTHhSCIAR/NpvNGBwcBAAMDg4i\nLy8v7uOJogBvhMdnFBVBp9Nh1qxZiZaquLNnz6q6/slo+fq0fG0Ar09tVq1ahVWrVgV/b2pqiuv9\nKZ/1NHbcoby8PFjwrl27QsYvpqqosBA9PT2yHY+IKFOlLCg8Hg9qamrQ1dWF2tpadHZ2Yv78+TAa\njaiqqoIkSViwYEHcx/VN0LU0o6gIFy5cmGrZRESa0dLSgg0bNsT9vpR1PUmShJ07d457vKGhYUrH\nFYUJup5mMCiIiMayWCywWCzp3/WUbI9dvRkAUFxUyKAgIpKB5oKi/pqRoGDXExGRPBSdHisHjlEQ\nEcWmpaUFLS0tqKysjOt9qm9RiGOm247FoCAiCmWxWBIazFZ9UExkRlERp8cSEclAs0FRkJ8Pq9UK\nj8ejdClERKqmuaDYeHwNAECn0yE/Px+9vb0KV0REpG6qD4rwweyGT9cA7iEMne9G3rRpOHfqC4Uq\nIyJKL4kuuFN9UEQazPYPWnHpf/4P5FzqQ9+5cwpURUSUfjiYHUGuQY9+3hKViGhKNB0UZr2EgYFB\npcsgIlI1TQdFrkHiLVGJiKZIcyuzA3s9AYCZ984mIgriyuxRgb2eAMBskGBl1xMREQAOZkdk1kuw\ncjCbiGhKtB0UBj0GBtmiICKaCk0HRa5eQj/HKIiIpkT1QTHRNuMA8PVpJhw9dpz7PRERgSuzgwJ7\nPQHAnOlZuPzrX8ef//znVJdFRJR2Eh3Mjnl67MDAAA4ePIjz588jJycHJSUluOyyy+I+YbI1fLom\nOPNJFAQ8uPq76OrqUrgqIiL1mjQo+vr68M4770Cv1+Pqq6/G5ZdfDpvNhpaWFjgcDixevBhz585N\nRa0JmT1rFpr2titdBhGRasUUFPfff/+4x2+55RYAwMcffyx/VTLKz8tDf3+/0mUQEanWpGMUV199\nNbZs2TLh89dee62sBcktPz8PfX19SpdBRKRaMY1RbNy4Ea+++iquuOIKVFRUoKKiAn19fSgtLU12\nfVOWn8egICKaipiCYuvWrbjyyitx+vRp7N27F88++yyamppw6NChZNcXt7F7PQFAvtmM/v5++P1+\nCBHuXUFERNHFFBRXXnklAGD27NlYvXo1Vq9ejY6OjmTWFbPwdRRj93oCAKPRCIPBAJvNhpycnFSW\nRkSUVlK+KWBZWVmib5VVpDvchSsoKGD3ExFlvJRvCvjUU0+hs7MThw8fTvQQKVNYWIje3l6lyyAi\nUqWE70dRWVmJ0tJSVfylzqAgIkpcTEHh9Xrx+uuv48svv0RlZSVuv/32YNdTQUFBUguUQ0FBAYOC\niChBMXU9rV+/HkajEfPmzcPWrVuxZs0auN3uZNeWkLF7Pfk9HrgH+pEjAhfOnFawKiIi9YqpRXHr\nrbfi3nvvBQCsXLkSfX19+OMf/4j6+vqkFpeIsXs9wWHDpXX/jmmfXUDP7MuVLYyISKVialGEtx4K\nCgpQXFyclIKSId+gRy+38SAiSkhMLYrNmzfjjTfeQF5eHsrKylBRUYHp06cnuzbZFBn1ONxzUeky\niIhUKaageO211zBz5kwMDg6io6MDW7duxT//+c9gd5SSot24KKDYZMC5CxdSUA0RUfpKdMFdTEEx\nc+ZMAEBubi6WLl2KpUuXYsWKFfFXmQSiIMA7yWuKTQacP82gIKLMZrFYYLFY0NTUFNf7Ym5RDA8P\nA0DIfklHjhyBwWDAd7/73bhOmkzhez0BQHGWAed7ehSohohI/WIKitWrVye7DtmE7/UEADmSDu7h\nYTgcDmRnZytQFRGResUUFK+//jp2794NURRxxx13YPny5cjKykp2bbIRBAGF+fno7e1lUBARxSmm\n6bHZ2dl47rnncP3118Nut+MnP/kJWltbk12brAoLCnDxImc+ERHFK6agKCoqwvr169Hd3Y26ujps\n2rQJZ86cSXZtsiosyGdQEBElIKaup4ULF+K2225DY2MjNm7cCJvNhrvvvjvZtcmqMJ8tCiKiREwa\nFB0dHSgrK4MkSVi+fDmWL18e8nxnZ2da3RJ14/E1EQe02aIgIkrMpEFxww034MUXX4TL5cLll1+O\nnJwc2O12nDp1CtnZ2bjnnntSUWfMQvZ6GoNjFEREiZk0KHJycvDDH/4QPp8Px44dg9VqRXFxMWpr\nayFJCd/OIuUKC/LxWfd5pcsgIlKdmO9wJ4oiWlpa8Ic//AHPP/887HZ7MuuSHccoiIgSE9etUK+/\n/nps2bIF//Vf/4WXX345WTUlRRG7noiIEhJX31FjYyN27NiBq666CkNDQ/B4PDhz5gzmzJmTrPom\nFcumgAAHs4mIkropYMCKFStw1VVXwev14oMPPsATTzyBjo4O7Ny5M66Tyil8U8BIez0BHMwmIkrq\npoABFRUVwZ9nz56N73znO7iQZtt3R5rxBAAFeSMtCr/fH7KxIRERRRfXGEUkqrjTnc8HSQRMBgPO\nf/GF0tUQEanKlINCFYZcuPTj+zFT8OLzjw8rXQ0RkapkRlCMmmUy4Ez3OaXLICJSlYwKiq9lGXCm\nu1vpMoiIVEVzQbHx+JoJn2NQEBHFT3NB0fDpxEEx02RA9/n0mqVFRJTuNBcU0Vxm1KP7PPd7IiKK\nR0YFxcwsA4OCiChOGRUUlxn16D7HoCAiikdGBUW+QYLd6YTT6VS6FCIi1dBcUEy01xMACIKAK2Z/\nHUeOHElhRURE6qa5oJhor6eA+/9tJTZt2pSiaoiI1E9zQTGZe79dh3fffVfpMoiIVCPjguKqKy5H\nd3c3hoaGlC6FiEgVMi4oJEnCrFmzcPr0aaVLISJShYwLCgCYM2cOvvzyS6XLICJShbQNiqamJvz2\nt7/FY489Ftf7ou31FPC1r30N3dzziYgoJikLiu7ubpSUlCA7Oxs+ny/4eH19PaqqqrB27dqQ11dX\nV+NnP/tZyGtjEW2vp4CZM2fi3DluN05EFIuUBUVhYSGam5tRVlYWfOzgwYOw2+1oa2uD2+3GgQMH\nsGPHDmzfvh0A0NDQgAcffFD2WmbOnIlt27Zh2bJlsh+biEhr4rpn9lQYDAYYDAb4/f7gYx0dHaip\nqQEw0oJob2/Ho48+CgDYvHkz9u/fD1EUsXDhwgmP6/P5J3wuErvdAZPJhPfffx8AcPbs2XgvJWUu\nXbqU1vVNlZavT8vXBvD61Gbv3r1ob28P/l5XVxfX+1MWFJFYrVbMmzcPAGA2m0NWTK9ZswZr1kze\njSSKArxxnHPatGxcd911wd/z8vKQnZ0dxxFS5+zZs5g1a5bSZSSNlq9Py9cG8PrUZtWqVVi1alXw\n96amprjen/LBbEEQgj+bzWYMDg4CAAYHB5GXl5fck/v9GHY4MN1lCz504QLvT0FEFE3Kg8Lv9we7\nn8rLy4PJtmvXrpDxi0RF2+sJPh9s6/4dWc/+KvhQb2/vlM9JRKRlKQsKj8eDmpoadHV1oba2Fp2d\nnZg/fz6MRiOqqqogSRIWLFgQ93F9/tAxisn2egKAAsNIj9s3v/lNBgURZYyWlhZs2LAh7velbIxC\nkiTs3Llz3OMNDQ1TOq4oxDdGAQA6QcBnnR342VMbGRRElDEsFgssFkv6j1Gki6LCQhQWFqKvr0/p\nUoiI0lrGBgUwsraDLQoiouhUHxThYxTxKCgoYFAQUcZIdIxC9UEhjpluC8S211MAu56IKJNYLJbM\nDIpwsez1FFBcXMzNAYmIJqG5oIjHddddh6NHjypdBhFRWlN9UCQ6RuEZ6Id5oBeDAwPo7++XuSoi\novST9usokiWRdRQAYPvf/wkAuGbuXHz88ccoLy+XtzAiojTDdRQJuvqqK3H8+HGlyyAiSluaC4qo\nez1FMGf21/H5558nqRoiIvXTXFDEstfTWPlmM6xWa5KqISJSP9UHxVQW3AGAOTeXg9lElBG44C5B\neQwKIsoQXHCXoOzsLDidTqXLICJKWxkfFCaDEU6nE42NjVi7dq3S5RARpR3NBUU8ez0BgNFohMvl\nwgsvvDDle2MQEWmR5oIinr2eACDLNBIUoqi5fwoiIlmofmX2VGc9GY0jXU86nU6mioiI0lNLSwta\nWlpQWVkZ1/tU/2f0VGc9ZRlNbFEQUUbgrKcEmYyGkFlPbrdbwWqIiNJPxgeFpNfD4/FgeHgYAGC3\n2xWuiIgovWguKOLd60nS6TA8PBxsSdhstmSURUSkWpoLinj3etJLUkiLgkFBRBRKc0ERL2k0KIaG\nhgAwKIiIwmX89FhBEKDT6eB0OiGKIoOCiDSL02MTNHT4ICRRhM1qRUFBAYOCiDSL02MT5Pq/z0Dy\neWF32JGfny9rUJw6dUq2YxERKUVzQRHvXk8AIIkCnE6XrEExPDyMOXPm4OTJk7Icj4hIKZoLinj3\negIASRDgdLlk7Xrq7e0FAN7rgohUT3NBkQhJEOBwOWVtUbhcLgDApUuXZDkeEZFSGBQYCYrhYQ+y\nbQPoO31almMGptsyKIhI7VQ/PVYOkjgyc2r60YOwTZsuyzHZoiAirWCLAoB+dIqtWa+DzeGQ5ZiB\nFkUgMIiI1Er1QRG+4C7evZ4AQDe6FMOs18HuiHz/7EuXLiErKwtPPvkkamtrJz1mICgCW4MQESmt\npaUloXUUqu96EgUB3jG/x7vXEwDoxUCLQoJ9ghbFvn374HK58Itf/GLS4x07dgzHjh0DwKAgovRh\nsVhgsVjQ1NQU1/tUHxRy0I3pepooKE6cODHpcfbt24dFixbB4/EEH+P9LYhI7VTf9SSHwCYg0yUd\nbA4nBgYGxr2mr69v0uN8+OGHISEBsEVBROrHoBhjmk5Eb38/8vLy8Mknn4Q819/fP+ntUoUI+04x\nKIhI7RgU+OoLPlvS4VxPDwDgs88+C3lNf38/5s2bF/U4gQHssdj1RERqp7mgSGSvJ//ozKlcSRd8\n7MyZMyGv6evrw9y5c6MeJ9KqbrYoiEjtNBcUiez1FGDSffXPYbVaQ56LpUUR6X7bDAoiUjvNBYVc\nBgcHQ36PpUXhCJsxZTQaNdX1VF9fn9AcbCJSN06PDdO85VU0Hvh/44Kiv79/0qAID4WcnBxNtSh+\n97vfwWw2MyyIMgxbFGFummbAbNGL3jOhmwP29/aiyHoh+LvP5xv33vBQmD59uqxBsXPnTmzZskW2\n4yUiOztb0fMTUeqxRRHGuem/oe/uw+CMr8YjPB4P7E4n8rY8H3zM6/WOmy4bHgq5ubmydj098MAD\n6OnpwX333SfbMeOVlZWl2LmJSBmaa1EkstdT+H23cyQdLtm+Gpi2Wq3InTYNZkmHV377axiNRni9\n3vDDRAyKqbQoPB5PSBdYOtzPO1JLioi0TfVBEb4pYCJ7PYWvk8vSiXCM2fW1r68P+bk5EAQB315c\nCUmSxq3ABuRvUTzxxBMwm83B3yOdM9UYFETqleimgKoPivDWQELHQOgxTKII55iguHDhAory84O/\n63S6mIJC8nrgvDQ47nWxOnz4cGidk6wMTwUGBZF6WSyWzAwKOYhhWWPSiXC5RlZZf/nll9iyZQtm\n5OcFn5ckKaauJ++HnXD3XUy4rqKiIgBffTnrdLpoL08Jf1gLjoi0j0GB8f8IJp0I5+h2HM888wz+\n9Kc/xdSiCO9mEgXAO4W/wAMtiO7ubgCR95JKNbYoiDIPgwLjv4BNogjX0EjXkySNTAwrMOcGn491\njGLeNNOUvlidzpGbKHV1dQFg1xMRKUP5bx6ZJbLXU/jf6SadEGxRBL78K+bfGnw+WtfTr3/9a7z0\n0ksYOH4EpQXT4fUm/sXqcDhQUVGB9957DwCDgoiUofw3j8wS2espfOaUSRThdI60KGw2G1544QVU\nly0MPh9tMLt8/i34t8oKeIc90AnClLqenE4namtr0dbWBiA9giIdaiCi1OKCOwC/v/Uq9Lu/+uI3\niAI8Xi+8Xi/sdjumT58e8vpoLQrBdgkDP74fwMi9uMNDKB7Dw8MoKSnBs88+CyA9BpIZFESZh0EB\nYKbJgJkmQ/B3QRCQZTTC6XTCbreP27Yi2hhFYEwDAARMrUXh8XgwY8YMWK1W+P3+tOj2SYcBdSJK\nLf55OIGsLBOcTieGh4dhMBhCnovW9aQfExQ6YWp9+h6PB9OnT4dOp4PT6WSLgogUwf/1E8gyGuFw\nOODxeMatX5io68ntdkOvHxsUU29RSJKEvLw8WK1WRVsUgZBii4Io82guKBLZ6ykS02jXU+DLeqxY\nu56muo7C6/VCkiRkZWXB5XIpGhSBYEyHVg0RpZbmgiKRvZ4iMZlGup7cDjtcpz6Dz9oHAPBZ+wH3\nEAZPnwp5vePcGbgdDoiur25epBME+BKcHutxu+EZcsFx4jhMRiNcLpeiX9KBacLpsN8UEaUWB7Mn\nEBjMHnbY4dr0f+DIH5n55Hj2SeD05xjqDr2n9nBvD9zWPrh//yvAqAcA6KYwmO3zeuGxXYJ7619g\nNJngcrlC/qpPdRdQICgidbkRkbZprkUhl0DXk9frgxT2payDEHl6rN8P/ZjXTrXryePzQxIFmEaD\nIvBlrUQXFIOCKHMxKCaQZRodo/B6oQsLCkkU4PWN/8L0+PwhoTLlwWyfDzpBhNFgCN6P22AwKNL9\n4/V6odfrGRREGYhBMQG9pIfb7YbX54UurJdHJwCeCF+YHr8felEIed2Upsf6fNDrBBhNJtjtduh0\nugkH0pPN6/UqFlJEpCzNBUUiez1FYtCPBEWkFoVO+KrryeVyYe/evQCA4bAWhTjFFoXX74dOEGEy\nmeBwOGIOikOHDuHQoUMJnzcSn8/HFgVRhtJcUCSy11Mker000qLweMePUQhCcLO/F198EYsWLRpZ\nOQ2EtD6m3vU0Mkah1+vhcrkgimJMQVFbW4ubb7454fNGEmhRMCiIMk/aBsW+ffvwzDPP4LHHHlPk\n/EaDAUNDQ/B4fePuoje268k1eic8r9cLEaEL0nRTGMz2+Xzw+f0QBQF6SYLL5YJOFKETRQw5HVHf\n29vbm9A5J6tHrx9ZaOjsuSDLMQMtMS3S8rUBvL5Mk7Kg6O7uRklJCbKzs0P67evr61FVVYW1a9eG\nvH7hwpHdWk0mU6pKDBHoevL6vJDC/pV0EMZt9uf1+cZ1UYmCkPAYhcfjgV4UIQgCpNGgEL1eiE47\nnIPRb6+ajKmzXq8XEkaux3HogCzHbG9vl+U46UjL1wbw+jJNyoKisLAQzc3NKCsrCz528OBB2O12\ntLW1we1248CBA9ixYwe2b98OAFi3bh30en2qSgyhDwSFd3zXkxBhkNrr9Y67paoOibcoPB4PdKMH\nlKSRbjBRACQB8Hiid/8k45apPp8POp0IEYA37AZNRKRtKVtwZzAYYDAYQlYXd3R0oKamBgBQXV2N\n9vZ2PProowCAN954AydOnFDsPtEGvTTa9eQd1/U0dpA6MF4w5HZPMOideFBIoxvw6aSRWkRBgCSI\n8Hqjj1EkY+M+n883MkYiJn5NRKROgj/F+0LceeedaGpqgiiK+M1vfoOSkhIsXboUTU1NaG9vx+OP\nPx7zsZqampJYKRGRdlVXV8f82pRv4TG2/9xsNmNwtL99cHAQeXl5cR0rngslIqLEpHzWk9/vD3Y/\nlZeXB1sFu3btChm/ICKi9JCyoPB4PKipqUFXVxdqa2vR2dmJ+fPnw2g0oqqqCpIkYcGCBakqh4iI\nYpTyMQoiIlKXtF1wN5mJ1l9owRdffIGZM2diyZIlqK2tVbocWURaR/P000+jsrIS3//+91W/4jvS\n9eXl5WHJkiVYsmQJrFarwhVOzb59+7Bo0SJUVVVh3bp1AICnnnpKE59fpGvT0md3+PBhLFq0CIsX\nL8YjjzwCIP7PTpVBMXb9xdDQEA4ckGcBWDpZunQpmpub8c477yhdiizC19H09PSgpaUF7733Hm6+\n+Wa8+eabClc4NZHWCd10001obm5Gc3Nz3BM10s2VV16J3bt3o62tDRcuXEBbWxtaW1s18fmFX9tH\nH32Em2++WTOf3bXXXos9e/agtbUVQ0ND6OzsjPuzU2VQjF1/cdddd2lyFWVzczMWL16MhoYGpUuR\nhcFggNlsDv6+f/9+WCwWAF+toVGzwPWN7ck9evQoFi9ejJ///OcKViaP4uJiGAwGACMLQI8cOaKZ\nzy/82nQ6HY4cOaKZz27sWjSHw5HQ/z1VBoXVakVubi6AkSm2am8ahps1axY++eQT7N69G01NTfjo\no4+ULkl2Wv0Mx07//vTTT9Ha2gqr1YodO3YoWJV8urq6cPHiReTl5Wnu8wtc23XXXae5z+6tt97C\nTTfdBJPJhPz8/Lg/O1UGxVTXX6Q7vV6PrKwsiKKIb3/725oMCq1/hgCC17RixQpNfIb9/f348Y9/\njL/85S/Izc3V1Oc39toA7X12y5Ytw6FDh5CTk4Ps7GwMDAwAiP2zU2VQaH39hc1mC/68Z88ezJs3\nT8Fq5BXomiktLUVraysAbX2GgXVCDocjOKithc/Q6/XiwQcfxNNPP40ZM2Zo6vMLvzatfXZutzv4\nc25uLqz8SGFfAAADwElEQVRWK9ra2gDE/tmpMii0vv7ivffew4IFC3DHHXdg9uzZKC0tVbqkKRu7\njuZb3/oWTp48iaqqKlRWVuLDDz/Evffeq3SJUxK+Tuijjz5CaWkpLBYLTp8+jVWrVild4pRs3boV\n+/fvx09/+lMsWbIEn332mWY+v/Br6+rq0tRn984778BisQSv58EHH0RlZWVcnx3XURARUVSqbFEQ\nEVHqMCiIiCgqBgUREUXFoCAioqgYFEREFBWDgoiIomJQEBFRVAwKIiKKikFBGa21tRXr16+P+NwX\nX3yB3bt3T/jea665Bq+99pqs5yRKRwwKynhjd3wd6+TJk2hubo74XFdXF+6880689dZbsp5zIvv2\n7UNjYyNef/31hM5HNBUMCqJRjzzyCO666y4sW7YMVqsVmzdvxiuvvBK898lY27Ztw8MPPwyXy4Xh\n4WG0trairq4Oy5cvR2VlJRwOB9xuN1asWIG6ujrcd999ePnll6OeM7CjZziPx4O33noLdXV1UVs4\nRMnCoCACsGPHDsyZMwe7du3Cj370I2zatAkPP/wwfvCDH2Dnzp3jXn/w4EHcdtttWLp0afB5o9GI\nf/zjH6irq8OuXbvw5ptvYtGiRWhsbER+fv6k53z++ecj1tbY2AhBEPD8889DFPlfllJPUroAIqX5\n/X4cPXoUW7Zswb/+9S94PB6Ul5djov0yT5w4gUOHDqGurg5DQ0O45ppr8L3vfQ833ngjgJEbT1mt\nVnR3d+OWW24BANx6660xnTOS999/H/X19WhsbMQdd9wh45UTxYZBQRlPEARce+21eOihh7B27VoA\nI/co+OCDD+DxeMa9ftu2bXjxxRdx5513Ahi5uY3P5xs37jB37tzgtuNdXV1YuHBhyPORzun1etHb\n24vi4uLg62w2G4qLi9HR0YHf//73sl47USzYjqWM5vf7odPpsGzZMnz++eeorq7GXXfdhbfffhs3\n3ngj9uzZg/vuuy/kPY2NjaioqAj+fsMNN+D9998fd+wVK1Zgz549uPvuu3H+/Hno9frgc4IgRDzn\nyZMn8fjjj4cc5+6778bf/vY3PPDAA8jKypL5X4BocrwfBWW0V199FXa7HQ8//HBSju/1eqHT6fDI\nI4/goYcewu233x719du2bUNBQQEsFktS6iFKBIOCMtZrr72G5557Dn//+99RVFSUlHPU1tbCZrPh\nG9/4Bl566aWknIMo2RgUREQUFccoiIgoKgYFERFFxaAgIqKoGBRERBQVg4KIiKJiUBARUVQMCiIi\niur/A68HLeFiVp4AAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure()\n", + "\n", + "rateunits = u.gigaparsec**-3*u.year**-1\n", + "data = np.loadtxt('bns_highspin_pycbc_unlam_lnvt.txt')\n", + "rate, pdf = (data[:,0])*rateunits, 1e3*data[:,1]*rateunits**-1\n", + "bns_prior = BNSDistribution(rate, pdf)\n", + "o1scenario = Scenario(bns_prior)\n", + "\n", + "posterior = BeamingAnglePosterior(o1scenario, efficiency_prior=DeltaDistribution(1.0), grb_rate=3/u.gigaparsec**3 / u.year)\n", + "chain = posterior.sample_theta_posterior(nburnin=200, nsamp=200)\n", + "# Calculate the bandwidth of the poserior for the KDE\n", + "theta_bw = 1.06*np.std(posterior.theta_samples)*len(posterior.theta_samples/100)**(-0.2)\n", + "posterior.get_theta_pdf_kde(theta_bw)\n", + "theta_bin_size = 3.5*np.std(posterior.theta_samples) \\\n", + " / len(posterior.theta_samples/10)**(1./2)\n", + "theta_bins = np.arange(posterior.theta_range.min(), posterior.theta_range.max(),\n", + " theta_bin_size)\n", + "plt.hist(posterior.theta_samples, bins=50, normed=True, log=True, histtype='stepfilled',\n", + " )\n", + "plt.plot(posterior.theta_grid,posterior.theta_pdf_kde, \\\n", + " color='k')\n", + "plt.axvline(np.percentile(posterior.theta_samples, 95.0), linestyle='--')\n", + "plt.ylim([1e-3, 1])\n", + "plt.xlim([0,30])\n", + "plt.xlabel(r'Jet Angle, $\\theta$')\n", + "plt.ylabel(r'$p(\\theta | D,I)$')" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$[0,~18.362073,~36.030198,~\\dots, 98.03755,~96.80477,~95.585731] \\; \\mathrm{}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 39, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rate*(pdf*rateunits**-1)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/grbeams-0.1.0-py2.7.egg/grbeams/beamingangle.py:175: RuntimeWarning: divide by zero encountered in log\n", + " + np.log(self.comp_efficiency_prob(efficiency))\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAESCAYAAADjS5I+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYU3e+P/B3YgBFJSAIiChWFLG4UfcFjCIUsUptraMd\nr+2dzmiv7UwrnS5zb7Xa6fZr1bHeaev1drljaxdxqdWiVpaAIggoLSriRgUX3AUUMEiS3x+UyBqS\nkJOTnLxfz5Onmpzlc0z1zfd8lyPT6/V6EBERtUEudgFERGTfGBRERGQUg4KIiIxiUBARkVEMCiIi\nMopBQURERjEoiIjIKAYFEREZpRC7gLbk5eVh9+7d0Ol0eOONN8Quh4jIadm8RVFWVoaRI0fC3d0d\nOp3O8H5CQgIiIyOxdOlSAMCOHTuwbNky+Pj44MqVK7Yuk4iIfmPzoPD29kZqairGjRtneC8/Px9V\nVVXIyMhAbW0t8vLyDJ9xhREiInHZ/NaTq6srXF1dmwRAdnY2oqOjAQBRUVHIzs7GrFmz8Pe//x33\n7t2Dn5+frcskIqLf2EUfRXl5OYKDgwEASqUShYWFGD16NEaPHm10v5SUFFuUR0QkOVFRUSZvK1pQ\nyGQyw6+VSiUqKysBAJWVlfD09DT5OIOTE6EtuH+rqnrBcxj29GVUV7+CTp06Wa9gEbz00ktYvXq1\n2GUIRsrXJ+VrA3h9js7cH7JFGx6r1+sNt5/Gjx9vKDw5OblJ/4W5fJRK1Nb+DceOHbNKnUREzs7m\nQVFXV4fo6GgUFBQgNjYWubm5CA8Ph5ubGyIjI6FQKDBq1CiTj6dro7M7MzPTWiUTEUmCWq3GihUr\nzN7P5reeFAoF9u3b1+L9tWvXWnQ8uUwGbSvvHzhwAEuWLLHomPZi/PjxYpcgKClfn5SvDeD1OSqV\nSgWVSuU4t56ElpOTI3YJHTZhwgSxSxCUlK9PytcG8PqcjWSDoqysDBUVFWKXQUTk8Bw+KFrro3jj\nDWDo0KEoKCgQoSIiIvtkaR+FwweFvNEw2wYrVgDh4eHIz8+3fUFERHZKpVI5Z1C0hUFBRGQdDAoi\nIjLK4YOirXkUYWFhOHnyJOrq6mxcERGRfWIfRTPu7u7w9/fHr7/+auOKiIjsE/soGmn4cxg8eDBO\nnDghai1ERI5OkkGxcmX9fxkUREQdJ8mgaMCgICLqOIcPirY6swEGBRFRY+zMbkVoaChOnDjBx6kS\nEYGd2a3y9vZG586dcenSJbFLISJyWJIMijfeuP/rAQMG4OzZs+IVQ0Tk4CQZFI1bVsHBwQwKIqIO\nkGRQNMYWBRFRxzh8UBgb9QSwRUFE1ICjntrAoCAiqsdRT21gUBARdYwkg6JxYPbs2RO1tbUoLy8X\nrR4iIkcmyaBoWOsJAGQyGVsVREQdIMmgaI5BQURkOacIiv79+6O4uFjsMoiIHJLDB0V7w2MBICgo\nCKWlpTaohojIfnF4rBFBQUEoKSmxQTVERPaLw2MbabzWEwD07duXLQoiIgtJMiiaByaDgojIcpIM\niuY8PT2h0+lQUVEhdilERA7HKYJCJpOxn4KIyEJOERQAbz8REVmKQUFEREZJMihaG/3FuRRERJaR\nZFA0XuupQd++fdlHQURkAYXYBXSUKTOzAd56IiJSq9VQq9WIiIgwaz+Hb1GYMjMb4OxsIiLOzG5H\nr169cPXqVWi1WrFLISJyKE4TFC4uLvDy8sLVq1fFLoWIyKFIMiiar/XUoHfv3rh48aJtiyEicnCS\nDIq2bsH17t0bly5dsmktRESOTpJB0ZaAgAC2KIiIzORUQcEWBRGR+ZwqKNiiICIyn1MFBTuziYjM\nJ8mgaKszOyAggLeeiIjMJMmgaG2tJ4AtCiIiS0gyKNri7e2N6upq1NTUiF0KEZHDcPigMHVRQKD+\nSXe9evXi7Scickpqtdo513oydVHABrz9RETOiosCmogd2kRE5pFkULS11hPAFgURkbkkGRTGWlYM\nCiIi80gyKIzx9/fH5cuXxS6DiMhhOGVQXLlyRewyiIgchlMGBVsURESmc7qg8PPzY1AQEZlBkkFh\nrDPb29sblZWVqK2ttVk9RESOTJJB0dZaTwAgl8vh6+vLZ2cTEZlIkkHRHt5+IiIynVMGBTu0iYhM\n57RBwSGyRESmccqg4K0nIiLTSTIojK31BPDWExGROSQZFO2tostbT0REppNkULSHLQoiItPZbVAc\nPXoUs2fPxqlTp6x+bPZREBGZzmZBUVZWhpEjR8Ld3R06nc7wfkJCAiIjI7F06dIm2w8dOhSzZ88W\npBbeeiIiMp3NgsLb2xupqakYN26c4b38/HxUVVUhIyMDtbW1OHz4MHbt2oUdO3YAAPR6PfRmPBPb\nVEqlEhqNBtXV1VY/NhGR1NgsKFxdXaFUKpv8w5+dnY3o6GgAQFRUFLKysvDII48gPj4epaWl+Omn\nn/Dll1/i3r17Zp2rvc5smUwGPz8/tiqIiEygEPPk5eXlCA4OBlD/U35hYaHhs759+2LTpk3tHkOn\na9riqKurw8qVwKJFxp+L7e3tjWPHjsHNzc2Cym3j9u3bkn6+t5SvT8rXBvD6HM3BgweRlZVl+H1c\nXJxZ+9s8KGQymeHXSqUSlZWVAIDKykp4enqafTy5XAZto98rFPWXFBAQYHS/Pn36QKvVtrudmC5d\numTX9XWUlK9PytcG8PoczZw5czBnzhzD71NSUsza3+ajnhr3O4wfP95QcHJycpP+C6Fx5BMRkWls\nFhR1dXWIjo5GQUEBYmNjkZubi/DwcLi5uSEyMhIKhQKjRo0y+7g6Czu7OZeCiJyNWq3GivY6cVth\ns1tPCoUC+/bta/H+2rVrO3RcuazprSdT9ezZE0VFRR06NxGRI1GpVFCpVPZ/68kW2lvrCQB8fX1x\n7do14YshInJwkgwKU1pWPXv25FPuiIhMIOrwWGuwtI+CLQoicjZqtRpqtRoRERFm7efwLQp5o+G2\n5mCLgoicjUqlsqgz2+GDwlLe3t4oLy+HVmtJVzgRkfNw2qBQKBRQKpW4ceOG2KUQEdk1SQaFqS0r\n9lMQEbXP4YOitc7slStN29fX15f9FETkNCydcOfwQWFpZzZQ36HNFgUROQtLO7Mdfnhsc3p9/UOR\n7pbfMrzXyc0NLl3cW2zLFgURUfskFxQ1//dPAHNw49U/Gd7r+vx/wXNoeIttOUSWiKh9kgsKVN8B\nAOgunLv/XqNHrzbm6+uLo0eP2qAoIiLH5fBB0Vpn9osDNpi0L1sURORMODO7kYQQ04KCw2OJyJlw\nZrYF2KIgImqfUwcFWxRERO1z6qDo0aMHKisrce/ePbFLISKyWw4fFJYuMw4AcrkcPXr04HpPROQU\nODO7kTWnFpm8PyfdEZGzEHxmdkVFBfLz83HlyhV0794dI0eOhJ+fn9kntIW1ZxaZPPKJy3gQERnX\nblDcvHkTe/bsgYuLCwYMGIA+ffrgzp07UKvVqK6uxuTJk9G/f39b1CoItiiIiIwzKSiefPLJFu8P\nHz4cAFBUVGT9qmyIQ2SJiIxrt49iwIAB+Oabb9r8PDQ01KoF2RqHyBIRGWdSH8WaNWvw1VdfoW/f\nvpgwYQImTJiAmzdvYvTo0ULXJ7iePXviyJEjYpdBRGS3TAqKxMRE9OvXDxcuXMDBgwfx3//930hJ\nSbHbBfVMXesJYIuCiKg9JgVFv379AACBgYGYO3cu5s6di+zsbCHrMllr8yhMHfEEsDObiJyHzRcF\nHDdunKW7WlVHnnAHcHgsETkPmy8K+MEHHyA3NxfHjx+39BB2gS0KIiLjLA6KiIgIjB49Gr169bJm\nPTbn6emJ6upqaDQasUshIrJLJvVRaLVabNmyBefPn0dERATGjh1ruPXUo0cPQQsUmkwmg4+PD65f\nv47evXuLXQ4Rkd0xqUWxbNkyuLm5ITg4GImJiVi0aBFqa2uFrs1i5qz1BHDSHRGRMSa1KEaMGIFH\nH30UADB79mzcvHkT//znP5GQkCBocZYyZ60ngENkiYiMMalF0bz10KNHD/j6+gpSkBjYoiAiaptJ\nLYoNGzZg+/bt8PT0xLhx4zBhwgR069ZN6Npshi0KIqK2mRQUmzdvhr+/PyorK5GdnY3ExET8+OOP\nhttRYurIg4sacIgsETkDSyfcmRQU/v7+AAAPDw/ExMQgJiYG8fHx5lcpALlMBm0Hj9GzZ08UFxdb\npR4iInulUqmgUqmQkpJi1n4mtyganistazQTurCwEK6urnjiiSfMOqnQzFnrCWCLgojIGJOCYu7c\nuULXYVXmjHgC2JlNRGSMSUGxZcsWpKWlQS6XY9KkSZg1axa6dOkidG02w85sIqK2mTQ81t3dHR99\n9BEefPBBVFVV4YUXXkB6errQtdkMWxRERG0zKSh8fHywbNkylJWVIS4uDuvXr8fFixeFrs1mPDw8\nUFtbi5qaGrFLISKyOybdehozZgweeughJCUlYc2aNbhz5w6mT58udG02I5PJDLef+vbtK3Y5RER2\npd0WRcMDihQKBWbNmoX3338fH3/8MWbOnAkAyM3NFbZCC5i71hPA51IQEbWl3RZFWFgYPvvsM9y9\nexd9+vRB9+7dUVVVhdLSUri7u+ORRx6xRZ1mMXetJ4BDZImI2tJuUHTv3h3PPPMMdDodTp48ifLy\ncvj6+iI2NhYKhUl3rhwCWxRERK0z+V96uVwOtVqNjIwMeHt74+2334ZSqRSyNptii4KIqHVmPeHu\nwQcfxDfffIO3334bGzduFKomUXCILBFR68y6d5SUlIRdu3bhgQcegEajQV1dHS5evIigoCCh6muX\nNRYFBOpbFCdPnrTKsYiI7JGgiwI2iI+PxwMPPACtVotDhw5h5cqVyM7Oxr59+8w6qTW1tiiguWs9\nAWxREJH0CbooYIMJEyYYfh0YGIjHH3/cLv9xNXfEE8BlPIiI2mJWH0VrpPKkO3ZmExG1rsNBIRUc\nHktE1DoGxW+6du0KvV6PqqoqsUshIrIrDIrfyGQydmgTEbVCkkFhyVpPADu0iYhaI8mgWHvGsqBg\ni4KIqCVJBoWl2KIgImqJQdEIh8gSEbXEoGiEQ2SJiFpiUDTCFgURUUuSDApL1noC2JlNRNQaSQaF\nJWs9AezMJiJqjSSDwlJsURARtcSgaKShM1tvpWdcEBFJAYOiEXd3dygUCty+fVvsUoiI7IZZz6Nw\nVHqZDNVlF5u85+bdE51cXVts29BP4eHhYavyiIjsmt0GRUpKCnJzc3H58mWsXbvWrH3XnFrUpEP7\nzpsvosrlfijIevqjx/I16OLT8lkaDUNkg4ODLS+eiEhCbHbrqaysDCNHjoS7uzt0Op3h/YSEBERG\nRmLp0qVNto+KisJrr73WZFtTtVjrqaYa+sry+6/bFW3uyw5tIqKmbBYU3t7eSE1Nxbhx4wzv5efn\no6qqChkZGaitrcXhw4exa9cu7NixAwCwdu1aLFiwwFYlAuAQWSKi5mx268nV1RWurq5NRhRlZ2cj\nOjoaQH0LIisrC88//zwAYMOGDcjLy4NcLseYMWPaPK5OZ9kIpYqKCtyqrWvxfpcuXXDmzBlcunTJ\nouNa0+3bt+2iDqFI+fqkfG0Ar8/RHDx4EFlZWYbfx8XFmbW/qH0U5eXlhr4ApVKJwsJCw2eLFi3C\nokXtLxcul8ugteDcSqWy1T6K/v374/z58wgICLDgqNZ16dIlu6hDKFK+PilfG8DrczRz5szBnDlz\nDL9PSUkxa3+bD4+VyWSGXyuVSlRWVgIAKisr4enpaetyWmAfBRFRUzYPCr1eb7j9NH78eEOyJScn\nN+m/6AhL13oC2EdBRNSczYKirq4O0dHRKCgoQGxsLHJzcxEeHg43NzdERkZCoVBg1KhRZh9X18os\nakvXegLqg+LKlSsW709EZK/UajVWrFhh9n4266NQKBTYt29fi/fNnSPRnFxmWR9FW/z8/BgURCRJ\nKpUKKpXK/vso7J2vry9u3LiBurqWI6KIiJwRg6IZhUIBHx8ftiqIiH7j8EHRWh9FR/Xq1QtlZWVW\nPy4RkZgs7aNw+KCQNxpu22DNqfbnXxjDoCAiKVKpVM4ZFK1psdaTmRgURET3STIoOqpXr16Smr5P\nRNQRDh8UQvRRBAQEsEVBRJJj9/MohGLteRRAfYti9+7dVj4q6XQ63C46Bn2tBgDQSaNBZXUlugcP\narK0CxEJw9J5FA4fFEJgH4UwdDodanZ8g3sH7k+8lP9xKboHDxKxKiJqj8PfempNR9Z6AhgURESN\nSTIoOrLWEwD4+/vjypUr0GqtfVOLiMjxOHxQCNGZ7erqCqVSievXr1v92EREYuGEOyvjyCcikhpO\nuLMy9lMQEdVjULSBQUFEVE+Sw2PXnFrU4Q5tzs4mKau6dAGaooIm77n0G4ju/QeKVBHZM0kGxdoz\nHQ+KwMBAHD161EoVUWM3a+uQcrUcPd1c8LAAgxGofbq7NahavbzJe93fWS9SNWTvHP7WkxCjngCg\nb9++KC0tFeTYzuxw6UVE7z+O1KuVeLfoIp75YB0fEkVkIxz1ZGV9+/bF+fPnBTm2szp//jz+sGkb\n3hsShE8e6o+dE0NRdfculi9f3v7ORNRhHPVkZWxRWN/y5csxf+RwRPt5AgBc5XKs+/NifPrppygq\nKhK5OiJqC4OiDZ6entDpdKioqBC7FEk4efIk9uzZgyWRY5q87+fliRdffBHvvfeeSJURUXskGRQd\nXesJAGQyGfr06cNWhZW8+eabeOGFF+DRuXOLz5YsWYKdO3fi4sWLIlRGRO2RZFB0dMRTA/ZTWEdh\nYSGSk5Px3HPPtfq5p6cnHn/8cXz55Zc2royITCHJoLAW9lO0TltXh8ozJ1F5psjwunP+XJvbr1y5\nEi+99BK6d+/e5jZPPfUU/vWvf0HP4bJEdsfh51EINTwWYFC0RXvvHm7/823oThca3nN/YTm69enX\nYtuCggKkp6fj888/N3rMCRMmQKPR4Oeff0Z4eLi1SyYi1A+PVavViIiIMGs/h29RCDU8FgD7KKxg\nxYoVeOWVV9C1a1ej28lkMsTHx2Pnzp02qozI+XB4rADYR9ExR44cQXZ2Np599lmTtp85cyZ27dol\ncFVEZC5JBsWaU4uschzeeuqY5cuX47XXXoO7u7tJ20dEROD06dNcjJHIzkgyKNaesU5QBAYG4tKl\nS1xiwgIZGRk4fvw4Fi9ebPI+Li4umDZtGvbt29f+xkRkMw7fmW0J/Z3bqDn+C2oadYTLfHzhGTYc\nskZ9Hm5ubvD390dpaSn69+8vRqkOSafT4eWXX8Zbb70FNzc3s/adMmUK0tLSsHDhQoGqI2dVc+0K\nak4cBXQ6w3uKwCB4DBgkYlWOwSmDAtV3UPX+35q81fnJxdA/OKxJUADAgAEDcObMGQaFGdatWwdX\nV1fMnz/f7H1VKhU++OADAaoiZ6fT3EXV2hWA5q7hvW5vrAXAoGiPJG89WVNDUJBpDhw4gLfffhtf\nfPEF5HLz//caPHgwampqcO7cOesXR0QWYVC0g0Fhuj179uCxxx7DV199hQEDBlh0DJlMBpVKhbS0\nNCtXR0SWcvigaG3CnTXWemoQHByMs2fPWu14UqTT67Hm68145plnsG3bNjz88MMdOl5ERAQOHjxo\nperImOI7d1FYWc0Z8U7C0udROHwfhVwmg7bZe9Za6wlgi6I91XVaLC04h2teVcjNzUVAQECHjzl2\n7FisX8+nrQlJr9fjnaIL2HrhBjp3kmOYsis+v3dP7LJIYCqVCiqVCikpKWbt5/AtCqEFBwejuLgY\nukYjJaherU6HPxw+Cze5HN+//5ZVQgIAhg0bhuLiYty+fdsqx6OWvtqyFalXK5AyOQxpk8NQo9Xh\nzXX/FLssslMMinZ07doVXl5euHTpktil2J23TlxAN4Uc/xjeD51dXa12XFdXV4wYMQJ5eXlWOybd\nV1FRgZWr1uDDEQ/A00UBV7kca4b3Q2LSHvzyyy9il+f0tFotbuZl4cZPPxhet47/IurtQQaFCUJC\nQnDixAmxy7ArORevYvflcqwa1g+dBFhva+zYsTh06JDVj0vAp59+isjx4xDmcX/GfA9XBf7y9EL8\n/e9/F7EyalCbugvVH75peNWdPCZqPQwKE4SFheH48eNil2E39Ho93sk8gv8M7Q1PF2G6ucaOHYvs\n7GxBju3MtFotPv74Y/zH00+1+Ozf5zyO9PR0FBcXi1AZ2TNJBoW11npqMGTIEAZFI+r0dNy6q8Gs\ngB6CnYMtCmFkZGTAw8MDo0cMb/GZe5fOWLBgAT777DMRKiN7JsmgsNZaTw3CwsJw7Ji4TT978tHH\nH2PxQw8KcsupQVBQEGpra3H58mXBzuGMtm7diieeeKLNz//4xz/iiy++gFbbfCwhOTNJBoW1hYWF\nobCwkGPNAVy+fBn7DxzArJB+gp5HJpNhxIgR7Fy1Ip1Oh+3bt+Pxxx9vc5uwsDD4+/tj//79NqyM\n7B2DwgTe3t5wd3fHhQsXxC5FdF999RXiZ85EN1cXwc81YsQI/Pzzz4Kfx1nk5OTAy8sLgwYZX9to\nzpw52Lp1q42qIkfAoDDRkCFDcPToUbHLEN13331n0WJ/lmBQWNe+ffsQGxvb7naPP/44tm3bxrlD\nZMCgaESv10On0xlejW81jRw50unH9V+4cAHFxcWImDTJJucbPnw4g8KKUlNTERUV1e52gwYNgpeX\nl6QHE2i0DEFzOPwSHq2xZK0nTfpe3NDUNHnPdZwKng8OAwCMGTMGn3/+uVXqc1Q//PAD4uLioFDY\n5n+b0NBQlJSUoKqqqt1nbpNxNTU1yMvLQ0REhEnbP/LII0hKSsL48eMFrsy2Ku/V4YW8M1Bfq8AY\nr27YWF4OL7GLcgAO36JobVFAS9Z60l88B83WjU1e+tsVhs/HjBmDnJwcp+7Q3rFjB+Lj4212PldX\nV4SGhnLEmRVkZmZi+PDh6Natm0nbT58+Hbt37xa4KtvS6/VIOHwaPV0VKIoJR5iHO/607E2n+jtt\n6aKADh8UcgGHaDbWu3dvKBQKlJSU2OR89kaj0SAzMxPR0dE2PS/7KawjNTUVU6ZMMXn7CRMm4MyZ\nM7h69aqAVdnW3tQ0FN+uwdtDguDWSY7/DA3EtZs38f3334tdmkF6ejoSEhKwLu0gbtZa/xHMKpXK\nOYPCVmQyGcaPH4/MzEyxSxFFbm4uQkNDoVQqbXpeBoV1ZGdnY5IZfUsuLi6IiorC3r17BazKtt79\nx4d4NSwILvL6Hy4VchlefuZprFq1StzCUD90OSEhAU8//TT8/Pxw7votzDhwAiXVGrFLA8CgMMu0\nadOwb98+scsQRUZGBiIjI21+Xs6K7zidTofDhw9j5MiRZu0npdtPubm5uH7zBqb1arqaQFzkJFy4\ncEH0H0ZWrlyJgwcP4siRI3j11Vex5okZ+OMDvlh8+Cxq7WD0GYPCDDExMfjpp5+c6p5mg/T0dEye\nPNmsffR6fZOXJfsNHjwYx48f51DNDjh9+jS8vb3h4+Nj1n6xsbH46aefJDFLe8OGDfjD759ssZqA\nQqHAvHnzsHnzZpEqA1JSUvD5559jx44d8PK637X+h36+8O3sgo0l10SrrYEkg8Laaz01CA4Ohru7\nu+g/fdhaXV0dsrKyzLp1AQDlJ47i+o5vDK9buzaj7lT7rYOKk4WGfTplp0J/rxYlhWxVWCovLw+j\nRo0ye7/AwED06tULOTk5AlRlO/fu3cP27dvxu9mPtvr53LlzkZiYKMoPgDU1NVi8eDHWr18PPz+/\nJp/JZDK8HhqIj89eRo2m1ua1NSbJ4bFrzyyy6lPuGvvd736HTZs2ITw8XJDj26MjR46gX79+8Pb2\nNms//eWLuPu/a8w+n/7mtSb7DXDrhKKiE3hgyFCzj2WP7t64jupjRwDt/c5KReAD8AgZLMj5cnNz\nMXr0aIv2jYuLw+7dux16mKxarUZwcDACAwJws5XPH3roIdTW1uLkyZMIDQ21aW2rV69GeHg4ZsyY\n0ernId27YLhnV2zJyMTSeU/btLbGJNmiENLChQuxadMm3L17V+xSbEas/okGIcquKDp1WrTzW5uu\nVoPq9e+javVyw0t7/Ypg57O0RQHcDwpHtnXrVqPrW8lkMkRFRZn9eNCOunXrFtauXYt33nnH6Hb/\nHuSLjftSbVRV6xgUZho0aBDGjBmDjz/+uM1t7t68jjslxU1emopbNqzSuizpn7CmEGU3nDh1SrTz\nO7K6ujr8/PPPZndkN2gYJnvlinBBJiS9Xo+kpCTMmjXL6HbTpk1DcnKy4PVoym8Z/k14Z9l/Yca0\nKPT1Nj7lb6JPd1y5VY6ioiLB62sLg8IC7733Ht59912kpaVBp9MhMzMTzz//PPr27Qs3NzeooqZh\n97w43HpunuGlueaYf9G0Wi0OHDhg8oxeIQz06IoiBoVFTpw4gcDAQHh4eFi0f8Mw2T179li5Mts4\nffo09Hp9uwshTp06FWq1WvCOe83VMtx6bh5KFz2B/92wAYuvFuFumfHFRjvJZHh04jhs2rRJ0NqM\nYVBYYPDgwdi0aRMWLlyIzp07Y9GiRejVqxf27t2LmzdvYsmCJ/Fc/ll8f+E6oNfVvxzU0aNH4efn\nB39/f9FqCFF2xYmTDApLdKR/okFcXBySkpKsVJFtJScnY9q0aZC1MzHX398ffn5+thmKrdch8cI1\nTPTujj5dTFuFOX7CWFEnBkqyM9uStZ7MFRMTg9LSUty9exddunRp8ll8zDQE7PgX5h06hQe6umG4\np+OuUyR2/wQA+HZ2RZ1Wi2vXrqFnz56i1uJoOtI/0SA2NhZ//etfUVdXZ7N1vqwlOTnZaP9EYw1P\nVRw2bJigNen0evzfuatYNayfyfuED+iPq1evoqSkBEFBQcIV1wZJtiiEGvHUnEwmaxESDQZ174KV\nD/bBXwvO2cWEGUtlZGSI2j8B1P85Dw4J4cQ7C1ijRREQEICgoCCHe4a5VqtFWlqaSSvmArZ7/G7a\ntUp0U3TCaC/Tf4CUy+WYPn06fvzxRwErM3J+Uc7qJGb28kLvLm749FfHXC9Hr9fbRYsCAEJDBjIo\nzKTRaFBoEj4EAAAMwElEQVRYWIgRI0Z0+FgzZszAjh07rFCV7RQUFMDf39/k26a2CorN569jQV+f\ndm+HNRcXF8egkCKZTIZlgwOx4dcruFVR0f4OdqaoqAjdunVDnz59xC4FgwYOFHXUhyM6evQoBgwY\nAHd39w4fa968efjuu+8caob8oUOHMG7cOJO3HzZsGH799Vfcvn1bsJpu3LyFzBu38UizpURMERMT\ng4yMDGg0tl//iUEhsOBunRHj54kP//czi/a/c/4cKk8eb/Kqu1vT/o5WkJ6ebhetCQAYGNwfpzjy\nySzW6J9oMGTIECiVSodaFPPQoUMYO3asydu7uLhg6NChyM/PF6ymLT8mYUpPD3i4dDJ7X09PTwwa\nNEiUmfIMChtYOqAX/u+7zSgrKzN739qzJ1GR8JThVblqGepqqgWosiW1Wi16/0SDgf374+TJk6Kd\nv6amxuEmWVqjf6KxJ598UtQhmuYyNyiA+lZFQUGBQBUBX2/fgScCzVtzq7EpU6ZArVZbryATSTIo\nhFrryVK9urji94/NxltvvSV2KSbT6/VQq9VmPcNASEF9+uDy5cuoqbFNa6qBXq/HW2+9BT8/P/j5\n+WHdunVWP8cNzT3864dd2Lhxo1Vve1izRQHUr0qwefNmlJeXW+2YQqmoqEBpaSmGDjVv2Zfhw4cL\nFhSFhYW4cu0aJvp0t/gYU6ZMQVpamhWrMo0kg2LtGfsKCgBIWLwI3377LYqLi8UuxSRFRUXo0qUL\n+vXrJ3YpAOpX+ezfvz9On7bNUh6aigrcvXEda9//f/j266/x88GDOJyXh3Xr1uHrr7+22nmyb9xG\n9P5CZB75Gdu3b0dISAhSUzu+XEN1dTXOnDlj9j+UxvTu3RvTp0/Hp59+arVjCiU3Nxfh4eFmD+cV\nskWRmJiI2dMfbrGCrTkmTZqEnJwcm7du7TYocnJysHr1arz44otil2IV3j288Oc//9msp0sdPHhQ\nuILakZaWZjetiQaDBg2yWT9F1S+5+GXJfLz5xnKsDwtA5y2fok8vfyQmJmLp0qXt/lRtyndXUq3B\nf+QXY92IB7BhxevYvn07Nm3ahHnz5nV43aH8/Hw8+OCDcHNz69Bxmnv55ZexevVqu38uiyW3nQBg\n6NChOHbsGA4cOGD1mrZu3YpZMR17QqSHhwfCwsJsMjqrMZsFRVlZGUaOHAl3d/cmIycSEhIQGRmJ\npUuXNtl+zJgxAIDOnTvbqkTBJSQkYO/evSY/AzorK8vic1VdKEXFiaNNXrWVpo+8SktLg0qlsvj8\nQggJCbFdP0WtBqsy8/FUHx8EVZdDf71+iHN4eDhmzZqF1atXG93dlO/u9WOleLa/Hyb53F9eY+rU\nqUhMTMT8+fM7FIq5ubmGv0PWNGLECMyYMQP/+Mc/rH5sa7I0KJRKJXx8fKy+ZMnJkydx/fp1jH2o\n46tOi3H7yWZB4e3tjdTU1CbD1fLz81FVVYWMjAzU1tbi8OHD2LVrl2G89ksvvQQXF9OmuDsCDw8P\nvPHGG/jTn/6EurrWn4d7584dfPfdd3jzzTdx/PhxVFZVWXSu2oslqPzrv99/Lf8z7t2uNGlfrVYL\ntVqNqVOnWnRuoQwaNMhmQVF++w5+ulKOZ/r5tvjslVdewf/8z/+gutryQQXpZTdQWq3BH/r5tfhs\n8uTJWLlyJebOnWvxLYacnBxBggIAVq1ahZKSEqMLY4pJr9dbHBRA/e2na9es+7CgrVu34rHHHoNc\n3vF/clUqlc07tG0WFK6urlAqlU0eDpKdnY3o6PqmWFRUFLKysvDII48gPj4e27dvx6pVq9Cpk/nD\nyOzZs88+i65du+LVV19t8mdRXV2NDz74AMHBwdi4cSM0Gg3Onj2LcQv/iP3XTfsHvjU6vR4Z1yrx\nrzMXkbZ/v0mLnmVlZaF3794IDAy0+LxCsGVQbE1LR2RPD3i6trzHPXDgQIwdOxbffvutxcf/5EQJ\nXhzYy/D85uaeffZZhISE4LXXXrPo+EIGhaenJ2bPno01a9ZgwYIFyM/Pt+ihP1XnS3Dr0P4mr+pL\nFztcX2lpKeRyucXzf8LCwnDzZmtPrrDcli1bTF5KpD0TJ05EXl6eTfspZHobP9ZpypQpSElJgVwu\nx7vvvouRI0ciJiYGKSkpyMrKwuuvv27ysWy9fjwRkVSYurQJIMKigI2nrSuVSlRW1v+0XFlZCU9P\nT7OOZc6FEhGRZWw+6kmv1xuaqePHjze0CpKTk82abk9ERLZhs6Coq6tDdHQ0CgoKEBsbaxjn7Obm\nhsjISCgUCqtODiIiIuuweR8FERE5FrudcNeetuZfSEFJSQn8/f0xdepUxMbGil2OVbQ2j2bVqlWI\niIjAv/3bvwn+CEqhtXZ9np6emDp1KqZOneoQy14Yk5OTg4kTJyIyMhIvvfQSAOCDDz6QxPfX2rVJ\n6bs7fvw4Jk6ciMmTJ2PJkiUAzP/uHDIoGs+/0Gg0OHz4sNglWV1MTAxSU1Md9lnFzTWfR3Pt2jWo\n1Wrs378fw4YNE/Uxj9bQ2jyhoUOHIjU1FampqWYP1LA3/fr1Q1paGjIyMnD16lVkZGQgPT1dEt9f\n82s7duwYhg0bJpnvLjQ0FJmZmUhPT4dGo0Fubq7Z351DBkXj+RfTpk3r0Axme5WamorJkydj7dq1\nYpdiFQ3zaBrk5eUZZn43zKFxZK3NEzpx4gQmT56Mv/3tbyJWZh2+vr5wdXUFUL/uVmFhoWS+v+bX\n1qlTJxQWFkrmu2s8F626utqiv3sOGRTl5eXw8Khf9kCpVDp807C5gIAAnD59GmlpaUhJSTF5yQ9H\nItXvsPHw7zNnziA9PR3l5eXYtWuXiFVZT0FBAa5fvw5PT0/JfX8N1zZ48GDJfXc7d+7E0KFD0blz\nZ3h5eZn93TlkUHR0/oW9c3FxQZcuXSCXyzFjxgxJBoXUv0MAhmuKj4+XxHd469Yt/OUvf8Hnn38O\nDw8PSX1/ja8NkN53N3PmTBw9ehTdu3eHu7s7Kn574qap351DBoXU51/cuXPH8OvMzEwEBweLWI11\nNdyaGT16NNLT0wFI6ztsmCdUXV1t6NSWwneo1WqxYMECrFq1Cj179pTU99f82qT23dXW1hp+7eHh\ngfLycmRkZAAw/btzyKCQ+vyL/fv3Y9SoUZg0aRICAwOt+pQysTSeR/Pwww/j3LlziIyMREREBH75\n5Rc8+uijYpfYIc3nCR07dgyjR4+GSqXChQsXMGfOHLFL7JDExETk5eXhlVdewdSpU1FcXCyZ76/5\ntRUUFEjqu9uzZw9UKpXhehYsWICIiAizvjvOoyAiIqMcskVBRES2w6AgIiKjGBRERGQUg4KIiIxi\nUBARkVEMCiIiMopBQURERjEoiIjIKAYFObX09HQsW7as1c9KSkqQlpbW5r4hISHYvHmzVc9JZI8Y\nFOT0Gq/42ti5c+eQmpra6mcFBQWYMmUKdu7cadVztiUnJwdJSUnYsmWLRecj6ggGBdFvlixZgmnT\npmHmzJkoLy/Hhg0b8OWXXxqefdLYtm3bsHjxYty9exf37t1Deno64uLiMGvWLERERKC6uhq1tbWI\nj49HXFwc5s+fj40bNxo9Z8OKns3V1dVh586diIuLM9rCIRIKg4IIwK5duxAUFITk5GQ899xzWL9+\nPRYvXoyFCxdi3759LbbPz8/HQw89hJiYGMPnbm5u+OGHHxAXF4fk5GR8//33mDhxIpKSkuDl5dXu\nOT/55JNWa0tKSoJMJsMnn3wCuZx/Zcn2FGIXQCQ2vV6PEydO4JtvvsHevXtRV1eH8ePHo631Ms+e\nPYujR48iLi4OGo0GISEhmDdvHoYMGQKg/sFT5eXlKCsrw/DhwwEAI0aMMOmcrTlw4AASEhKQlJSE\nSZMmWfHKiUzDoCCnJ5PJEBoaiqeeegpLly4FUP+MgkOHDqGurq7F9tu2bcNnn32GKVOmAKh/uI1O\np2vR79C/f3/DsuMFBQUYM2ZMk89bO6dWq8WNGzfg6+tr2O7OnTvw9fVFdnY2PvzwQ6teO5Ep2I4l\np6bX69GpUyfMnDkTv/76K6KiojBt2jTs3r0bQ4YMQWZmJubPn99kn6SkJEyYMMHw+7CwMBw4cKDF\nsePj45GZmYnp06fjypUrcHFxMXwmk8laPee5c+fw+uuvNznO9OnT8e233+L3v/89unTpYuU/AaL2\n8XkU5NS++uorVFVVYfHixYIcX6vVolOnTliyZAmeeuopjB071uj227ZtQ48ePaBSqQSph8gSDApy\nWps3b8ZHH32ErVu3wsfHR5BzxMbG4s6dOxg4cCC++OILQc5BJDQGBRERGcU+CiIiMopBQURERjEo\niIjIKAYFEREZxaAgIiKjGBRERGQUg4KIiIz6/27aT7spFsUCAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "posterior = BeamingAnglePosterior(o1scenario, efficiency_prior=UniformDistribution(), grb_rate=3/u.gigaparsec**3 / u.year)\n", + "chain2 = posterior.sample_theta_posterior(nburnin=1000, nsamp=400)\n", + "\n", + "plt.figure()\n", + "theta_bw = 1.06*np.std(posterior.theta_samples)*len(posterior.theta_samples)**(-0.2)\n", + "posterior.get_theta_pdf_kde(theta_bw)\n", + "theta_bin_size = 3.5*np.std(posterior.theta_samples) \\\n", + " / len(posterior.theta_samples)**(1./3)\n", + "theta_bins = np.arange(posterior.theta_range.min(), posterior.theta_range.max(),\n", + " theta_bin_size)\n", + "plt.plot(posterior.theta_grid,posterior.theta_pdf_kde, \\\n", + " color='k')\n", + "plt.hist(posterior.theta_samples, bins=theta_bins, normed=True, log=True, histtype='stepfilled',\n", + " )\n", + "plt.axvline(np.percentile(posterior.theta_samples, 95.0), linestyle='--')\n", + "plt.ylim([1e-3, 1])\n", + "plt.xlim([0,30])\n", + "plt.xlabel(r'Jet Angle, $\\theta$')\n", + "plt.ylabel(r'$p(\\theta | D,I)$')" + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 81, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAESCAYAAADjS5I+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwW9XZLvBHW1uSJVuyk5iQQoCQhAy3ACFxyAU7IsGp\nMZBAm3JKC6Uz9ISZlCkkmenXzkemdHp62imXphe+8DFDOwfooUPCJZAa2sSOY0hscFIfnAuUcklC\nEmMndmRZkq3Lls4f0pYlWZIlW9LW3np+M53Gti7vjgY/WWvt9S5dKBQKgYiIKAVB6QKIiKi4MSiI\niCgtBgUREaXFoCAiorQYFERElBaDgoiI0mJQEBFRWgwKIiJKS1S6gFQOHjyIt99+G8FgED/72c+U\nLoeIqGQVfETR09ODhQsXwmKxIBgMRr+/adMm1NXVYePGjQCAnTt3YsuWLaiurkZvb2+hyyQiooiC\nB8W0adPQ0tKCJUuWRL/X1dUFt9uNtrY2+Hw+HDx4MPozdhghIlJWwaeejEYjjEZjXAB0dHSgvr4e\nALBq1Sp0dHRgzZo1+MUvfgG/348LL7yw0GUSEVFEUaxROBwOzJkzBwBQWVmJY8eOoaamBjU1NWmf\n19zcXIjyiIg0Z9WqVRk/VrGg0Ol00T9XVlbC6XQCAJxOJ6qqqjJ+nav2bMcBlx/PHPkce/fuBQA8\n/nj4f2q3efNmPPXUU0qXkTdavj4tXxvA61O7bP+RrdjtsaFQKDr9tHTp0mjhe/bsiVu/yITJYIDX\n641+/fOf565OIqJSV/CgCAQCqK+vR3d3NxoaGtDZ2YkFCxbAZDKhrq4Ooihi0aJFGb9eMBSC8OXn\nGHa78lg1EZH6tba24vEJTLcUfOpJFEXs3r17zPe3bt06odcTdDqYPC54vb7JllZ0li5dqnQJeaXl\n69PytQG8PrWy2+2w2+3qmXrKJZMgxE09acWyZcuULiGvtHx9Wr42gNdXajQRFGWiHp7hYaXLICLS\npKK4PXYygqEQzKIewzFBwY4fRERjtba2orW1FbW1tVk9T/UjCkGng1kfP6LQwq2xRES5ZrfbJ7SY\nrfqgAADD8BCCwSD8fr/SpRARaY4mgkIXCMBSVgaPx6N0KUREmqP6oAhGNu2Zy0wMCiKiNCa6j0L1\nQSFEWoFYzGYGBRFRGiW9RgEA5pipJy5mExHljmaCwlJWBrfbDYC9noiIckkzQWHmYjYRUV5oYsMd\nwDUKIqLxlPSGOyB+6omIiMbiYjZvjyUiygsNBcXoiIK9noiIckczQWG1WTE0NASAt8cSEeWSdoJC\nh+i520RElDuqDwr5rier0YjBwUGFqyEiKl4l38LDajFzREFElEbJ3/VUUW5hUBAR5YFmgsIqitGg\n4GI2EVHuaCYozKc+h8PhAMBeT0REuaSZoCj3OOEcdChdBhGR5mgnKM72wOkcUroMIiLN0UxTQKuo\nh3NoCKFQCIBO2aKIiIpQyTcFNOkF6HQ6eL1ehSsiIipOJX97LABYKyrgdDrZ64mIKIc0FRQ2qxVO\np5O3xxIR5ZB2gsLvg7XMxDYeREQ5pp2gAGAtL+fubCKiHNNYULCNBxFRrmkqKGwVVgYFEVGOaSoo\nKsrCrca5mE1ElDvaCgpBgNPpZK8nIqIc0szObACwcY2CiCilkt+ZDfDwIiKidLgzG4DVbOY+CiKi\nHNNWUFg49URElGvaCorIGgV7PRER5Y62gsJsZq8nIqIc01ZQWLhGQUSUaxoLCq5REBHlmqaCotwz\nCKfTGTnljoiIckFTQSF89CEEQeApd0REOaSpoAAAm82Gxx4LKF0GEZFmaCooQoMO2Coq8NRTFUqX\nQkSkGdoKit7TqLBYlC6DiEhTNNUUEABsFeUKVUJEVNzYFDDCWsFpJyKiZNgUMMJmtSpdAhGRpmgu\nKCqqL0BDw/tKl0FEpBmaCwqr1Yrly3crXQYRkWaofjE7UfmZExjU65Uug4hIMzQ3oqi2VeDs2bNK\nl0FEpBkaDAob+vr6lC6DiEgzNBgUVo4oiIhySJNB8ckn31G6DCIizdDcYna1zQqn82GlyyAi0gzN\njSjKDOHs83q9PJeCiCgHNBcUgQN7AQDnenrg6u1RuBoiIvXTXFAgKAEAXC4XRxRERDmgvaCIcLvd\nSpdARKQJmgyKGTP+G0Mul9JlEBFpgiaDYt68/ws3g4KIKCc0GRQVFRVwuoaULoOISBOKNigOHz6M\nu+++G5988klWzwv298Ec8MHt8eSpMiKi0lKwoOjp6cHChQthsVgQDAaj39+0aRPq6uqwcePGuMfP\nnz8fd999d/ZvNOyBxe2E28XFbCKiXChYUEybNg0tLS1YsmRJ9HtdXV1wu91oa2uDz+fDoUOHsGvX\nLuzcuRMAEAqFJnSLq8VogNvDoCAiyoWCBYXRaERlZWXcL/6Ojg7U19cDAFatWoX29nbccccdWLt2\nLU6ePIl//OMfePHFF+H3+7N6rw9P/k/817P/jQc3/DCn10BEVIoU7fXkcDgwZ84cAEBlZSWOHTsW\n/dmll16Kv/zlL+O+RjA4dsTR8dkDAL6PHadO4cyZMzmrt9CGhoZUXf94tHx9Wr42gNenNgcOHEB7\ne3v068bGxqyeX/Cg0Ol00T9XVlbC6XQCAJxOJ6qqqrJ+PUHQQUrz80qTEeXTqrN+3WJw5swZXHTR\nRUqXkTdavj4tXxvA61ObdevWYd26ddGvm5ubs3p+we96il13WLp0abTgPXv2xK1f5IpXV7Q3dhER\nqULBfosGAgHU19eju7sbDQ0N6OzsxIIFC2AymVBXVwdRFLFo0aKsXzc4zmL3wMDAREsmItKU1tZW\nPP7441k/r2BTT6IoYvfu3WO+v3Xr1km9rqBLPfVkrajAV199hblz507qPYiItMBut8Nutxf/1FMh\n/I+b3gIAXHv11ejpYatxIqLJ0GRQbLnrfVxzzdWYf918HO3uVrocIiJVU31QJFujqLZZ8cH+A7hp\n8U049H6HAlURERWfia5RqD4ohJjbbWWBDzvh945gxowZ6OvvV6AqIqLiY7fbSzMokgm5hxAUBEy/\noBpnzzEoiIgmQ5NBAQAIBmGdMpUHGBERTZImg+Lp7vsAABaLBcMjIwpXQ0SkbqoPimSL2b89fD8A\nwGw2wzM8PKEOtEREWsPF7EQGI/R6PYwGA7xeb2GLIiIqQlzMThAKhvdrm81meHjaHRHRhGk2KBAZ\nRVjMZrjdPMSIiGiitBsUgQCCPV/CYi7jiIKIaBJUHxTJFrM3zn8x/IeRYZjLOPVERARwMTvOpute\nAgB4fvVjmEU9p56IiMDF7OQ8bljKOPVERDQZ2g4KAOXmMo4oiIgmQfNBYTaZGBRERJOg+qAY7yhU\n+a6nkUFHgSoiIipOXMyOIfd6AgBLWXgfxXCQbTyIqLRNdDE74zOzBwcH0dXVhd7eXlitVixcuBAX\nXnhh1m9YCL89fD8ewwkAgKWMU09ERJMxblAMDAzgnXfegcFgwNy5c3HJJZfA5XKhtbUVHo8HK1as\nwOzZswtRa3akAADAUsbFbCKiycgoKL7zne+M+f71118PAPj4449zX1UOuP/jBwDCaxRfRW6P9bpc\nMFVUKFkWEZHqjLtGMXfuXLz88sspf37llVfmtKCccQ8BGF2jAACPj11kiYiyldEaxdNPP42XXnoJ\nl156KZYtW4Zly5ZhYGAANTU1+a5v0syXXY7z//o0/AXXs4mIspZRUGzfvh2zZs3CqVOncODAAfzh\nD39Ac3MzDh8+nO/6JuTRuc9F//y1mZfgixMn0NfXh+nV1QpWRUSkThkFxaxZswAAM2fOxD333IN7\n7rkHHR0d+awrY8n2UWyaNxoUV111Fbq6ulBTU4MTX3xRyNKIiIpKa2srWltbUVtbm9XzJryPYsmS\nJRN9ak6lPOEuYtq0aXjvvfdQXl5eoIqIiIpTwZsCPvHEE+js7MTRo0cn+hIFc8UVV+DcuXMIBoNK\nl0JEpDoTDora2lrU1NTga1/7Wi7ryQuDwYCysjK4XC6lSyEiUp2M1igkScKOHTvw5Zdfora2Fjfd\ndFN06mnq1Kl5LTBXKisrMTg4qHQZRESqk9GIYsuWLTCZTJgzZw62b9+O9evXw+fz5bu2CXv6k/Vj\nvmez2eB0OhWohohI3TIaUdxwww246667AAB33303BgYG8Mc//hGbNm3Ka3ETtfXT9XF3PgEMCiKi\nicpoRJE4epg6dSqmT5+el4LyhUFBRDQxGY0onnvuObz++uuoqqrCkiVLsGzZMlSorGcSg4KIaGIy\nCopXXnkFM2bMgNPpREdHB7Zv346//e1v0ekoJY13cJHMZrNhkEFBRCVsohvuMgqKGTNmAAj/sl29\nejVWr16NtWvXZl9lHgg6HaQMHldZWckRBRGVNLvdDrvdjubm5qyel/GIwu/3AwB0MTuhjx07BqPR\niG9961tZvWm+xfZ6ktlsNjjOn1egGiIidcsoKO65555815FTiXc8AeGgOHnihALVEBGpW0ZBsWPH\nDuzduxeCIODmm2/GmjVrYDab811bTtlstjEb7gKBAHyOAViq1XUHFxFRIWV0e6zFYsEzzzyDq6++\nGm63G4888gj27duX79pyau7cufgo4TS+YDAYnVIjIqLkMgqK6upqbNmyBT09PWhsbMSzzz6L06dP\n57u2nJo9ezZOffml0mUQEalORlNPixcvxo033oimpiY8/fTTcLlcuO222/JdW05ZrVaMeL3w+Xww\nGo1Kl0NEpBrjBkVHRweWLFkCURSxZs0arFmzJu7nnZ2dRXck6tOfjG3hodPpMGXKFAwMDERv9yUi\novGNGxTXXHMNnn/+eYyMjOCSSy6B1WqF2+3GyZMnYbFYcMcddxSizqwk6/UEADarFU6nk0FBRJSF\ncYPCarXiwQcfRDAYxL/+9S84HA5Mnz4dDQ0NEMWMZq4UF/J5oTOaYCkvh9vtVrocIiJVyfjgIkEQ\n0Nrait///vfYtm2ban7hhrwjGP5fmwGE795SS91ERMUiqxPurr76arz88sv45S9/iRdeeCFfNeWe\nzwsgHBQ85Y6IKDtZzR01NTVh165duPzyy+H1ehEIBHD69Glcdtll+apvXJk2BQSA8vJy9nsiopKV\n16aAsrVr1+Lyyy+HJEl4//338fOf/xwdHR3YvXt3Vm+aS8maAibr9QSEg4IjCiIqVXltCihbtmxZ\n9M8zZ87EN7/5TfT19WX1hoWQ7I4ngGsUREQTkdUaRTJqOumOaxRERNmbdFCoicVigcfjSflz/8hw\nAashIlIHzQdF6FwfQsEggPAaRbqpJ9fwSKHKIiJSDc0Hhec/HkTIMYCQzztuUBAR0ViaDIqnP1kf\n/w2PC96XtnExm4hoAjQZFFs/XZ/0+wwKIqLsqaNZU47IQREKhRDKYqMeEVEp0+SIIpWKigp4PB6E\nQiE4+8/h2W3b8I17v6t0WURERa3kRhSx+yje2LkTe1V2pCsRUaGV1IgicY1Cp9MpWA0RkTpoMihS\n9XqyWCzwuN2QpMTuUERElIomgyJVr6eKigq43G74PbzziYgoU5oMilTkFh6844mIKHMlFRQGgwHm\nsjI4BgeVLoWISDVKKigAYN68efj4448BgCMLIqIMlFxQVF9wAQYGBpQug4hINYo2KJqbm/HrX/8a\njz76aNbPHdPrKUaVzQaHI37qiSMLIqLUChYUPT09WLhwISwWC4KRtt8AsGnTJtTV1WHjxo1xj1+1\nahV+8pOfxD02U6l6PQFApc0Gh8MBAPD7/XH/T0REYxUsKKZNm4aWlhYsWbIk+r2uri643W60tbXB\n5/Ph0KFD2LVrF3bu3AkA2Lp1K+67776c1mG12eAcGgIA+Hw+AAwKIqJ0CtbCw2g0wmg0xk3zdHR0\noL6+HkB4BNHe3o6HH34YAPDcc8/h4MGDEAQBixcvTvm6wWCG00aRx5VbLDh37hykYDAaFKdOnYLV\naoXZbMaZM2cmcnl5MTQ0VFT15JqWr0/L1wbw+tTmwIEDaG9vj37d2NiY1fMV7fXkcDgwZ84cAEBl\nZSWOHTsW/dn69euxfn3qKSSZIOiQ0T5rIdyuw2q14kzPGegFAV6vF0D43O8pU6bg/PnzuOiii7K+\njnw5c+ZMUdWTa1q+Pi1fG8DrU5t169Zh3bp10a+bm5uzen7BF7Nj+ytVVlbC6XQCAJxOJ6qqqvL+\n/uU2G9yu8M5secqJLT2IiFIreFDEngWxdOnSaLLt2bMnbv1iMlL1egLCbTyGIh1k5RFFIBDIyfsS\nEWlRwYIiEAigvr4e3d3daGhoQGdnJxYsWACTyYS6ujqIoohFixZl/brBJLe2pur1BAAWgyHaQZYj\nCiIqJa2trXj88cezfl7B1ihEUcTu3bvHfH/r1q2Tel1Bl+EaRURFRTmGXC6EImsU5rKyggWFPJJi\ne3MiUoLdbofdbi/+NQqlVZSXh0cUgh4+nw9ms7lgQeHs/YqjFyJSndIJipFhAPIaxeg+CnNZ2Zg1\nilAohJFBR85LCIE7wIlIfVQfFMnWKJIJHGgBEB5RuIZcCIVC8Pv9KDOPnXoKBoMYieyxICLSiomu\nUag+KIQk8/3pej2Vl5fDMTiI+fPnQxRFiHox+XRQlv/497qGMHyezQaJqHjZ7fbSDIpk0vV6MhqN\nAMIbagwGA/R6ISfrBpLPBylyuy0RkZZoMihSkT7uBgA88sgjABAJChGBQACSJGHq1KkYivSBIiKi\nMNUHRaZrFADg3fF/AABbtmwBEL5NVa/XQ5IknD17FgDw5Zdf5r5IIqIiwDWKbJ4jhC/b5/NBiEw9\nnT1zGgDg8XhyWh8RUbGY6BqFok0Bleb1eqFHeGf2wEB4yml4eFjZooiIiozqRxTJpOr1FOrvg3Rw\nP4DwqCIYDEIUw2sUcnPCkZGRgtU5xA14RKQCmgyKVL2eQme/gq9pOwDAZDIBAARJQiAQwOBg+HjU\nQo4opJgGiURExUr1QZHNYnYs+TZZvRBeo5BPvctlULgH+nP2WkREk8XF7CxFRxSRxWw5IHIZFD6d\n6v96iUhDuOEuSwaDAQAgRm6PHR4Or0342LojJ0KhEILBoNJlEFEOlGxQyCMKvSAgEAhgxBsOCvmM\nCpocSZLgPNundBlElAOaDIp0vZ5ko2sU4RHFCEcUecCFeiIt0GRQpOv1JJOnnuQ1CrfLBaPRGB1R\nhAwGTp0QEUEDG+4meteTvH9BXqMYcTlhtVrjpp546yoRaUlraytaW1tRW1ub1fNUP6KY6F1Pv/nN\nb/DMM8+MrlGMeFFRUcE1CiLSLN71lKVly5bh3nvvhV7QY8Q1BI8UhNVqhc/nQygUwrvvvqt0iURE\nRaFkg0Im6AX4fV54/f7o1FNvby/Wrl2Lc+fOKV0eEZHiNBkUqXo9JRPemR3EyMhINCj6+8M7qh2O\n3J+bTUSkNpoMilS9npLRCwKCkZ3Z8tSTI3L/v9z/iYiolGkyKLIhH1w0MjKCiooKeL1eeLzhvRQM\nCiIiDQTFRG+PlQmCACkYP/Ukb7rjIUZEpCVsCjhBcvdYr9cbDQqv1wsgu13afefO8dAjIipqvD12\ngkS9HoFAAD6fD+Xl5QgEAtGgkP8/EyvqV2Pzf/5nvsokIlKMJoMik15PMr1+dOqpvLwcPp9vQkFx\n8uRJHPxnV9a1EhEVO00GRSa9nmSCoEdQkuDz+aI7s+Upp2wbBLLlBxFpkep7PU2WHqGxaxSRgMhm\nRAGE76DKhM/lhCAas66ViEgJpR0UZ05CcA/BHwhAkiSYzebwiGICU08AIAiZLax7XS6IprKsyyUi\nUoImp56yoRcEeDwemEwmGAyG8BrFREcUQmYjCiIiNWFQCAKGh4ejQRE7osh2jSLTqSciIjVRfVAk\n23CXTa8nIWFE4ff7MeL1oqysDF6vd8zZz8FgMOWitU5I/9fZ9f8+xJkzZzKuTTfJPSJERLG44S5G\nul5PwdMnETj4XvRrvSDAMzwMo9EIURSjeypsNht8Pl/47Oe+3ujjBx3nU762Xp/+r3PlbY34/ve/\nn/YxscrKuI5BRLnDDXeZkgII9Y+2D5ennsrKykbXKLzhQ4wynXqSRxg6jD8CyHbdg4hIaaUXFAn0\ngg6DziGYzebomdk+nw9WqzXjX+ry1FQgEBj3sSaTaVL1EhEVWskHhSAI+PyLL3DppZdCFMVor6fy\n8vKMj0WVAyKToDCbzZOql4io0Eo+KPSCgLPnzuHyyy+Pm3qS+z5lIhoUkjTuY0WxlLau6CBl8HeS\nrcQbDApFkiTuvqeSpMmgyKbXkxi5pXXWrFnRu57koPC6hjJ6DTkoYkcgoVAIra2tYx5bSr9oQnoR\nQwP9OX9d/7AHQ709OX/d8Th54iGVKE0GRTa9nvSR3dRVVVVjpp4mM6I4fvw4vvGNb4x5DSX+Jaw9\nvG2YqJA0GRTZECJ7H8rLy6OL2XJQSIHMpk2iQREzohgaCo9GEs+oKKURBRFpQ8kHhT4mKOQ1Crnl\neEDKdkQx+viRkREAY2+HZVAQkdowKCJBYbFYIAhCdKd2eOopsxGFJEmRaavRoJDDIzEoYqee/r5n\nD44ePTrZSyAiyisGRcyIAgAMBgN6+/rCt8cGwlNJIVGE5E+9+S4QCKCsrCxuPSKToPj2/d/Dhg0b\nkr6mf5jndRNRcdBkUGTb6wkYDQqv14uzZ8+iurp6dEShFzEyOJjyNZIFhXwH1HhTT6mmolwud8bX\nQESUT6oPimRNAdP1ekoUO/UU6+KLL4aUxRpFqhFFYhuQxGBIdWeVDlzLIKLcYlPACRITRhQAsP4H\nP4DRaMx4jUIOCn8GU0+JeLssERUKmwJOkHyGhNE4ejSp1WYLd5LNYkRhMpnSTj1FGwcmBBuDgoiK\nXckHhdxiIvYXuLWiItJyPPMRhRw0iQ0C5f+X3ycxGLQaFG+++Sa+8a1vKV0GEeVAyQfFzAumjfne\n7NmzodfrMx5RSJIEvV4f3dkNjA0I+etSCYpXX301aQsTIlIfTQZFNr2errpkZtzBREcOH8btt98O\nURQhBST09/ejvz+mX5Eu/FcWCoWi00mBQACiKMIQOfgIGJ16koMi1YhiIg3u5Pcu5pDJRzNAIlKG\nJoMim15PugsvRkA0AJFf8BddfDF08uggEMAtK+pQX18/+gSnAz6PByOO83BHAkYOinQjCvnup8Rf\noKFQCIMDA1nt2JZP3Rss4iZ1xRxiRJQdTQZFJoK9p8N/kJcmQpFfbJGNdXq9HoFAAKdOn8Hx48ej\nz2tr3Yve3l4gYUQhTz0lrk0krlkknXrSZd/aI4QQUMRnagvjnB9OROrB/5pTMBgMCEgSzOb4c6vv\n+v6D+MUvfxn3vUxGFKW2RsGgINIO/tecgjw6KDdbxvysoqI87uvoGoUh+zUKrQaFfNsxEakfgyIF\nURQhSVLSo0vNCeERnXrSjx1RJE49Ja5RMCiIqNhpMiiy6fWUirxGYbGMDYrEX/aBQACCIEA0GFLu\nnyi1qScGBZF2aDIosun1lEp4Z7aEMpMJQPwvdH9CJ9lgMBheo9DroyOKxKmnVEGRuFNb1t/fj4sv\nvnjS16EUrlEQaQf/a05BXqOQf7G73KPdXH0+f9xjg8EgBEEIL4CnuOsp1RpFqqD4qrcXw8PD4/aK\nKlbyXVw8qIlI/RgUKej1ekiSFN3/4HQ6oz9LNn0kCELciCLTqadUQSG3GZePVFWbxIAkIvUqyaAI\nHP0ngp9+HPe94L+PxX99rAuiKMLtCR8gNBgTFIn/SpYkCYIgQG80RkcAnsj5FbFBYTQaxw0KV+9X\nOH/+PIYiZ23LR6qOZ/h8P0Yc5zN6bCEk3v2VjvNsX77LIaJJKMmgkP7ZDunwobjveX77s7ivvS/8\nF0S9HsNuNyoqKjA4OJhyOiUYDEIYGozbR+GLnI4X+y/rZEExprZgOFjckamuTH7RAoDk90MKZPbY\nQhg9R3z8EYUkGvJdDhFNgiaDIpteT+mIoojhES+qq6vhdA6NmUaShaeedAlrFPH7JwKBAAwGQxZT\nTy4AmQdFsYlet0rrJ6JRmgyKbHo9pSOKerhHhjFt2jQ4h5wx/0pOEhQ6AfrYNQpp7BqF0WiEJElw\nud3jjizUHhTZTD0RUXEr2qD44IMP8NRTT+HRRx9VrAa9XsRIZEQx6HSOuaNJlmxE4feHN+HFbriT\nRxTX3bQU//s3TwDQ7ogik6mnAwcOFKqcgtPytQG8vlJTsKDo6enBwoULYbFY4v41vWnTJtTV1WHj\nxo1xj1+8eDEAoKwsvtdSIYmiHgFJCo8oYqae/P7kI4rYNQopcupdbAsPeY3iXH8//vXJJ2nfW+1B\nkcnUU3t7e6HKKTgtXxvA6ys1BQuKadOmoaWlBUuWLIl+r6urC263G21tbfD5fDh06BB27dqFnTt3\nAgA2b94Mg0G5hU5RLwJAZI1idEQhJRlR6ARdfPdYaXSqCSi9NQr5Ov0Jf1dEpD5iod7IaDTCaDTG\n3THU0dERPeth1apVaG9vx8MPPwwAeP311/HZZ58p2gpCFMPvXV1djZNffJ5y3n10RCHE7KOQ4kYU\nsWsU8nOA1BvSsr3rqdhEg0Kl9RPRKF2owFtnb7nlFjQ3N0MQBPzqV7/CwoULsXr1ajQ3N6O9vR2P\nPfZYxq/V3Nycx0qJiLRr1apVGT+2YCMKWexUS2VlZXTHs9PpRFVVVVavlc2FEhHRxBT8rqfYs6aX\nLl0aHRXs2bMnbv2CiIiKQ8GCIhAIoL6+Ht3d3WhoaEBnZycWLFgAk8mEuro6iKKIRYsWFaocIiLK\nUMHXKIiISF2KdsPdeFLtv9CCEydOYMaMGVi5ciUaGhqULicnku2jefLJJ1FbW4v7779f9V1mk11f\nVVUVVq5ciZUrV8LhcChc4eR88MEHWL58Oerq6rB582YAwBNPPKGJzy/ZtWnpszt69CiWL1+OFStW\nYMOGDQCy/+xUGRSx+y+8Xi8OHTo0/pNUZvXq1WhpacE777yjdCk5kbiP5uzZs2htbcW7776L6667\nDm+88YbCFU5Osn1C8+fPR0tLC1paWrK+UaPYzJo1C3v37kVbWxv6+vrQ1taGffv2aeLzS7y2I0eO\n4LrrrtNw4pkmAAAFUElEQVTMZ3fllVdi//792LdvH7xeLzo7O7P+7FQZFLH7L2699VZN7qJsaWnB\nihUrsHXrVqVLyQmj0YjKysro1wcPHoTdbgcwuodGzeTri53J/eijj7BixQr89Kc/VbCy3Jg+fTqM\nRiOAcLPMY8eOaebzS7w2vV6PY8eOaeazi92L5vF4JvTfniqDwuFwwGazAQjfYqv2oWGiiy66CP/+\n97+xd+9eNDc348iRI0qXlHNa/Qxjb//+9NNPsW/fPjgcDuzatUvBqnKnu7sb586dQ1VVleY+P/na\nrrrqKs19dm+99Rbmz5+PsrIyTJkyJevPTpVBMdn9F8XOYDDAbDZDEATcfvvtmgwKrX+GAKLXtHbt\nWk18hufPn8ePfvQj/OlPf4LNZtPU5xd7bYD2Prs777wThw8fhtVqhcViwWDkYLVMPztVBoXW91/I\nfZ4AYP/+/ZgzZ46C1eSWPDVTU1ODffv2AdDWZyjvE/J4PNFFbS18hpIk4b777sOTTz6JCy64QFOf\nX+K1ae2zk49zBgCbzQaHw4G2tjYAmX92qgwKre+/ePfdd7Fo0SLcfPPNmDlzJmpqapQuadJi99F8\n/etfx/Hjx1FXV4fa2lp8+OGHuOuuu5QucVIS9wkdOXIENTU1sNvtOHXqFNatW6d0iZOyfft2HDx4\nED/+8Y+xcuVKfP7555r5/BKvrbu7W1Of3TvvvAO73R69nvvuuw+1tbVZfXbcR0FERGmpckRBRESF\nw6AgIqK0GBRERJQWg4KIiNJiUBARUVoMCiIiSotBQUREaTEoiIgoLQYFlbR9+/Zhy5YtSX924sQJ\n7N27N+Vz582bh1deeSWn70lUjBgUVPJiO77GOn78OFpaWpL+rLu7G7fccgveeuutnL5nKh988AGa\nmpqwY8eOCb0f0WQwKIgiNmzYgFtvvRV33nknHA4HnnvuObz44ovRs09ivfbaa3jooYcwMjICv9+P\nffv2obGxEWvWrEFtbS08Hg98Ph/Wrl2LxsZG3HvvvXjhhRfSvqfc0TNRIBDAW2+9hcbGxrQjHKJ8\nYVAQAdi1axcuu+wy7NmzBz/84Q/x7LPP4qGHHsL3vvc97N69e8zju7q6cOONN2L16tXRn5tMJrz5\n5ptobGzEnj178MYbb2D58uVoamrClClTxn3Pbdu2Ja2tqakJOp0O27ZtgyDwP1kqPFHpAoiUFgqF\n8NFHH+Hll1/G3//+dwQCASxduhSp+mV+9tlnOHz4MBobG+H1ejFv3jx8+9vfxrXXXgsgfPCUw+FA\nT08Prr/+egDADTfckNF7JvPee+9h06ZNaGpqws0335zDKyfKDIOCSp5Op8OVV16JBx54ABs3bgQQ\nPqPg/fffRyAQGPP41157Dc8//zxuueUWAOHDbYLB4Jh1h9mzZ0fbjnd3d2Px4sVxP0/2npIkob+/\nH9OnT48+zuVyYfr06ejo6MDvfve7nF47USY4jqWSFgqFoNfrceedd+KLL77AqlWrcOutt+Ltt9/G\ntddei/379+Pee++Ne05TUxOWLVsW/fqaa67Be++9N+a1165di/379+O2225Db28vDAZD9Gc6nS7p\nex4/fhyPPfZY3Ovcdttt+Otf/4rvfve7MJvNOf4bIBofz6OgkvbSSy/B7XbjoYceysvrS5IEvV6P\nDRs24IEHHsBNN92U9vGvvfYapk6dCrvdnpd6iCaCQUEl65VXXsEzzzyDV199FdXV1Xl5j4aGBrhc\nLlxxxRX485//nJf3IMo3BgUREaXFNQoiIkqLQUFERGkxKIiIKC0GBRERpcWgICKitBgURESUFoOC\niIjS+v8FOgUHEJ3jiQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "posterior = BeamingAnglePosterior(o1scenario, efficiency_prior=UniformDistribution(), grb_rate=3/u.gigaparsec**3 / u.year)\n", + "chain2 = posterior.sample_theta_posterior(nburnin=1000, nsamp=1000)\n", + "\n", + "plt.figure()\n", + "theta_bw = 1.06*np.std(posterior.theta_samples)*len(posterior.theta_samples)**(-0.2)\n", + "posterior.get_theta_pdf_kde(theta_bw)\n", + "theta_bin_size = 3.5*np.std(posterior.theta_samples) \\\n", + " / len(posterior.theta_samples)**(1./3)\n", + "theta_bins = np.arange(posterior.theta_range.min(), posterior.theta_range.max(),\n", + " theta_bin_size)\n", + "plt.plot(posterior.theta_grid,posterior.theta_pdf_kde, \\\n", + " color='k')\n", + "plt.hist(posterior.theta_samples, bins=theta_bins, normed=True, log=True, histtype='stepfilled',\n", + " )\n", + "plt.axvline(np.percentile(posterior.theta_samples, 95.0), linestyle='--')\n", + "plt.ylim([1e-3, 1])\n", + "plt.xlim([0,30])\n", + "plt.xlabel(r'Jet Angle, $\\theta$')\n", + "plt.ylabel(r'$p(\\theta | D,I)$')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/grbeams-0.1.0-py2.7.egg/grbeams/beamingangle.py:176: RuntimeWarning: divide by zero encountered in log\n", + " + np.log(self.comp_efficiency_prob(efficiency))\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/matplotlib/font_manager.py:1288: UserWarning: findfont: Font family [u\"'New Century Schoolbook'\"] not found. Falling back to Bitstream Vera Sans\n", + " (prop.get_family(), self.defaultFamily[fontext]))\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAESCAYAAADjS5I+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXX+P/DXvVxBSbxsIpArWmNoKgEKKHiVZZRUdHQc\nnRydmean8yuzxGkeNWXZb75TTZlj05Tp49tui9m4IlosAooQSyQiZpILqYgoXFZZ7vL7g7ixXu5+\n7j28no9Hj+Dcc895n0756nM+y5FotVotiIiI+iAVugAiIrJvDAoiItKLQUFERHoxKIiISC8GBRER\n6cWgICIivRgURESkF4OCiIj0kgldQF8KCgpw9OhRaDQaPP/880KXQ0Q0YNm8RVFRUYHg4GC4urpC\no9HoticmJiIqKgobN24EABw8eBCbN2+Gt7c3KisrbV0mERH9xOZB4eXlhfT0dISFhem2FRUVobGx\nEVlZWWhtbUVBQYHuM64wQkQkLJs/enJ2doazs3OXAMjNzUVsbCwAIDo6Grm5uVi0aBH+/ve/o62t\nDSNGjLB1mURE9BO76KNQKpUYP348AEAul6O0tBShoaEIDQ3V+720tDRblEdEJDrR0dEG7ytYUEgk\nEt3PcrkcdXV1AIC6ujq4u7sbfJz7UvdCc60cXlvfxRAfXzzxxBMYPXo0EhMTsWULsGWLhQu3oU2b\nNuG1114TugyrEfP1ifnaAF6fozP2f7IFGx6r1Wp1j5/Cw8N1haempnbpvzDWiBEjdJ3fL7xgfp1E\nRAOdzYNCpVIhNjYWxcXFmDdvHvLz8xEUFAQXFxdERUVBJpMhJCTE4ONpunV2dw4KIiL6WUZGBraY\n8JjF5o+eZDIZUlJSemzfvn27SceTSiTQdPrd29sbt2/fNrE6+xIeHi50CVYl5usT87UBvD5HpVAo\noFAoHOfRk7V4eXmJJigiIiKELsGqxHx9Yr42gNc30IguKDw9PUUTFERE9sAuhseao3sfhZeXF6qr\nqwEAXPmDiOhnGRkZyMjIQGRkpFHfc/gWhbTTMFugvUWhVCqh0WgcemgsEZGlKRQKkzqzHT4oupPJ\nZLjrrrtQW1srdClERKIguqAA2E9BRGRJouujALr2UxARUTv2UXQipiGyRESWwj6KTjqCgp3ZRETm\nE2VQeHp6orq6mms9ERFZgCiDgo+eiIgsx+GDoq/ObAYFEVFXpi4K6PBB0VtnNofHEhH1xM7sTjg8\nlojIckQbFLdv3+ZaT0REFiDqoODwWCIi84kyKDw9PVFTUyN0GUREouDwQdHbqKdhw4ahoaEBKpVK\ngIqIiOwTRz113iaVQi6XQ6lUClAREZF9Gtijntpa0VJTjbrvStCibH/k1DE7m4iIzCOKoNDWKVGf\nuBq1f/kj2mp/DopXXnEVuDIiIscniqDojYeHB955Z6TQZRAROTzRBoWnp6fQJRARiQKDgoiI9BLl\nG+4ABgURUXd8w103Hh4eNq6EiMi+Dezhsb3w9PTElCn7hS6DiMjhiTooRo9+V+gyiIgcnqiDghPu\niIjMJ9qg8PDwYFAQEVmAaIOCLQoiIssQbVB4eHigpqYG2j6GzxIRkWFEGxQuLi4AtqCxsVHoUoiI\nHJpogwIA2tr+xsdPRERmcvig6GtmdgcGBRFRO764qA98JSoRUTvOzO4DWxREROZhUBARkV6iDoqI\niBQGBRGRmRx+mXF9Fi4sRE2NUugyiIgcmqhbFFzGg4jIfKIOCi7jQURkPgYFERHpJfqg4DwKIiLz\niDooPvhgHFsURERmEvWop9dfd8fQoQwKIiJziLpFAQDNzc1oa2sTugwiIofl8EHR36KAHe+lICIa\n6LgoYB84l4KIqB0XBewDh8gSEZlH1EHx/PMMCiIic4k6KLZs4VwKIiJziTooAPZREBGZS/RBwUdP\nRETmYVAQEZFeogsKrVQKtVoN7U/zK9hHQURkHnEt4aHVoHHPu2hyGYwh85di2+5fIDSUfRREROYQ\nXYui7XgyWo/tg7atFS+8wEdPRETmEl1QdMdHT0RE5hkQQcEWBRGR6UQfFO7u7qipqYFGoxG6FCIi\nhyT6oBg0aBBcXV1RX18vdClERA5J1EHx/PPtf2c/BRGR6UQdFB2r6XIZDyIi04k6KDqwQ5uIyHR2\nGxRnzpzBkiVL8P3335t9LAYFEZHpbBYUFRUVCA4Ohqura5cRSImJiYiKisLGjRu77H///fdjyZIl\nFjk3+yiIiExns6Dw8vJCeno6wsLCdNuKiorQ2NiIrKwstLa2orCwEElJSTh48CAAQKvV6tZsMgf7\nKIiITGezoHB2doZcLu/yB39ubi5iY2MBANHR0cjJycGCBQuQkJCA8vJyfPXVV/joo4/Q1tZm0jk7\nOrP56ImIyHSCLgqoVCoxfvx4AIBcLkdpaanus9GjR+Pjjz/u9xgaTe8tjubmZrzwArB27XVIpVJc\nvXoV169ft0zhNlJfX+9wNRtDzNcn5msDeH2O5tSpU8jJydH9Hh8fb9T3bR4UEolE97NcLkddXR0A\noK6uDu7u7kYfTyqVQN3L9sGDBwMA/P39ERAQgFOnTsHf39+kmoVy/fp1h6vZGGK+PjFfG8DrczTL\nli3DsmXLdL+npaUZ9X2bj3rq3O8QHh6uKzg1NbVL/4UlsY+CiMh0NgsKlUqF2NhYFBcXY968ecjP\nz0dQUBBcXFwQFRUFmUyGkJAQo4+rMaCzm30URERARkYGtnR03hrBZo+eZDIZUlJSemzfvn27WceV\nSnp/9NQZg4KICFAoFFAoFPb/6MmWuNYTEZH5RB0UHS2sjkl+TU1NgtZDROSIHD4o+uqj0Gi1aCi/\nhOZbVZBIJPDz80NFRYWNqyMish+m9lE4fFBIOw237azhqf+Dmsd+i+byiwDah8mKaVw0EZGxFAqF\nfXdm25xaDUANoL3F4efnx6AgIjKBw7coDOXv789HT0REJhB1UGz7fq3uZz56IiIyjcMHhb4Jd9vL\nGBRERB3Ymd0P9lEQ0UBname2wweFodhHQURkmgEVFGxREBEZb8AEhVwuh1qt1i1rTkREhnH4oNDX\nmf3EhF26nyUSCcaOHYvLly/boCoiIvvDzuxeJN67q8vv48aNw8WLF61dEhGRXWJntgHGjh2LS5cu\nCV0GEZFDGVBBMW7cOAYFEZGRGBRERKSXwy8KaMirUDswKIhoIMvIyEBGRgYiIyON+p7Dtyj0dWZ3\nXusJ+DkotEaECxGRWFh9mfHa2loUFRWhsrISbm5uCA4OxogRI4w+oS1tL1uL55Gr+10ul+Ouu+7C\n9evXcffddwtYGRGR4+g3KKqrq3Hs2DEMGjQIEyZMwKhRo9DQ0ICMjAw0NTVh9uzZCAgIsEWtFjFp\n0iSUlpYyKIiIDGRQUPz2t7/tsX3q1KkAgO+++87yVVlRYGAgSktLERsbK3QpREQOod8+igkTJuDT\nTz/t8/OJEydatCBr6wgKIiIyjEF9FNu2bcPu3bsxevRoREREICIiAtXV1QgNDbV2fRYXGBiIzz77\nTOgyiIgchkFBsXfvXowdOxZXr17FqVOn8MYbbyAtLQ1nzpyxdn1maV/raUqXbYGBgTh79iy0Wi0k\nBr7LgohoIDMoKMaOHQsAGDlyJJYvX47ly5cjNzdX/5dsRN88iva1nv7TZZuPjw+kUikqKyvh6+tr\n5eqIiOyHzedRhIWFmfpVizL0DXcdJBIJpk6diqKiIitVRERkn2y+KOCrr76K/Px8nD171tRDCCY4\nOBiFhYVCl0FE5BBMDorIyEiEhobCz8/PkvXYRHBwMAoKCoQug4jIIRgUFGq1Gnv27MHWrVvx9ddf\nA/j50ZOnp6f1qrOSkJAQtiiIiAxkUFBs3rwZLi4uGD9+PPbu3Yu1a9eitbXV2rWZbdv3awEnGdRq\nNTQajW57QEAAGhoacPPmTQGrIyJyDAaNepo2bRoWL14MAFiyZAmqq6vxn//8B4mJiVYtzlzby9bi\nyd0Po9l7BFziFsMjaDqA9g7tBx54AIWFhZg/f77AVRIR2TeDWhTdWw+enp7w8fGxSkGWpi49jbas\nr4CG+i7bQ0JCkJeXJ1BVRESOw6AWxa5du7B//364u7sjLCwMERERGDp0qLVrsyypVPf4SSqVIiIi\nAm+99ZbARRER2T+DguLzzz+Hr68v6urqkJubi7179+LIkSO6x1FCMvTFRXeS9qC15BsMmhEFz2mh\niIiIwOrVq6FWq+Hk5GTlKomIhGfqhDuDgqJjBvOwYcMQFxeHuLg4JCQkGF+lFUglEqgN2E9dXAB1\ncQFkY+8BAAwfPhx+fn4oKSnRrYRLRCRmCoUCCoUCaWlpRn3P4BZFW1sbAHRZH6m0tBTOzs749a9/\nbdRJbaV9rae+zZw5EydPnmRQEBHpYVBQLF++3Np1WEX7Wk99mzVrFlJSUvDoo4/aqCIiIsdjUFB8\n8cUXOH78OKRSKWbNmoVFixZhyJAh1q7N6mbOnGnSuidERAOJQcNjXV1d8eabbyIwMBCNjY14/PHH\nkZmZae3arO6ee+7BnTt3UF5eLnQpRER2y6Cg8Pb2xubNm1FRUYH4+Hi8/fbbuHbtmrVrszqJRII5\nc+bg+PHjQpdCRGS3DAqK6dOn4/nnn0dISAi2bduG9evXw83Nzdq12URMTAxSUlKELoOIyG71GxQd\nLyiSyWRYtGgRXnnlFbz11ltYuHAhACA/P9+6FZph2/dre2xrPZWKW4c/R3XJt9BqtYiOjkZqaiq0\nBs7HICIaaPrtzJ40aRLeeecdNDc3Y9SoUXBzc0NjYyPKy8vh6uqKBQsW2KJOk2wvW9tj5JOq4BRU\nBacw+PfrgUlTERAQAFdXV5w9exaTJ08WqFIiIvvVb1C4ubnh4Ycfhkajwfnz56FUKuHj44N58+ZB\nJjNo0JTdi4mJQWpqKoOCiKgXBr+4SCqVIiMjA//+97+xY8cONDY2WrMum2I/BRFR34x6w11gYCA+\n/fRT/OMf/8CHH35orZpsLiYmBidOnEBTU5PQpRAR2R2jgiI5ORlPPvkkPv74Y7S0tEClUuHKlSvW\nqs0ghi4KqI+npyeCg4ORmppqgYqIiOxTRkaGSZOMjepkSEhIwLhx46BWq/H111/jhRdeQG5urqCP\nbfQtCtjfWk+dJSQk4MCBA1i0aJFlCiMisjNWXRSwQ0REhO7nkSNHYunSpXb9OtH+1nrqLCEhAS++\n+CKXHSci6saoR0+9cZQ33fVn3Lhx8PPzQ05OjtClEBHZFbODQkyWLVuGzz77TOgyiIjsCoOik4ce\negh79uzp8Y5wIqKBjEHRSUBAACZOnIhjx44JXQoRkd0QdVD0ttZTf1avXi2qOSJEROYSdVBsL9Mf\nFBqNBhqNpsuCgMuXL0daWhoqKiqsXR4RkUMQdVDo05KShFv/ux3VGV9Co9HotsvlcqxYsQI7d+4U\nsDoiIvsxYINCe+0yWg99ClVxz2XS169fj507d7JTm4gIAzgo9Jk0aRICAwPxxRdfCF0KEZHgGBR9\n2LBhA9544w2hyyAiEpyog8KYtZ66W7BgAW7cuIG8vDwLVkRE5HhEHRTGrPXUnZOTE9avX89WBREN\neOJ4RZ0ZtLVKNJSdB7Qa3DVhYpe39v3xj39EQEAAbty4AV9fXwGrJCISjqhbFIZQ5WWhLnE1Gna9\n1mWYLAB4eHhwqCwRDXgDPij689hjj+Htt9+GSqUSuhQiIkEwKPoRGBiIMWPGGP2iDyIisbDboEhL\nS8PLL7+MJ554wuRjGLvWk+pOE1rqaqHu1npYuXIlPvnkE5PrICJyZDYLioqKCgQHB8PV1bVLX0Bi\nYiKioqKwcePGLvtHR0fjqaee6tFvYIz+1nrqTFNWipqn1uH2P/+G1prqLp8tX74cBw8exJ07d0yu\nhYjIUdksKLy8vJCeno6wsDDdtqKiIjQ2NiIrKwutra0oLCxEUlISDh48CADYvn07Vq1aZZsC1Wpo\nyn+A5tqVHh/5+fkhKCgIX331lW1qISKyIzYbHuvs7AxnZ+cuK7Xm5uYiNjYWQHsLIicnB+vXrwcA\n7Nq1CwUFBZBKpZg+fXqfx9VotH1+ZqraujrUtHV9/DR79mzs2bMHoaGhFj9fX+rr63H9+nWbnc/W\nxHx9Yr42gNfnaE6dOtXlNc/x8fFGfV/QeRRKpRLjx48H0L5qa2lpqe6ztWvXYu3a/h8dSaUSqC1c\nl3zYMAwZ3vVd4KtWrUJkZCR8fX0hldqmIXb9+nX4+/vb5FxCEPP1ifnaAF6fo1m2bBmWLVum+93Y\nwTk278yWSCS6n+VyOerq6gAAdXV1cHd3t3U5BpswYQLc3d1RWFgodClERDZl86DQarW6x0/h4eG6\nZEtNTe3Sf2EJ5qz11JuFCxfi8OHDFj0mEZG9s1lQqFQqxMbGori4GPPmzUN+fj6CgoLg4uKCqKgo\nyGQyhISEGH1cjbbvPgpz1nrqzcKFC3Ho0CGLHpOIyFYyMjKwZcsWo79nsz4KmUyGlJSUHtu3b99u\n1nGlEsv3UfQlPDwc165dQ3l5OUaPHm2jsxIRWYZCoYBCobD/PgpHJpPJ8OCDD7JVQUQDCoOiD1qt\nFspzJajJTkfDtR+hVquh0WiQkJCgm+dBRDQQOPwy4/r6KEyhbb6D5h8voaX8B7ScSEFryiHIQmeh\nwU2OIfN/hdjYWKxZswa1tbWQy+UWPTcRkTVlZGQgIyMDkZGRRn3P4VsU0k7Dbbszdq0nAEB9LRo2\nP4r65x5Da0r7IyZV/km0pR+BtqUZQ4cORVRUFI4ePWpqyUREglAoFCZ1Zjt8UOhjzFpPxkhISGA/\nBRENGKIOCmtZuHAhjh49iqamJqFLISKyOocPCkv3URjC19cXM2bMYKtCAHcqK1B/+Qc0VVwTuhQi\nh2PqPAqHDwp9fRTWtHr1ary7823cOrwXNWdPd1nskKyn+XIZlI/+Bs3nzwpdCpHDYR+FjS1evBh5\n3xThyuv/A9X5EqHLISKyGlEHhaXXeurM1dUVC2Pm4sC16v53JiJyYKIOCkuv9dTdbxYtwr5rt616\nDiIioTl8UNiyM1srleHO7Sq01rcvjT5regiUbSqcvVxusxqIiEzFzmwbaPznU7j95J/QWHq6/dxS\nKRb7e2Fv5kmb1UBEZCp2ZtuAtrYG2sprgKpNt+1Xd3ti38kcqNW2WsOWiMi2GBRmutdtCHzc5cjM\nzBS6FCIiqxB1UJi01pMJ4meEcPIdEYmWqIPCWms9dffLkAdw+PBhTrqzkUuNzXjx3Q9QXc2hyUS2\n4PBBIcQSHt0FjhkFjUaDs2c5W9jatFot1n1zEam5edi0aZPQ5RA5FI56EpBEIuH7tG3k5DdF0Gi1\n+OLVl7B//34olUqhSyJyGBz1JLD58+fjyy+/FLoM0UvOOonF/p7wlA/DjBkzkJGRIXRJRKLHoLCQ\nqKgofPPNN6ivrxe6FFHLyCtApPcwAEBMTAxSU1MFrohI/EQdFNZc66m7u+66C6GhoRwma0W3bt1C\nRdUtTJa7AgDCwsJQUFAgcFVE4ifqoLD2Wk/dxcbGIiUlRe8+dRcvoLa4EI3XuOyHsUpKSnBfwDg4\n/dQvNWXKFJw5c4aTHYmsTNRBYWtxcXH9BkVb2TnUPb0ObVWVNqpKPEpKSnDf+HG63+VyOXx9fXHh\nwgUBqyISP4cPCnsYHtshKCgIN2/exNWrV4UuRZQ6WhQAoG1uQl3ZeQTeMwFnzpwRuDIix8DhsXZA\nKpUiOjq631YFmaY9KAIAAE1v/A9qH38IY4cOQVlZmcCVETkGDo+1E4b0U5DxtFotSkpKMDFgXJft\n4/x9GRREVibqoLDVWk+dxcbGIjU1FRqNxubnFrNr165h8ODB8PZw77I9wN+PQUFkZaIOClut9dTZ\nmDFj4OnpiaKiIpufW8xKSkowefLkHtsZFETWJ+qgEMqDDz6II0eOCF2GqPQVFH5enqiurkZjY6MA\nVRENDAwKK1iwYAGSkpKELkNU+goKqVSKUaNG4ccffxSgKqKBgUFhBbNmzcKFCxdw48YNoUsRjbNn\nz/YaFAAYFERWxqCwgkGDBiEuLg7JyclClyIKGo0G586dQ2BgYK+fjxw5kkFBZEUOHxT6JtzZcq2n\n7hYsWIADBw4Idn4xuXTpEry9vTFs2LBePx81ahQnORIZgBPuemHrtZ46W7x4MU6cOIHKyp+X6mhp\nacFjW1/H9PRi7DlyVLDaHE1JSQkmTZrU5+dsURAZhhPu7IybmxuWLFmC9957D0D7hLE1a9agtrER\nO4IC8Let27hGkYH66sjuwBYFkXUxKKwoMTER27ZtQ2VlJV588UVcvHgRu57+C4I9hmLdit/g5Zdf\nFrpEh9BfULBFQWRdDAoL0mq1UKlUUKlU0Gq1mDx5Mh5//HGMHz8eu3fvxr59+zDY2RkAsGbpYuzb\nt4/j/w3AFgWRsBgUFtR8qwpV27bg1v9uh0qlAgA888wzOH/+PIqLizFy5EjdviO8vREWFsb5Fv1o\na2tDWVkZJk6c2Oc+Hh4eaGtr49sFiaxE1EEhxFpP6tJvoSk712Xb3XffjUGDBvXYd8GCBRYZQqtW\nq6FSqUS5vtR3332HMWPGYMiQIX3uI5FI4O/vj4qKChtW1j+VSsWXKpEoiDoohFjryRjz58/HsWPH\n9P4Br1ap0FB+CQ3ll6Bqael1n7rThah65VnUnvnGWqUK5ttvv8XUqVP73W/EiBF2NcFRVl2Fqlee\nhTL7OLR29M4UIlOIOijsXUBAADw9PfHNN33/Aa9WqfDDS8/gZOJaqJoaet+psQ6q7FSgSXz9HadP\nnzYoKHx9fe0qKKBRQ5WdCk0FO9nJ8TEoBNbf61MvXrqEmA/2Y03yCfzrjf/YsDL7cPr0aUybNq3f\n/ewuKIhEhEEhsLi4OHz11Vd9fv7SSy9hzdRf4HCCAlu3v47y8nIbVicsrVZrVIui8+RGIrIcBoXA\nZs+ejfz8fDQ09Hys1NjYiIOHDuGhyffCf6grVj/0W/zrX/8SoEphXLt2DVqtFv7+/v3uyxYFkfWI\nOiiEXOvJUEOHDkVoaCgyMzN7fJaUlIRZs2bB5672ET9//tPD2L17N1r66NQWm+zsbMycORMSA96L\nbm+d2URi4vBBoW9RQCHXejJGX/0Ue/bswZrf/U73+5jRozF58uQBM/eiIygMwRYFUf+4KKADi42N\n7dFPUVZWhrKyMsyfP7/L9t///vd4//33bVidcE6cOMGgILIgLgooJLUazdW3oWprNenrQUFBuHnz\nZpf1it5//30sWbIEzj8t+dFh6dKlyM7OFv2SFdevX8eVK1cQGhpq0P4+Pj6oqqoS5aRDIqExKCyg\n+dNduP3Xh1H77KPQVhn/f7VOTk6IiYnRPX5Sq9X44IMPsGLFih77Dh06FCtWrMA777xjdt32LCkp\nCfPnz+91RntvXFxc4ObmhurqaitXRjTwMChMoNVo0Xj1CtQd/SMtzdDeuAZt5TWTjxkXF4djx44B\nAI4ePQo/Pz/cd999ve775z//GTt37hTNgoLqtjY0Xr3S/s+0rQ0A8Mknn2DJkiVGHYcd2kTWIeqg\nsNZaT42vPoPqDQ+h4SnLHT8hIQGpqam4evUqtm7divXr1/e575QpUzB79mw888wzolgeQtV8B8rX\nnoPyteegar6DwsJCXLx4EQkJCUYdh/0URNYh6qCw2lpPahXQ0gyY2CfRGy8vL2zYsAHBwcGor6/H\nypUr9e7/xhtvICsrC7NmzcKeL1PQonbwZ/OtLe1/Adi2bRs2bNhg8GOnDgwKIusQdVAIR4vWWiXu\nVFVCbUSYPPfcc9izZw/S09P7/UPS29sbeXl5SExMxEdJR/F/iy6KonXx49WrOHr0KFYt/RXuVFWi\n7U6Twd9lUBBZB4PCCtTflaDmL39AzSvPoLWu1uDvSaVSKBQKyOVyg/aXyWRYunQp9v/rFVQ0t+LQ\n8Z6T9hzNjp07sXr1akjzs3D7yYfRXGF4v8+IESO4jAeRFTAorEGrgbbqBjS3qyx6WA3aR0R1bzkM\nksmwYYIfdu3dZ9Hz2dqdNhU++Gg3NmzYADQ1GD2CjOs9EVkHg8JBaKsqUPf+m7i1cyuaKnrOoYjx\nccf5y5cdetHAAxfKETZjBsaOHWvS9znqicg6RB0UQq/1pK2rQWNOBm4fO4DGG9fNO1hbG9pSD6E1\nLQnoZVLZIKkE0TOm4+jRo+adR0CfnbuIP/3hDyZ/ny0KIusQdVAIvtbTnSbc2fFPNL39CrStlhsh\n1ZeY8Bm6uRiOprLxDsqU9YieO8fkY7BFQWQdog6KgSZ82lScOnXKIUc/pVy+DsVo3x5Llhhj+PDh\nqK6u5nuqiSyMQSEwjUYDtVptkT/c7vYZDplMhsuXL5tfmI0dL69AzJj+3zuhj0wmg4eHB6qqLDuI\ngGigY1AIrO7cGVRt/3+o/Tbf7GNJJBKEhYUhJyfHApXZjlarRX7FLczwH272sThElsjyGBRCa7mD\ntvQjQM0tixwuPDwcubm5FjmWrVwoK8MQmRP8h7r2+Kz16mXUnMpAw4+XDToWO7SJLE/UQWGttZ7s\nWVBQEE6fPi10GUY5lZuLED/vXj9r+ufTaPjHX6CuVRp0LHZoE1meqIPCams92bEpU6aguLjYoTq0\nv87LQ4hv70FhLHtpUWi1WnxcXoXScnG/N4QGBrsNiry8PLz22mt44oknhC5FMKdOnTL6O8OHD8fg\nwYMd4sVGHddXXHwG93t7WOSY9tKi2P3Ff/HK+WtYvfXfUKlUQpdjcab8u+lIxH59xrJZUFRUVCA4\nOBiurq5d3kKWmJiIqKgobNy4scv+06dPBwAMHjzYViXaFbVajcqKCqh+ej+DMTpaFfYuJycHGo0G\n586fxy+8DFvfqj/20qL4/PAR/H3SaMhdXZGdnS10ORbnaAMmjCX26zOWzYLCy8sL6enpCAsL020r\nKipCY2MjsrKy0NraisLCQiQlJeHgwYMAgE2bNhm91LRYtN6uwpLbP6Bm83poLp436ruOEhQAcPHi\nRXh5ecHN2TL32R5aFK2trThz/nsohsvxy+BpSE5OFrQeInPJbHUiZ2dnODs7d3l2npubi9jYWABA\ndHQ0cnJ0HVM1AAAHP0lEQVRydC/s2b9/P3744Qc4OTnZqkSb0jTUo7a0GOqGegCAuqYataXFkAx2\nhdu48e37XPwe2ts3jT72lClTcOTIEYvWa4z6S2XQ3GnCIE9vuPrqnxtRUlKCyYGBFju3PbQoCgsL\nMdLXF8MGOSHk3gnYcSpP0HqEotFo0FD2HbQqFQbfPQoucss8XiTbk2ht3Os5Z84cpKWlQSqV4qWX\nXkJwcDDi4uKQlpaGnJwcPPvsswYfKy0tzYqVEhGJV3R0tMH72qxF0UEikeh+lsvlqKurAwDU1dXB\n3d3dqGMZc6FERGQam4960mq1usdP4eHhulZBampql/4LIiKyDzYLCpVKhdjYWBQXF2PevHnIz89H\nUFAQXFxcEBUVBZlMhpCQEFuVQ0REBrJ5HwURETkWu51w15++5l+IwZUrV+Dr64u5c+di3rx5Qpdj\nEb3No9m6dSsiIyPxu9/9zuGXBu/t+tzd3TF37lzMnTsXSqVhS5DYq7y8PMycORNRUVHYtGkTAODV\nV18Vxf3r7drEdO/Onj2LmTNnYvbs2XjkkUcAGH/vHDIoOs+/aGlpQWFhodAlWVxcXBzS09Md9kVE\n3XWfR1NVVYWMjAycOHECU6ZMwYEDBwSu0Dy9zRO6//77kZ6ejvT0dKMHatibsWPH4vjx48jKysLN\nmzeRlZWFzMxMUdy/7tdWUlKCKVOmiObeTZw4EdnZ2cjMzERLSwvy8/ONvncOGRSd51/ExMSIchZl\neno6Zs+eje3btwtdikU4OztDLv959nVBQQEUCgWAn+fQOLKO6+v8JPfcuXOYPXs2nn76aQErswwf\nHx/dS6VkMhlKS0tFc/+6X5uTkxNKS0tFc+86z0Vramoy6b89hwwKpVKJYcOGAWgfYuvoTcPu/P39\nceHCBRw/fhxpaWkoKSkRuiSLE+s97Dz8u6ysDJmZmVAqlUhKShKwKsspLi7GrVu34O7uLrr713Ft\n9913n+ju3eHDh3H//fdj8ODB8PDwMPreOWRQmDv/wt4NGjQIQ4YMgVQqxYMPPijKoBD7PQSgu6aE\nhARR3MOamhps2LAB7777LoYNGyaq+9f52gDx3buFCxfizJkzcHNzg6urK2prawEYfu8cMijEPv+i\noaFB93N2djbGjx8vYDWW1fFoJjQ0FJmZmQDEdQ875gk1NTXpOrXFcA/VajVWrVqFrVu3Yvjw4aK6\nf92vTWz3rrW1VffzsGHDoFQqkZWVBcDwe+eQQSH2+RcnTpxASEgIZs2ahZEjRyI0NFTokszWeR7N\nL3/5S1y+fBlRUVGIjIzE6dOnsXjxYqFLNEv3eUIlJSUIDQ2FQqHA1atXsWzZMqFLNMvevXtRUFCA\nv/71r5g7dy4uXrwomvvX/dqKi4tFde+OHTsGhUKhu55Vq1YhMjLSqHvHeRRERKSXQ7YoiIjIdhgU\nRESkF4OCiIj0YlAQEZFeDAoiItKLQUFERHoxKIiISC8GBRER6cWgoAEtMzMTmzdv7vWzK1eu4Pjx\n431+995778Xnn39u0XMS2SMGBQ14nVd87ezy5ctIT0/v9bPi4mLMmTMHhw8ftug5+5KXl4fk5GR8\n8cUXJp2PyBwMCqKfPPLII4iJicHChQuhVCqxa9cufPTRR7p3n3S2b98+rFu3Ds3NzWhra0NmZibi\n4+OxaNEiREZGoqmpCa2trUhISEB8fDxWrlyJDz/8UO85O1b07E6lUuHw4cOIj4/X28IhshYGBRGA\npKQkjBkzBqmpqXj00Ufx9ttvY926dVi9ejVSUlJ67F9UVIQHHngAcXFxus9dXFxw6NAhxMfHIzU1\nFQcOHMDMmTORnJwMDw+Pfs+5Y8eOXmtLTk6GRCLBjh07IJXyP1myPZnQBRAJTavV4ty5c/j000/x\n5ZdfQqVSITw8HH2tl/nDDz/gzJkziI+PR0tLC+69916sWLECkydPBtD+4imlUomKigpMnToVADBt\n2jSDztmbkydPIjExEcnJyZg1a5YFr5zIMAwKGvAkEgkmTpyINWvWYOPGjQDa31Hw9ddfQ6VS9dh/\n3759eOeddzBnzhwA7S+30Wg0PfodAgICdMuOFxcXY/r06V0+7+2carUat2/fho+Pj26/hoYG+Pj4\nIDc3F6+//rpFr53IEGzH0oCm1Wrh5OSEhQsX4tKlS4iOjkZMTAyOHj2KyZMnIzs7GytXruzyneTk\nZEREROh+nzRpEk6ePNnj2AkJCcjOzsb8+fNRWVmJQYMG6T6TSCS9nvPy5ct49tlnuxxn/vz5+Oyz\nz/DQQw9hyJAhFv4nQNQ/vo+CBrTdu3ejsbER69ats8rx1Wo1nJyc8Mgjj2DNmjWYMWOG3v337dsH\nT09PKBQKq9RDZAoGBQ1Yn3/+Od58803897//hbe3t1XOMW/ePDQ0NOCee+7Be++9Z5VzEFkbg4KI\niPRiHwUREenFoCAiIr0YFEREpBeDgoiI9GJQEBGRXgwKIiLSi0FBRER6/X99405mUqmnggAAAABJ\nRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "posterior = BeamingAnglePosterior(o1scenario, efficiency_prior=JeffreyDistribution(), grb_rate=3/u.gigaparsec**3 / u.year)\n", + "chain2 = posterior.sample_theta_posterior(nburnin=1000, nsamp=1000)\n", + "\n", + "plt.figure()\n", + "theta_bw = 1.06*np.std(posterior.theta_samples)*len(posterior.theta_samples)**(-0.2)\n", + "posterior.get_theta_pdf_kde(theta_bw)\n", + "theta_bin_size = 3.5*np.std(posterior.theta_samples) \\\n", + " / len(posterior.theta_samples)**(1./3)\n", + "theta_bins = np.arange(posterior.theta_range.min(), posterior.theta_range.max(),\n", + " theta_bin_size)\n", + "plt.plot(posterior.theta_grid,posterior.theta_pdf_kde, \\\n", + " color='k')\n", + "plt.hist(posterior.theta_samples, bins=theta_bins, normed=True, log=True, histtype='stepfilled',\n", + " )\n", + "plt.axvline(np.percentile(posterior.theta_samples, 95.0), linestyle='--')\n", + "plt.ylim([1e-3, 1])\n", + "plt.xlim([0,30])\n", + "plt.xlabel(r'Jet Angle, $\\theta$')\n", + "plt.ylabel(r'$p(\\theta | D,I)$')" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "0.75237100333564122" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.percentile(posterior.theta_samples, 95.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/daniel/.virtualenvs/jupyter/lib/python2.7/site-packages/matplotlib/font_manager.py:1288: UserWarning: findfont: Font family [u\"'New Century Schoolbook'\"] not found. Falling back to Bitstream Vera Sans\n", + " (prop.get_family(), self.defaultFamily[fontext]))\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEYCAYAAAC3LjroAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHXiPvDnzIXLMMwMcjMxdTVTM0kTClAGHMRMBc2w\nTUtKsr5d6aeutWqX3VLzVttubX5Tu62h37KLa1rGcokh75opXta00hJQLjLMDMhlZs7vD3SSAAEd\nPMzwvF8vXw5zzvnMcw4DD3PmzDmCKIoiiIiILiGTOgAREXU+LAciImqC5UBERE2wHIiIqAmWAxER\nNcFyICKiJjyqHN577z307dsXqampUkchInJrbl8OM2bMcN6eOHEisrKyJExDROQZ3L4cBEFw3u7W\nrRvkcrmEaYiIPIPblwM/4E1E5HoKqQNciZMnTzp3Jx07dgwGgwEDBgzAypUrJU5GROQZ3LIc+vTp\ng9zcXABAWloa3n33Xec0URT5aoKI6Cq1abfS7NmzodfrMWvWrEb3FxcXIyEhASNHjkR2djYAwGq1\nIjk5GbGxsVi7di0AwG63IzU1FXq9HsuWLQMA7N69GyNGjIBer8ecOXOcY65YsQKxsbGYPn067HZ7\nu1Zmy5YtmD59OnJycjBlypR2LUtERL9ptRz279+PqqoqGI1G1NbWYt++fc5pS5YswaJFi5CZmYmF\nCxcCAFavXo2pU6fCaDRizZo1sNls2LRpEwYNGgSj0Yj8/HyUlJQ4//o3Go04e/YsDh8+jNLSUnzz\nzTfIz89HeHg4Nm7c2OoKXPqqYfz48cjPz0dhYSE2bNhwJduDiIjQhnLYuXMnEhMTAQCjR4/Gjh07\nnNMKCgoQFRUFlUoFjUYDi8XinF8QBAwdOhRHjx5tNMaoUaOwe/duhISEwMvLCwCgVCohl8uxd+9e\nxMfHAwASEhIaPRYREV07rb7nYDKZ0K9fPwCAVqvFkSNHnNMcDofztlarhclkgslkgkajAQBoNJom\n912c76KDBw+irKwMAwcOxP79+1uc76KLu6+IiKh9EhIS2jxvq+Wg1WphNpsBAGazGTqdzjlNJvvt\nhUdlZSUCAgKg0+lgNpsRFBQEs9nc6L6LY/Tv3x8AcO7cOaSnpzt3AWm1WhQWFjb7WFe6gnR5c+bM\nwauvvip1DI/B7ela3J6u094/rFvdrRQdHe0cNCsrC1FRUc5p4eHh2LlzJ6qqqmCxWKBWqxEVFYWs\nrCzY7XYcOHAAAwcORFRUlHOM3NxcREZGwm63Y/r06VixYgWCg4MBAJGRkcjLy2v2sYiI6NpptRyG\nDRsGb29v6PV6KJVKREREID09HQAwd+5cLFiwAGPGjMH8+fMBADNnzkRGRgbi4uKQlpYGhUKBpKQk\nFBQUQK/XY8SIEQgNDcWGDRuwd+9ePPPMMzAYDNi1axeCg4MRGxuL2NhYHDhwAJMmTerYtSciomYJ\n7nYN6ezsbO5WcqFPPvkEKSkpUsfwGNyersXt6Trt/d3plqfPcFRXSR3BY8TExEgdwaNwe7oWt6d0\n3LIcqvO+ljoCEZFHc8tysH7d+ofjiIjoyrllOThM51B34r9SxyAi8lhuWQ5+d0yC9evPpY5BROSx\n3LMcEpNw3vgfOGrOSx2FiMgjuWU5KIJC4XXTLag2ZkodhYjII7llOQCA+s7JqNrKXUtERB3BbcvB\nZ3g07GUlqPv5uNRRiIg8jtuWgyBXwG/MRL56ICLqAG5bDgDgN2YiqvO+hqOmRuooREQexa3LQRHS\nHV4DBuP8tiypoxAReRS3LgcA8Bs7GVbuWiIicim3LwffyJGwny3mG9NERC7k9uUgKBTwu3MyrJs/\nljoKEZHHcPtyAAD1HZNQnZ8Fh8UsdRQiIo/gEeUg7xYE38gRqMr6QuooREQewSPKAQDUE+6BdcsG\niA6H1FGIiNyex5SD18AhEFRq1Hy3Q+ooRERuz2PKQRAEqJPugfULvjFNRHS1PKYcAEClH4O6Hw6j\nvuhXqaMQEbk1jyoHmbcP/BKTYd2yQeooRERuzaPKAQDU41NQnb2FFwIiIroKHlcOitAeDRcCyv1K\n6ihERG7L48oBAPyT7oF188cQRVHqKEREbskjy8F76O0Q7XbUHtgjdRQiIrfkkeUgCAL877oPls8z\npI5CROSWPLIcAMBv1J2oO/Ff1P/ys9RRiIjcjseWg+DlDfX4u2HZyFcPRETt5bHlAADq8VNQ/W02\n7BXlUkchInIrHl0Ocm0AVLGjYf3yE6mjEBG5FY8uBwDwn3QfrF9+CkdtjdRRiIjchseXg/L6PvC6\ncTCqc76UOgoRkdvw+HIA0HBY68YMXuuBiKiNukQ5eA8ZDsHbBzV7t0sdhYjILXSJcmj4UNz9sHz+\nodRRiIjcQpcoBwBQxSbCVvwran84LHUUIqJOr8uUg6BQwH/ydFg+fl/qKEREnV6XKQcA8BszCbVH\nD6L+l5+kjkJE1Kl1qXKQ+fjAf+K9MG94X+ooRESdWpcqB6DhlBo1e7fBdqZQ6ihERJ1WlysHmZ8a\nfmMnw/zpWqmjEBF1Wl2uHADAf+JUnDdmwn6uTOooRESdUpcsB7muG1SGcbwYEBFRC7pkOQCA/+T7\nUZX5b9gtlVJHISLqdLpsOSiCu8M3Oh7WLz6SOgoRUafTZcsBAPynPADrFx/DUV0ldRQiok6lTeUw\ne/Zs6PV6zJo1q9H9xcXFSEhIwMiRI5GdnQ0AsFqtSE5ORmxsLNaubTgiyG63IzU1FXq9HsuWLXMu\nO3z4cKhUKjguOVuqTqeDwWCAwWCAyWRyyUq2RBnWGz63RvHVAxHR77RaDvv370dVVRWMRiNqa2ux\nb98+57QlS5Zg0aJFyMzMxMKFCwEAq1evxtSpU2E0GrFmzRrYbDZs2rQJgwYNgtFoRH5+PkpKShAY\nGIicnBxERUU1erzw8HDk5OQgJycHOp3OxavblGbqTFg2roOjytrhj0VE5C5aLYedO3ciMTERADB6\n9Gjs2LHDOa2goABRUVFQqVTQaDSwWCzO+QVBwNChQ3H06NFGY4waNQq7d++Gl5cXtFotRFFs9HhH\njhxBXFwc5s2b58r1bJGyZx/4RMTAsun/rsnjERG5g1bLwWQyQaPRAAC0Wm2jXT2X7g66OO3S+TUa\nTZP7fj+GIAiNHu/EiRPIy8uDyWTC5s2br2LV2k5z70xYN62Hw2q5Jo9HRNTZKVqbQavVwmw2AwDM\nZnOjXT0y2W/dUllZiYCAAOh0OpjNZgQFBcFsNje67+IY/fv3b/HxLo4/ceJEfP/995gwYUKTeebM\nmeO8HR0djZiYmNZW4/IEBcSbI1D04duQJU+7urHcjMViQVFRkdQxPAa3p2txe1657du3N9rTM27c\nuHYt32o5REdHY9WqVUhJSUFWVhZmzJjhnBYeHo6dO3diyJAhsFgsUKvViIqKQlZWFqZMmYIDBw5g\n4MCBiIqKQnZ2NiIiIpCbm4tp0377BSyKonPXUnV1NXx8fCCTybBt2zaEh4c3m+nVV19t10q2hS3t\nKZyd/QBC73sYcn+ty8fvrIqKitCjRw+pY3gMbk/X4va8cikpKUhJSXF+ffGgobZqdbfSsGHD4O3t\nDb1eD6VSiYiICKSnpwMA5s6diwULFmDMmDGYP38+AGDmzJnIyMhAXFwc0tLSoFAokJSUhIKCAuj1\nesTExCA0NBQ2mw2JiYk4ePAgxo4diz179uD48eOIjIxEfHw8Tp8+3WjFOpriup7wjYqHlZ+aJiKC\nIP7+HeFOLjs7GwkJCR0ytu1sEc4+PR3dV30Kuabjj5TqDPiXmWtxe7oWt6frtPd3Z5f+ENzvKUJ7\nwHdkAiyf8VrTRNS1sRx+R3NPGqq2fg676ZzUUYiIJMNy+B1FSHeoDONg/r93pI5CRCQZlkMzNH9M\nQ3XeVtiKT0sdhYhIEiyHZsi1AVAnT0Xl2pVSRyEikgTLoQX+k6ah5uBe1J34r9RRiIiuOZZDC2S+\nKmjunYnKD96UOgoR0TXHcrgM9di7YDtTiJr9u6SOQkR0TbEcLkNQKKCd/hhM778J8ZKTDBIReTqW\nQyt8R44GAJzf1r7zkhARuTOWQysEmQy6GU+h8oN/QqyvlzoOEdE1wXJoA5+ht0HRoxesWzZIHYWI\n6JpgObSRbub/g/mjd2Gv7NjrWhMRdQYshzZS9uoLlX4MzOveljoKEVGHYzm0g+a+R1Bt/A/qT/0o\ndRQiog7FcmgHuUYHzb0PoWL1a3Czy2AQEbULy6Gd1OOnwF56BjV7tkkdhYiow7Ac2klQKKCbOQum\nNa9BtNmkjkNE1CFYDlfAJ2IEFKFhsG7+WOooREQdguVwBQRBgO7hWRcOba2QOg4RkcuxHK6Qsldf\nqEbdicr33pA6ChGRy7EcroL2/v9BzXc7UHvkgNRRiIhciuVwFWQqNbRpT6PirSUQ7Xxzmog8B8vh\nKqni7oBMo4N1M8+7RESeg+VwlQRBQMBjz8L8f2tgLy+VOg4RkUuwHFxAeX0f+I2dDNM7r0sdhYjI\nJVgOLqL5Yxpqjx5EzYE9UkchIrpqLAcXkfn4IuCROQ1vTtfXSR2HiOiqsBxcyCcqDooevWDe8L7U\nUYiIrgrLwYUEQUDAE3+G9YuPeVpvInJrLAcXUwSFQjv9UZz7x0KIdrvUcYiIrgjLoQP4jZ0MQaHk\nNaeJyG2xHDqAIJMhIH0BzOtXw3a2SOo4RETtxnLoIMqw3vC/635UvLmYV40jIrfDcuhA/pOnw246\nh+qcLVJHISJqF5ZDBxIUCnRLfx6md/8Be0W51HGIiNqM5dDBvPoPgt+YiTj3xiLuXiIit8FyuAa0\n0x6GvaQY1VmbpY5CRNQmLIdrQFB6oducl2B69++wlRRLHYeIqFUsh2vE6w/94X/X/Tj3t79CdDik\njkNEdFksh2vI/+7pEOvqYN38sdRRiIgui+VwDQlyObrN/gvM61ej/teTUschImoRy+EaU4b1gua+\nR3HutRd43Wki6rRYDhJQj0+BTK2Fed0aqaMQETWL5SABQRDQbdaLsGZuRM3BvVLHISJqguUgEXm3\nIHR7+gWce/VF2M0mqeMQETXCcpCQb0QMVLGJOPf6S/z0NBF1KiwHiWkfeAL28lJe+4GIOpU2lcPs\n2bOh1+sxa9asRvcXFxcjISEBI0eORHZ2NgDAarUiOTkZsbGxWLt2LQDAbrcjNTUVer0ey5Ytcy47\nfPhwqFQqOC75UNiKFSsQGxuL6dOnw94FrqQmKJUIfGYRzBmrUPfzcanjEBEBaEM57N+/H1VVVTAa\njaitrcW+ffuc05YsWYJFixYhMzMTCxcuBACsXr0aU6dOhdFoxJo1a2Cz2bBp0yYMGjQIRqMR+fn5\nKCkpQWBgIHJychAVFeUcr7S0FN988w3y8/MRHh6OjRs3dsAqdz7KsF7QPTwL5Uvnw1FzXuo4RESt\nl8POnTuRmJgIABg9ejR27NjhnFZQUICoqCioVCpoNBpYLBbn/IIgYOjQoTh69GijMUaNGoXdu3fD\ny8sLWq220b72vXv3Ij4+HgCQkJDQ6LE8nZ9hPLxuvIkXByKiTqHVcjCZTNBoNAAArVYLk+m3I2su\n3R10cdql82s0mib3/X6Mtj5WVxDw+DzU/3wcVV99KnUUIuriFK3NoNVqYTabAQBmsxk6nc45TSb7\nrVsqKysREBAAnU4Hs9mMoKAgmM3mRvddHKN///7O5QRBaPRYhYWFzT7WpebMmeO8HR0djZiYmDat\nrDsQZ/4JFUufRaU2GMIf+re+wFWyWCwoKuJ1rl2F29O1uD2v3Pbt2xvtfRk3blz7BhBb8d1334mP\nPvqoKIqi+Pjjj4t79uxxTnv66afFHTt2iFarVRw1apQoiqL42muvievXrxdtNpsYFxcn1tfXi599\n9pm4ZMkSURRFMSkpSTxz5oxzjPj4eNFms4miKIolJSXihAkTRFEUxWXLlokbNmxokicrK6u1yG6v\naluOWPjgBNFWWdHhj1VYWNjhj9GVcHu6Fren67T3d2eru5WGDRsGb29v6PV6KJVKREREID09HQAw\nd+5cLFiwAGPGjMH8+fMBADNnzkRGRgbi4uKQlpYGhUKBpKQkFBQUQK/XIyYmBqGhobDZbEhMTMTB\ngwcxduxY7NmzB8HBwYiNjUVsbCwOHDiASZMmta/pPIQqZhRUsaNxbvnzELvAEVtE1PkIouhe735m\nZ2cjISFB6hgdTrTbUDr/cXjfEgHttEc67HGKiorQo0ePDhu/q+H2dC1uT9dp7+9OfgiukxLkCgQ+\nuxhVWz/H+d3fSh2HiLoYlkMnJu8WhMB5S3Hu9b+i/pefpY5DRF0Iy6GT8x4UDu2Mp1D28mw4LGap\n4xBRF8FycAPqxGT43BaLsqXzeIEgIromWA5uQpeWDkEQYHrn71JHIaIugOXgJgS5AoHPLEbNnm2w\nZv5b6jhE5OFYDm5E5q9B0AuvofL9N1F7aL/UcYjIg7Ec3Izy+j4InPsyyl75M+pPn5Q6DhF5KJaD\nG/IZFgXtg0+g9MWnYTedkzoOEXkgloObUicmwy/+TpT9dRYcNTVSxyEiD8NycGOa+/8Hip69cW7F\nczwHExG5FMvBjQmCgG7pz8NRZYXpndeljkNEHoTl4OYEpRJBC5ajZv8umD9dK3UcIvIQLAcPIFP7\nI/jlN2Dd/DGs/9kkdRwi8gAsBw+hCApF8MI3UfnBP1G9PVfqOETk5lgOHkQZ1hvBf/k7Kt5cjJoD\ne6SOQ0RujOXgYbxuGIjAeUtQvnQ+6o4fkToOEbkploMH8hkyHN3Sn0PpX2eh/pefpI5DRG6I5eCh\nfKPioEt7GqXPPcHTbBBRu7EcPJifYRw00x9D6fzHUV/4i9RxiMiNsBw8nDoxGZppD6N0wWOwFZ+W\nOg4RuQmWQxegHnsXNFNmoGT+Y7CdLZI6DhG5AZZDF6EenwL/u+5DybzHYCspljoOEXVyLIcuxD/5\nXvhPvBclzz6C+qJfpY5DRJ0Yy6GL8Z84FZp7ZqD0z//Dw1yJqEUshy5IfedkaB98EiXzH4PIgiCi\nZrAcuig/wzgEPDoXjtdfRO2xQ1LHIaJOhuXQhalGjobswXSU/XUWag7ulToOEXUiLIcuTgiPROCz\nr6B8yTxUb8uWOg4RdRIsB4LPLREIfvkNmN5eAeuWT6SOQ0SdAMuBAABe/QYiZOlqWDZmoHLt/0IU\nRakjEZGEWA7kpLiuJ0KWv4Pz+7aj4o3FEO02qSMRkURYDtSIXNcNIa/8L+ylxSh7aQ4c1VVSRyIi\nCbAcqAmZrwpBL74OeVAISp6ZCVvJGakjEdE1xnKgZgkKBQKenA+/hPEo+dMM1P5wWOpIRHQNsRyo\nRYIgwP+u+xHw2LMoe/FpVG/LkToSEV0jCqkDUOfnGx2P4ODuKHt5Dmy//gz/P6ZBEASpYxFRB+Ir\nB2oTrxsGIvS193F+dz7KX3kWjvPVUkciog7EcqA2kwcGI2TpKsj8/HF2zoM87TeRB2M5ULsISi8E\npD8H//H3oGTuQzi/d7vUkYioA/A9B2o3QRCgHp8CZZ8bUL50HurG3Q3NPWkQZPxbg8hT8KeZrpj3\n4KEIee0D1H63C2UvpsNeWSF1JCJyEZYDXRVFUAiCX1kJZb+BOJt+P2oPfy91JCJyAZYDXTVBroDu\nwScR8OQ8lC1+BuYN70N0OKSORURXgeVALuMbORKhf/sXzu80ouzFp2E/VyZ1JCK6QiwHcilFSHeE\nLF0FrxsH48xT9+H8zjypIxHRFeDRSuRygkIB7fRH4XNrFMpffQHn922H7qFZkPn4SB2NiNqIrxyo\nw3gPHorub6yDWF2Ns0/fj7oT/5U6EhG1UZvKYfbs2dDr9Zg1a1aj+4uLi5GQkICRI0ciO7vh+sNW\nqxXJycmIjY3F2rVrAQB2ux2pqanQ6/VYtmzZZcfV6XQwGAwwGAwwmUxXvYIkLZmfGoFzX4Zm6kyU\nvvAUKteuhFhfJ3UsImpFq+Wwf/9+VFVVwWg0ora2Fvv27XNOW7JkCRYtWoTMzEwsXLgQALB69WpM\nnToVRqMRa9asgc1mw6ZNmzBo0CAYjUbk5+ejpKSkxXGHDBmCnJwc5OTkQKfTddBq07XmFz8W3d9c\nj7qfj+Ps09NRd/yI1JGI6DJaLYedO3ciMTERADB69Gjs2LHDOa2goABRUVFQqVTQaDSwWCzO+QVB\nwNChQ3H06NFGYxgMBuzatavFcY8ePYq4uDjMmzfP5StL0pJ3C0LQ86/C/54ZKP3L/4Ppg3/yVQRR\nJ9VqOZhMJmg0GgCAVqtttKvHccmx7BenXTq/RqNp032XjnvixAnk5eXBZDJh8+bNLlpN6iwEQWh4\nFfHGOth+/Rln0u/nhYSIOqFWj1bSarUwm80AALPZ3GhXj+ySc+lUVlYiICAAOp0OZrMZQUFBMJvN\nje67OEb//v1htVqbHffi/xMnTsT333+PCRMmNMk0Z84c5+3o6GjExMS0e8WpgcViQVFRkSSPLc6Y\nBXFPPkpeSIdwexyEidMg+KgkyeIqUm5PT8TteeW2b9/eaE/PuHHj2rV8q+UQHR2NVatWISUlBVlZ\nWZgxY4ZzWnh4OHbu3IkhQ4bAYrFArVYjKioKWVlZmDJlCg4cOICBAwciKioK2dnZiIiIQG5uLqZN\nm4awsLAm41ZXV8PHxwcymQzbtm1DeHh4s5leffXVdq0ktayoqAg9evSQLkDYvbCPugOV7/0DNX9J\nh+7h2fAdmeC2FxOSfHt6GG7PK5eSkoKUlBTn1xcPGmqrVncrDRs2DN7e3tDr9VAqlYiIiEB6ejoA\nYO7cuViwYAHGjBmD+fPnAwBmzpyJjIwMxMXFIS0tDQqFAklJSSgoKIBer0dMTAxCQ0MbjatQKBAR\nEYHjx48jMjIS8fHxOH36dKMVI88l1wag2/97EYHPLELl+lUoeyGd14ogkpggiqIodYj2yM7ORkJC\ngtQxPEZn+8tMtNlg+fc6WD75AOqke6FJSYXg5S11rDbrbNvT3XF7uk57f3fyQ3DUqQgKBTR3pyL0\n7xmo/+kYzjxxL6p3fAM3+xuGyO2xHKhTUoR0R9BzK6B79BlU/ustlM5/DHU/HpM6FlGXwXMrUafm\nOzwaPkMjUbV1I0pfSIfvbSOhnf4Y5N2CpI5G5NH4yoE6PUGugHp8Cq57+xPI1BqceeKPMH/0Lhy1\nNVJHI/JYLAdyGzK1P3QPPY3Q1z5A3YmjOPPI3bBu/RyizSZ1NCKPw3Igt6O4rieCFixH4LylqDZm\n4sxjU1D1zVZefY7IhVgO5La8B96MkMUrEfDEPFj/vR5n0+/D+V1GHtlE5AJ8Q5rcns/Q2+B9SyRq\nduah8oN/wvzx+9De9wi8h93utp+0JpIay4E8giAI8I2Oh89tsajO/w8qVr0Kma8vNH98CD63xUKQ\n8UUyUXvwJ4Y8iiCXN5z19a2P4J/yACoz3sbZp6ah2pgJ0W6XOh6R22A5kEcSZDKoRiQg9B8Z0D7w\nJCz/Xt/wxvV/voBYXy91PKJOj7uVyKMJggDf20bCJ3IEag/uhfmj91C59i2oJ9wDvzsnQ+6vlToi\nUafEcqAuQRAE+NwSCZ9bIlH30w+wbMxA8UOT4DdqLNQTp0HZ43qpIxJ1KiwH6nK8+t6IwNl/hb28\nFJbNH6Nkzgx4Dx4G9aSp8B48jEc4EYHlQF2YPDAYugeegOaPaajK+gIVbyyCIFfAb1wK/Ax3QqZS\nSx2RSDIsB+ryZD6+8J9wD9Tjp6D24F5Yt2yAee1K+OrHQD0+BV59bpA6ItE1x3IguuDS9yVsZSWo\nytyIsheegrx7GNRj74JvTAJkPj5SxyS6JngoK1EzFEEh0E57BNe9+wX8J05Fdd7XKH5gHM69uRi1\nxw7xFB3k8fjKgegyBIUCqhEJUI1IgK3sLKqzt+Dc8ucApRf8EpPhZxgHua6b1DGJXI7lQNRGiqBQ\naP6YBv97ZqD28H5U/ecLFD8yGd6DwqGKGwvfqDipIxK5DMuBqJ0EQYDPzbfC5+Zb4Xh0Ls7vMqL6\nm62oWLkUuGkYzo+dBJ+IGAhKL6mjEl0xlgPRVZD5quAXPxZ+8WNhN5tQvOUzWDauw7m/vwzf6Hio\n4sfC++ZbIcjlUkclaheWA5GLyDU6yOLGImRqGmylZ1Bt/A9M77wOe1kJfKPiGs4aO/Q2vqIgt8By\nIOoAiuDu0Nw9HZq7p8NWfBrnd34D84b3Ub78OfgMj4YqxgCf4TGQqfykjkrULJYDUQdTXNcT/nfd\nD/+77oe9ohznd+ah6j9f4NzfF8L75mHwjY6Hb8QIyAODpY5K5MRyILqG5AGBUN85Geo7J8NRZcX5\nPd/i/K48VL73BuRBofCJHAHf4THwGjQEgpw/niQdPvuIJCLzUzvfzBbtNtQdO4Tze7ahYtUK2M8W\nw3vY7fCNGAGf4dGQBwRKHZe6GJYDUScgyBXwvmkovG8aCjzwBOzlpTi/bzvO785HxaoVUIRcB+/w\niIbrZd98K9+roA7HciDqhOSBwVCPmQj1mIkNryqOH0Xt97th+TwD5UvnQ9mnP3yGRsL7lkh4Dwrn\nEVDkciwHok5OkCvgPXAIvAcOgebeh+CorUHdkQOoObAHle/9A/W/noTXjYPhPXgovAcPg9fAIZD5\n+Eodm9wcy4HIzci8feAz7Hb4DLsdAOCwmFF79ABqD3+Pyg//F/U//QBl777wGjwM3oOHwfumoZBr\ndRKnJnfDciByczJ/DXxvi4XvbbEA0PDK4vgR1B3aj6qvPsO5v/0F8oAgeA24ueEVxoCboexzAwSl\nUuLk1JmxHIg8jMzbx3nuJwAQ7XbUn/oRdccOoe6Hw6j68hPYzhRC+Ycb4TVgsLMw5N3DeIlUcmI5\nEHk4QS6HV98b4dX3RuDOyQAAR3UV6k4cbTh8dls2Kt/9B8T6OnjdeBOUfQfAq98AKPsOgKJ7GAQZ\nL/vSFbEciLogmcoPPuER8AmPcN5nKytB3fEjqP/pGKpyvkT9mtfhsFqg7NsfXn0HQNn3xobS6NWX\nR0d1ASwHIgLQcPU7RVAIEB3vvM9uNqH+5+Oo+/EYag/uhWVjBuxnCqHo0QvK3v2g7NW34f/e/SAP\n7cFXGR7BP1lvAAAQPklEQVSE5UBELZJrdJBfuK72RY7aGth++Qn1J39E/S8/wfrVp6g/9RMcZhMU\nPftA2edCafTqB2XvvpAHd2dpuCGWAxG1i8zbB179b4JX/5sa3e+otqL+l59Rf6qhNGq+3w3bLz/D\nYa2EontPKHr2hiKsN5Q9ekER1guKsN48xLYTYzkQkUvIVGrnh/Uu5ThfDVvRL7AV/oL6wl9Qc2A3\nbF9+gvrCUxBk8gtF0ctZGvLQHlB0D4PMXyvRmhDAciCiDibzVcGr30B49RvY6H5RFOGorLhQGqdg\nK/wF1cZM2M4WwXamEHCIEAODUXZ9H2dhKELDoOjeA/LQHpB5+0i0Rl0Dy4GIJCEIAuS6bpDrusF7\n8NAm0x0WM4oPfQ+VvQ62s0UNu6p2fwvb2ULYSs5ApvaHIuQ6yINCIQ8OhSIo1HlbHhQKeUAgL896\nFVgORNQpyfw1EHrfAFWPHk2miQ4H7OdKYS85A3vZWdhKz8JWUozaw9/DVnYW9tKzcFgrIQ8Ialoc\ngSGQdwtq+BcQCMHLW4K16/xYDkTkdgSZDIqghl/6LRHr62AvK4GtvAT20obCsJ0+hZrvd8NRUQZ7\nRTnsFeWQ+aggu1AU8oALpdEtELKAoEu+DoKg8utSnyBnORCRRxKUXlBc1xOK63q2OI/ocMBhMcNe\nUdZQGOca/tlKz8J+7DAcFeUN91WUQayvh1wbAJlWB5k24MLt3/5vuH1hmiYAgtrfrcuE5UBEXZYg\nk0Gu1TUcUtvnhsvO66itgaPSBEdlBeyVFXCYKy7cNsF2prDhvsoK53SxrhYyfy3kum6Q+WshU/v/\n9r9aC5m/BjK1puH/S24Lvp3jFQrLgYioDWTePpCFdAdCurdpfrG+DvYLZeKwmuGwVMJhMV+4bYbt\nzGnnbef9VjPE2toLRXJJeag1EFRqyPzUkKnUEFR+Dbd//7VKDUGlcsn1x1kOREQdQFB6NZyOJCik\nXcuJ9fVwVFkulEYlHFYLHNZKOKqsEKusDSVSUtRwu7oKjuoL91c3fC2er4bg5Q2Zn7qhUFR+CH7p\njXbnZzkQEXUiglLpPMT3SogOB8Sa6oaiqLLCUWWF4Nv+KwO26YQns2fPhl6vx6xZsxrdX1xcjISE\nBIwcORLZ2dkAAKvViuTkZMTGxmLt2rUAALvdjtTUVOj1eixbtuyy465YsQKxsbGYPn067HZ7u1eI\n2mf79u1SR/Ao3J6uxe3ZfoJMBplKDUVQKJS9+8H7pluuaDdTq+Wwf/9+VFVVwWg0ora2Fvv27XNO\nW7JkCRYtWoTMzEwsXLgQALB69WpMnToVRqMRa9asgc1mw6ZNmzBo0CAYjUbk5+ejpKSk2XFLS0vx\nzTffID8/H+Hh4di4cWO7V4jaZ8eOHVJH8Cjcnq7F7SmdVsth586dSExMBACMHj260TeroKAAUVFR\nUKlU0Gg0sFgszvkFQcDQoUNx9OjRRmMYDAbs2rWrybjbt2/H3r17ER8fDwBISEjgE4OISCKtloPJ\nZIJGowEAaLVamEwm5zSHw+G8fXHapfNrNJo23afValFZWXnZxyIiomun1R1RWq0WZrMZAGA2m6HT\n/XaKXdkl52ivrKxEQEAAdDodzGYzgoKCYDabG913cYz+/fvDarU2GVen0+H06dPNPtalLr6/QVdv\n3Lhx3J4uxO3pWtye0mm1HKKjo7Fq1SqkpKQgKysLM2bMcE4LDw/Hzp07MWTIEFgsFqjVakRFRSEr\nKwtTpkzBgQMHMHDgQERFRSE7OxsRERHIzc3FtGnTEBYW1mTc3r1746233sLcuXORlZWFqKioJnkS\nEhJcuwWIiKiJVncrDRs2DN7e3tDr9VAqlYiIiEB6ejoAYO7cuViwYAHGjBmD+fPnAwBmzpyJjIwM\nxMXFIS0tDQqFAklJSSgoKIBer0dMTAxCQ0MbjatQKBAREYHg4GDExsYiNjYWBw4cwKRJkzp27YmI\nqFmCKIqi1CGIiKhzcasLu7b0eQtqv1OnTqF79+4wGAwYO3as1HHcVnFxMYYPHw6VSuU8QIOf1bky\nzW1LnU4Hg8EAg8HAA1Taaffu3RgxYgT0ej3mzJkDAFi+fHmbn5tuUw6X+7wFXZkxY8YgJycHW7du\nlTqK2woMDEROTo7z/TF+VufK/X5bAsCQIUOQk5ODnJycFg9Qoeb16dMHubm5MBqNKCkpgdFoRF5e\nXpufm25TDpf7vAVdmZycHMTFxeH111+XOorb8vLyglb727WO+VmdK3dxW166p/vo0aOIi4vDvHnz\nJEzmnkJCQuDl5QUAUCgUOHLkSLuem25TDvwMhGv16NEDx48fR25uLrKzs3Ho0CGpI3kEPk+v3qWn\nqz5x4gTy8vJgMpmwefNmCVO5r4MHD6KsrAw6na5dz023KYfLfd6C2k+pVMLX1xcymQzjx49nObgI\nn6eudXH7TZw4kc/RK1BRUYH09HS8++670Gg07Xpuuk05REdHOz8M09JnIKjtrFar8/a2bdvQr18/\nCdO4v4u7QiIjI5GXlweAz9MrJYoiRFFEdXW1841pPkfbz2634/7778eKFSsQHBzc7uem25RDc5+L\noCuXn5+PiIgIjBw5Ej179kRkZKTUkdySzWZDYmIiDh48iDvuuAMnT56EXq/nZ3WuwKXbcuzYsTh0\n6BAiIyMRHx+P06dPIyUlReqIbmXDhg3Yu3cvnnnmGRgMBvz000/tem7ycw5ERNSE27xyICKia4fl\nQERETbAciIioCZYDERE1wXIgIqImWA5ERNQEy4GIiJpgORARURMsByJqVn5+PtauXYsTJ064fOzP\nPvvM5WOSa7EcqFl5eXno06cPDAYDYmNjcezYMeTl5TU6edeMGTPw008/wWKxYMKECRg1ahRiYmLw\n3XffNRkvOzsbo0aNQnx8PO6++25UVFS0Ocfzzz/v0nVzhYvbJyEhAYmJiSgvL292vlOnTiE3N7dD\nsyxfvhznz593+bhnzpzBH/7wBxQXF7t03F9//RU5OTnOr7Ozs/Htt9+69DHo6rEcqEWpqanIycnB\n8uXLsXLlSgBAr169sGbNGgC/nVr5X//6F+6++27k5ubi22+/xYABAxqNU1ZWhpdffhlbtmzBN998\ng6VLl6Kurq7NOS49hXNnkpqaiuzsbDz44INYt25ds/OcPHmy0S/CtioqKsLIkSPx8ccf48MPP8SS\nJUtanDcsLAy+vr7Iy8vDjTfeiNLSUue0Bx98EM8++yxsNlu7M0yePBlHjhxxnrDNVWNff/316NGj\nh/PrhIQEnDp1qt35qGOxHKhFF0+7ZTKZoNVqIQgCkpOT8cUXX8DhcDinq1Qq7NixA+Xl5ZDJZPDz\n82s0zpdffonU1FSoVCoAwA033ICgoCBMmzYNcXFxuO++++BwOFBcXAyDwQC9Xo8nn3yySZ7HH38c\no0ePRlJSEiorK7Fjxw5ERUUhISEB7733HkRRxMMPP4z4+HiMHz++xeXy8vIwbtw4JCcnIzY2FtXV\n1c0u+/vlLrd9fHx8ms2/atUqrF271nmhqtbGvKhHjx4YMGAA7rnnHowfP/6yVz68mCMuLg4zZszA\nxx9/DKChlAMCAjB+/HgoFIoWlweAwsJC7Nq1C7t27UJBQYEzu16vx7vvvuvysX9/SrfO+gdAV3b5\n7yp1aWvXroXRaMSJEyeQmZmJsrIyyOVyJCUl4bPPPnP+QKempuL06dMYNWoUunfvjg8//BAhISHO\ncYqLixEeHt5o7M8//xyDBw/GunXrsHjxYnz66ae46667kJWVBZlMhunTp+PHH390zr9582b07t0b\nb731FrZu3YqVK1fi/PnzWLZsGfR6PQBg48aNCA0NxerVqy+7XHR0NLy9vfH555/jlVdeQXZ2Nux2\ne6Nlm1vuz3/+c5Pts3XrVlRWVuLQoUMQBKFJ/kceeQT9+vXDSy+91KYxLyorK4NCocDmzZuxf/9+\nZGRktPh9uvQXa3R0NFavXo0nnngCJ0+exPXXXw8AeO+995Cbm4tp06bh8OHDSEpKQnZ2NsrLy/Hw\nww8jLCwMYWFhjcYdMGCA8zrElxv7jTfewH//+1+88MILmDdvHkaMGIEffvgBSqUSCxcubDL26dOn\ncejQIRQWFjZ5TOo8WA7UotTUVLz00ksoLS3FQw89hD/96U8QBAEzZ87ElClTnD/Ycrkczz//PJ5/\n/nl89NFH+Nvf/oZXXnnFOc51112HwsLCRmP/+OOPuPXWWwEAw4cPx3fffYeysjI89thjMJlMOHXq\nFIqKigA0/JV59OhRrF+/Hl9//TVsNhuio6Mxa9YsvPzyy1izZg2eeuop/PDDD4iJiWn0OM0tBwA3\n33wzgIa/0CsqKnDmzJlGy7a0XHPb54EHHkBubi5uueWWZvO3Z8yLtm/fjkmTJuHOO+/Ehg0bnJd7\nbI0gCOjTpw9OnTrV6K9zg8GAX3/9FWPHjsWWLVswa9YsbNy4ERUVFQgMDGx2LIPBcNmxL47/5JNP\n4pFHHkFtbS1eeukl2O121NbW4vHHH2923J49e2L9+vVtWh+SDncrUYsu/vD7+fnBYrE479doNBgw\nYAB27doFAPjll1+c+52Dg4Ob7DIYN24cMjIynBcY+vHHHxEWFoa9e/cCaLjucr9+/bBu3Trcdddd\nyM3NRUxMTKNxBg4ciAceeAA5OTkwGo1YvHgxdDod/vnPf2Lp0qV48cUXMXDgQOd1cS8uO2DAgEbL\nLVq0CEDT3Ri/X/b3yy1evLjF7bRgwQKsXLkS69evb5JfqVQ6t01zY9rtdpSUlDQZc/v27Rg2bBgA\n4NixYwDQ4ry/d++992LJkiXo06dPk+8FAMhkMgiCAFEUERISgpqamlbHbG7siwRBQLdu3bB//370\n7NkTAJrsWiT3w3KgFn344YcwGAwYPXo05s6d2+gXTXp6uvOX1vfff4+RI0fCYDBg6dKlSE9PbzRO\nUFAQnn/+eUyYMAHx8fGYO3cu7rjjDhw+fBhxcXE4dOgQ7r77bhgMBqxYsQKTJ09GdXW1c3lBEJCU\nlISff/4ZCQkJGD16NL766iu8/fbbiIuLQ1JSEtLS0pCUlITi4mLnfQCQnJzcaLmtW7c2u67JycmN\nlv39cl999VWL2+nGG2+E1WpFTExMk/w333wztm3bhqlTpzY75smTJ/Hcc881Gm/v3r3Iyspyvtrq\n378/Pv/882bnvVRmZibeeecd9OrVC2q1GgCwY8cOfPnll6irq8ORI0eQmZmJAQMGYPny5XjrrbeQ\nmZnZ4niXG1sQBOfYNpsNKSkpzldK27Ztg9FobPQ9JPfDi/0QSeizzz5Dt27dEB8ff8Xzrlu3DtOm\nTbvssqdOncL777+PF1988SrSNq+8vBzl5eXQarUIDQ29ojEyMjJw3333uTgZXQ2+50AkocmTJ1/1\nvG35+2779u3Yv38/6urq2vz+RVu99dZbuP322zFmzJgrHoNHK3U+fOVA5OZWrFiBJ554Ar6+vlJH\nuSLZ2dnw8fFxHhVFnQPLgYiImuAb0kRE1ATLgYiImmA5EBFREywHIiJqguVARERNsByIiKgJlgMR\nETXBciAioib+P4cMyU8PivqwAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "plt.plot(o1scenario.bns_prior.rates.to(u.megaparsec**(-3)/u.megayear), np.exp(o1scenario.bns_prior.pdf_data))\n", + "plt.xlabel(r'BNS Coalescence Rate, $R$, ${\\rm [Mpc^{-3} Myr^{-1}]}$')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "?np.percentile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "np.percentile(posterior.theta_samples, 99.99)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "posterior.theta_samples" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "GRBeaming (Python 2)", + "language": "python", + "name": "grbeaming" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + }, + "latex_envs": { + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 0 + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/tests/__init__.pyc b/tests/__init__.pyc index 07a4f59..a87de62 100644 Binary files a/tests/__init__.pyc and b/tests/__init__.pyc differ diff --git a/tests/test_distributions.pyc b/tests/test_distributions.pyc index 5335a9b..bb39a19 100644 Binary files a/tests/test_distributions.pyc and b/tests/test_distributions.pyc differ diff --git a/tests/test_scenarios.py b/tests/test_scenarios.py index 23f6c65..85bced9 100644 --- a/tests/test_scenarios.py +++ b/tests/test_scenarios.py @@ -2,10 +2,10 @@ # -*- coding: utf-8 -*- """ -test_distributions +test_scenarios ---------------------------------- -Tests for `distributions` module. +Tests for `scenarios` module. """ import unittest @@ -20,7 +20,7 @@ def setUp(self): self.rateunits = u.megaparsec**-3*u.year**-1 # Create a BNS distribution which is uniform rates = np.linspace(0, 100, 100)*self.rateunits - pdf = np.ones(len(rates)) * (1.0/100) + pdf = np.ones(len(rates)) * (1.0/100) * self.rateunits**-1 self.bnsuniform = scenarios.BNSDistribution(rates, pdf) # Create an observing scenario with the same BNS distribution self.uniformscenario = scenarios.Scenario(self.bnsuniform) diff --git a/tests/test_scenarios.pyc b/tests/test_scenarios.pyc index b39fd14..8163c88 100644 Binary files a/tests/test_scenarios.pyc and b/tests/test_scenarios.pyc differ diff --git a/theta_rates.py b/theta_rates.py deleted file mode 100755 index b1de105..0000000 --- a/theta_rates.py +++ /dev/null @@ -1,151 +0,0 @@ -#!/usr/bin/env python -# -*- coding:utf-8 -*- -# Copyright (C) 2014-2015 James Clark -# -# This program is free software; you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation; either version 2 of the License, or -# (at your option) any later version. -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - -""" - -theta_rates.py - -Monte-Carlo demonstration of impact of jet angle on relative rates - -Procedure for Nbns events -0) Specify jet distribution -1) Draw a jet angle -2) Draw an inclination angle -3) If inclination < 0.5 * jet angle, GRB is seen; Ngrb += 1 - -""" - -from __future__ import division -import os,sys -import numpy as np -import scipy.stats as stats - -import matplotlib -from matplotlib import pyplot as pl - -def truncparms(low,upp,mu,sigma): - a = (low - mu) / sigma - b = (upp - mu) / sigma - return a, b - - -def compute_efficiency(k,N,b=True): - - if b: - # Bayesian treatment - epsilon=(k+1)/(N+2) - stdev_epsilon=1.64*np.sqrt(epsilon*(1-epsilon)/(N+3)) - else: - # Binomial treatment - if N==0: - epsilon=0.0 - stdev_epsilon=0.0 - else: - epsilon=k/N - stdev_epsilon=1.64*np.sqrt(epsilon*(1-epsilon)/N) - return (epsilon,stdev_epsilon) - - -# ------------------------ -# Setup -Nbns = 1e5 - -# ------------------------ -# Distributions - -# -# Inclination angle distribution: uniform in cos(iota) -# -iotas = np.arccos( 0 + 1*np.random.rand(Nbns) ) -# to degrees -iotas *= 180/np.pi - -# -# Jet angle distribution: truncated normal -# -theta_mu=np.arange(5,35,5) -theta_sigma=np.arange(1,16) -theta_low=0.01 -theta_upp=90 - -FracGRB=np.zeros(len(theta_sigma)) -deltaFracGRB=np.zeros(len(theta_sigma)) - -f = pl.subplots() - -linestyles=['-', '--', ':'] -markers=['s', '^', 'o'] - -l=0 -mk=0 -for m,mu, in enumerate(theta_mu): - print m - for s,sigma in enumerate(theta_sigma): - - #a, b = truncparms(theta_low, theta_upp, theta_mu, theta_sigma) - #thetas = stats.truncnorm.rvs(a, b, loc=theta_mu, scale=theta_sigma, size=Nbns) - a, b = truncparms(theta_low, theta_upp, mu, sigma) - thetas = stats.truncnorm.rvs(a, b, loc=mu, scale=sigma, size=Nbns) - - - # ------------------------ - # Count GRBs - Ngrbs=sum(iotas <= thetas) - FracGRB[s], deltaFracGRB[s] = compute_efficiency(Ngrbs,Nbns) - - if l