2026 IPTA Student Week
Author: Bjorn Larsen
Credits to past IPTA Bayesian inference lecturers for inspiration and materials, particularly Anuradha Samajdar (IPTA 2021) and David Wright (IPTA 2025).
Note this tutorial is broken up first into some pen and paper exercises practicing Bayesian statistics and then an in-depth data analysis application covering parameter estimation and model selection techniques. The full tutorial including all exercises may take longer than the alloted time, so budget your time wisely!
# imports for computational exercises
from IPython.display import display, Math
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import tqdm
import emcee
import nautilus
import corner
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) /tmp/ipykernel_3024/1209946719.py in <cell line: 0>() 5 import numpy as np 6 import tqdm ----> 7 import emcee 8 import nautilus 9 import corner ModuleNotFoundError: No module named 'emcee' --------------------------------------------------------------------------- NOTE: If your import is failing due to a missing package, you can manually install dependencies using either !pip or !apt. To view examples of installing some common dependencies, click the "Open Examples" button below. ---------------------------------------------------------------------------
Definitions:
All probabilities must be normalized. If $X$ is a discrete variable, then the normalization condition for the probability $P(X)$ is
\begin{align} 1 = \sum_XP(X). \end{align}If $x$ is a continuous variable then the normalization condition for the probability density $p(x)$ is
\begin{align} 1 = \int p(x)dx. \end{align}Marginalization amounts to a reduction in dimensionality via integration over one or more variables. For example, for two continuous variables $A$ and $B$, one may recover $P(A)$ from the joint probability $P(A,B)$ as such,
\begin{align} P(A) = \sum_B P(A,B), \end{align}or for the continuous case,
\begin{align} p(a) = \int p(a,b)db. \end{align}Usually when we do inference, we have some observations or data $D$ which we want to use to understand the parameters $\theta$ of some underlying model $\mathcal{M}$, or whether an event $X$ has occurred, or how likely our model $\mathcal{M}$ is to explain the observations $D$. For the following definitions we will assume we are attempting to infer the parameters $\theta$ of a model $\mathcal{M}$. Note that the fancy choice of symbols chosen here is a common convention but not universal (oftentimes, inluding in this tutorial, simply $p$ will be used, with the items in the parentheses used to contextualize which probability is being used):
Prior: $\pi(\theta) \equiv$ probability density for the parameters $\theta$, prior to making observations $D$
Posterior: $\mathcal{P}(\theta | D) \equiv$ probability density for the parameters $\theta$, after we have observed $D$
Likelihood: $\mathcal{L}(D | \theta) \equiv$ probability of obtaining the observations $D$ assuming the model parameters $\theta$
Evidence: $\mathcal{Z}(D) \equiv$ probability of obtaining the data $D$ regardless of the value of $X$ is true or false. The "evidence" is better known as the "model evidence" (or a model likelihood) since in practice we are always assuming some model, i.e. $\mathcal{Z}(D) \doteq \mathcal{Z}(D | \mathcal{M})$. This is also sometimes referred to as the "marginal likelihood", since $\mathcal{Z}(D)$ can be obtained via the marginalization formula as \begin{align} \mathcal{Z}(D) = \int \mathcal{L}(D | \theta)\pi(\theta) d\theta. \end{align} For a discrete variable $\theta_i$, the marginalization is \begin{align} \mathcal{Z}(D) = \sum_{i} \mathcal{L}(D | \theta_i)\pi(\theta_i). \end{align} Also I apologize there are so many names for this thing but I didn't come up with that!
Bayes theorem: This is how to compute the posterior probability, starting from the prior, likelihood, and evidence, \begin{align} \mathcal{P}(\theta | D) = \frac{\mathcal{L}(D | \theta)\pi(\theta)}{\mathcal{Z}(D)}. \end{align} Note here that the evidence plays the role of a normalization factor.
Bayes Factor: Bayes Factors are computed for the purpose of comparing two models $\mathcal{M}_1$ and $\mathcal{M}_2$ as the ratio of model evidences \begin{align} \mathcal{B}^{\mathcal{M}_2}_{\mathcal{M}_1} = \frac{\mathcal{Z}(D | \mathcal{M}_2)}{\mathcal{Z}(D | \mathcal{M}_1)} \end{align}
Odds ratio: This is just the Bayes Factor with a prior on the two models under consideration explicitly considered, \begin{align} \mathcal{O}^{\mathcal{M}_2}_{\mathcal{M}_1} = \mathcal{B}^{\mathcal{M}_2}_{\mathcal{M}_1}\frac{\pi(\mathcal{M}_2)}{\pi(\mathcal{M}_1)}, \end{align} Odds ratios may be interpreted as a "betting odds", i.e., what are the odds I should bet on the description of reality offered by model $\mathcal{M}_2$ as opposed to $\mathcal{M}_1$.
These definitions underlie nearly all of Bayesian parameter estimation and model selection.
You can start here with some Bayesian statistics exercises, or skip to the computational problem. These exercises should be helpful to get comfortable with conceptually with priors, posteriors, evidences, and likelihoods if you are not already familiar. You can work on your solutions in the markdown boxes below each problem definition, or work on separate paper. Look at the hints only if you are stuck!
Derive Bayes' theorem for the posterior $p(x | y)$ as a function of the prior, likelihood, and evidence. Start from the joint probability $p(x, y)$. This is just a warm up and should take only 2-3 lines.
Your solution here:
Add blockquote
A classic problem where a seemingly unexpected answer can become very intuitive using Bayesian statistics. This will help give some practice setting up Bayes' theorem, marginalization, and computing posteriors and evidences.
A rare disease $X$ has been going around, right now it is only affecting 1 out of 100,000 people in the population but you decide to take a test for it just in case. The test reports either True ($T$) or False ($F$). The test is not perfect but still highly accurate; it has a 1 in 5000 chance to report a false negative (test returns $F$ when it should be $T$) and a 1 in 1000 chance to report a false positive (test returns $T$ when it should be $F$).
Suppose you take the test and it returns $T$. Oh no! What is the probability that you have contracted the disease $X$? (Use $!X$ to indicate the case of not having the disease). Click the solution box if you want to check your numerical answer, but first ask yourself if your answer makes sense and why.
Your solution here:
##### use and modify this space to perform the numerical calculation
# define variables (priors, likelihoods)
prior_X = 1
prior_not_X = 0.5
likelihood_T_given_X = 0.99999
likelihood_T_given_not_X =0.2
# compute evidence, posterior
evidence_T = 0.5
posterior_X_given_T = 0.5
print(posterior_X_given_T)
0.5
This is a famous problem in statistics. The problem is simple, but the answer may not be what you expect! There are multiple ways to get the right answer, Bayesian inference is one way to get there.
(Image not available in this export)
You are a contestant on Monty Hall's game show, Let's Make a Deal, and you have made it to the final round! In this final round there are three doors -- behind one door is a brand new Ferrari sports car, and behind the other two there are just goats. At first, all three doors are closed, so each is equally likely to have either the car or one of the goats. Monty lets you pick a door (for ease of the problem setup, suppose you first pick door #1). Now, after you pick the door but before you open it, Monty says "but wait! are you sure you don't want to consider this door...?!" and then Monty opens another door you didn't pick (let's say, door #3) and there is a goat behind it. Monty then gives you an opportunity to change your guess. Now... a cute goat might be great prize since you love animals and are on a safari, but let's assume you are going for the Ferrari, so we will rule out picking door #3. Given this opportunity, should you switch your guess from door #1 to door #2, or stick with door #1, or does it not matter? Justify your answer by computing the posteriors $\mathcal{P}(d1 | D)$ and $\mathcal{P}(d2 | D)$, where $d\#$ indicates the outcome of the car being behind door $\#$, and $D$ represents our observation of the goat behind the door Monty just opened. Alternatively, you can also solve the problem using the odds ratio.
Your solution here:
This part of the tutorial is adapted from the version presented by David Wright at the IPTA 2025 student week. Various sections will ask you to write the code yourself -- it is encouraged to attempt each but if you get stuck or want to move on, you can copy from the code box which should include a working implementation.
Time for a problem in cosmology! Here we are going to use the method of standard candles. Suppose we use the luminosity of our standard candles to directly measure the radial comoving distance $\chi$ (in units of Mpc) as a function of the source redshift $z$, which is used to measure the universal expansion rate,
\begin{align*} \chi &= \int dz\frac{c}{H(z)}, \end{align*}where $H(z)$ is the Hubble parameter. A $\Lambda\rm{CDM}$ cosmology yields our model for the Hubble parameter,
\begin{align*} H(z) &= H_0\left[\Omega_{r,0}(1+z)^4 + \Omega_{m,0}(1+z)^3 + \Omega_{k,0}(1+z)^2 + \Omega_{\Lambda,0}\right]^{1/2}, \end{align*}where $H_0$ is the (present-day) Hubble constant in units of km/s/Mpc, $c$ is the speed of light in km/s, and $\Omega_{r,0}$, $\Omega_{m,0}$, $\Omega_{\Lambda,0}$, $\Omega_{k,0}$ are, respectively, the fractional energy densities for radiation, matter, dark energy, and curvature at present day. Note that these are normalized such that $\Omega_{r,0} + \Omega_{m,0} + \Omega_{\Lambda,0} + \Omega_{k,0} = 1$. For notational convenience later, we'll bundle all of these into a parameter vector $\vec{\Omega} = \{\Omega_{r,0}$, $\Omega_{m,0}$, $\Omega_{\Lambda,0}$, $\Omega_{k,0}\}$.
Dave compiled us some nice Hubble parameter data at last year's student workshop, so we will go ahead and use that directly and then we don't have to worry about solving any comoving distance integrals. We have many observations so going forward we will write the vector of hubble parameter measurements as $\vec{H}_{\rm obs}$ rather than just $D$. To clarify, $\vec{H}_{\rm obs}$ will be our measured data and $H(\vec{z})$ will be the predictions of our model.
df = pd.read_feather("data/h_z_measurements.feather")
df
Visualizing the data:
plt.figure()
plt.errorbar(df["z"], df["H_z"], yerr=df["H_z_err"], fmt="o", markersize=4)
plt.xlabel("z")
plt.ylabel(r"$H_{\rm obs}\ \left[km / s/Mpc\right]$")
plt.tight_layout()
plt.show()
Okay great, so we are basically fitting a model to something that looks roughly linear with some slight additional trend.
Next we are going to start modeling. As we go through, consider what assumptions we are making as we set up the model and how those assumptions may or may not impact the results.
First we should write a function that gives us $H(z)$ as a function of redshift and the model parameters. To start with the simplest version of this problem, it will be nice to assume all parameters take their $\Lambda\rm{CDM}$ values and infer just on $H_0$.
Here is the model again so you don't have to scroll up:
\begin{align*} H(z) &= H_0\left[\Omega_{r,0}(1+z)^4 + \Omega_{m,0}(1+z)^3 + \Omega_{k,0}(1+z)^2 + \Omega_{\Lambda,0}\right]^{1/2}, \end{align*}def H_z_model(z, H0, Omega_m=0.3, Omega_lambda=0.7, Omega_r=9e-5, Omega_k=0):
'''Return H(z) given the input redshift and parameters, assuming Lambda-CDM'''
# YOUR CODE HERE
Now for the Bayesian inference elements! The goal is to obtain estimates of the parameter $H_0$ with accurate uncertainties as given by our model and the observations. To do so in a Bayesian way, we should infer the posterior $\mathcal{P}(H_0 | \vec{H}_{\rm obs}, \vec{\Omega}, \Lambda\rm{CDM})$. Note we include all three elements $\vec{H}_{\rm obs}$ (our data), $\vec{\Omega}$ (the other model parameters) and $\Lambda\rm{CDM}$ (our model) in the conditional to remind ourselves about our assumptions, but it is very common to drop some of these terms from the conditional for notational convenience. We notice that since the posterior must be normalized by necessity, we do not need to evaluate the evidence $\mathcal{Z}(\vec{H}_{\rm obs} | \vec{\Omega}, \Lambda\rm{CDM})$, we can just draw samples directly following the proportionality
\begin{align} \mathcal{P}(H_0 | \vec{H}_{\rm obs}, \vec{\Omega}, \Lambda{\rm{CDM}}) \propto \mathcal{L}(\vec{H}_{\rm obs} | H_0, \vec{\Omega}, \Lambda{\rm{CDM}})\pi(H_0 | \vec{\Omega}, \Lambda{\rm{CDM}}), \end{align}and numerically normalize the samples later.
First let's define our prior on $H_0$. Priors must be chosen with care. Too narrow a prior can result in our pre-conceived notions about the problem impacting our inference results more than the data does. When first starting a problem it is good to use an uninformative prior, which is to say we want the likelihood itself to control our posterior distribution (likewise there are also times when informative priors are useful, e.g. if we have some other astrophysical measurement we'd like to account for). A uniform distribution is a oftentimes a good choice of uninformative prior, i.e.,
\begin{align} H_0 &\sim \mathcal{U}({\rm min}, {\rm max}) \\ \pi(H_0 | \vec{\Omega}, \Lambda\rm{CDM}) &= \begin{cases} \frac{1}{\rm{max} - \rm{min}} & {\rm{where}} \; {\rm{min}} < H_0 < {\rm{max}} \\ 0 & \rm{elsewhere} \end{cases}. \end{align}Implement this uniform prior below for some reasonable choice of min and max value of your choosing (referencing the above plot may be useful to decide what min and max values to choose). Note also for some reasons of numerical and analytic convenience it is convention to define the log of the prior, $\log\pi(H_0 | \vec{\Omega}, \Lambda\rm{CDM})$, in code.
def log_prior(H0):
"""Return the prior probability for a given value of H0"""
min_val = # set min bound of the prior
max_val = # set max bound of the prior
# YOUR CODE HERE
File "/tmp/ipykernel_3024/2931051460.py", line 3 min_val = # set min bound of the prior ^ SyntaxError: invalid syntax
Next, we notice there is some error around the $H(\vec{z})$ measurements. We will denote these errors via the vector $\vec{\sigma}$. As such, we need to update our data model from
\begin{align} \vec{H}_{\rm obs} = H(\vec{z}) \to \vec{H}_{\rm obs} = H(\vec{z}) + \vec{n}, \end{align}where $\vec{n}$ indicates some random noise vector deviating the observations away from the predicted model. We will model this noise as random draws from a multivariate Gaussian distribution given by the measurement errors $\vec{\sigma}$,
\begin{align} p(\vec{n} | \vec{\sigma}) = \prod_{i=1}^{N}\frac{1}{\sqrt{2\pi}\sigma_i}\exp\left(-\frac{n_i^2}{2\sigma_i^2}\right), \end{align}where $N$ is the number of datapoints. We can rewrite this in terms of our observables using $\vec{n} = \vec{H}_{\rm obs} - H(\vec{z})$ to convert into our likelihood function,
\begin{align} \mathcal{L}(\vec{H}_{\rm obs} | H_0, \vec{\Omega}, \vec{\sigma}, \Lambda{\rm CDM}) = \prod_{i=1}^{N}\frac{1}{\sqrt{2\pi}\sigma_i}\exp\left(-\frac{(H_{{\rm obs},i} - H(z_i))^2}{2\sigma_i^2}\right), \end{align}The log-likelihood function is then
\begin{align} \log\mathcal{L}(\vec{H}_{\rm obs} | H_0, \vec{\Omega}, \vec{\sigma}, \Lambda{\rm CDM}) = \sum_{i=1}^{N}\left[-\frac{(H_{{\rm obs},i} - H(z_i))^2}{2\sigma_i^2} - \frac{1}{2}\log(2\pi\sigma_i^2)\right]. \end{align}Now we are ready to implement.
def log_likelihood(df, H0):
"""Gaussian likelihood -- return the log probability to obtain the data given the model parameters"""
z = np.array(df['z'])
Hobs = np.array(df['H_z'])
sigma = np.array(df['H_z_err'])
# YOUR CODE HERE
Note this tangent can be skipped without loss of continuity in this tutorial:
The multivariate Gaussian is by far the most common choice of likelihood function for astronomical data, including for PTA analyses. This is partially for mathematical convenience but also this is motivated by the concept of the Central Limit Theorem as well as its relationship to least-squares regression.
This likelihood can be made even more compact notationally if we introduce the concept of a variance matrix $N_{ij} = \sigma_i^2\delta_{ij}$ which only includes diagonal elements due to the Kronecker delta $\delta_{ij}$ (indicating the measurement noise is uncorrelated between observations). This is convenient as it allows us to do away with the sum and express the likelihood as a matrix equation,
\begin{align} \log\mathcal{L}(\vec{H}_{\rm obs} | H_0, \vec{\Omega}, \vec{\sigma}, \Lambda{\rm CDM}) = -\frac{1}{2}(\vec{H}_{\rm obs} - H(\vec{z}))^T\mathbf{N}^{-1}(\vec{H}_{\rm obs} - H(\vec{z})) - \frac{1}{2}\log\det(2\pi\mathbf{N}), \end{align}where $\det$ indicates the determinant. While this extension goes beyond the scope of this tutorial, the matrix form of the Gaussian likelihood allows one to implement noise models where the noise is correlated between observations, which is done by adding elements to the off-diagonals of the matrix $\mathbf{N}$, at which point $\mathbf{N}$ has become a covariance matrix. Covariance matrix models are also known as Gaussian processes and these are usually how stochastic signals and noise processes are modeled in PTA gravitational wave searches.
We are now ready to infer $H_0$. Since we are only inferring one parameter, a very simple way to get our posterior and figure out the best value of $H_0$ and its uncertainty is to setup a simple grid search over different $H_0$ values. This can be implemented via the following algorithm:
N = # set number of gridpoints
H0_grid = # set up grid here
log_posterior = # initialize posterior array here
# SETUP A GRID SEARCH OVER H0 VALUES HERE
# note that exponentiating the unnormalized log posterior could lead to nonsense values numerically
# since this is not normalized anyways we will subtract off the max value
# to keep posterior values either finite or zero
posterior = np.exp(log_posterior - np.max(log_posterior))
plt.figure()
plt.plot(H0_grid, posterior, '.')
#plt.plot(H0_grid, log_posterior, '.') # option to plot log posterior
plt.ylabel(r'PDF (unnormalized)')
plt.xlabel(r'$H_0$')
plt.show()
The above plot should effectively give you an (unnormalized) probability density for the Hubble constant. It's worthwhile pausing here to see what you can figure out about your results
You can compute the best fit and credible error interval on $H_0$ most easily by generating samples from the posterior distribution. The fastest way to do this from a posterior array is using inverse transform sampling via the cumulative distribution function (CDF). A function to do this for you is written below
def inverse_transform_sample(H0_grid, posterior, N_samples=10000, make_plot=False):
CDF = np.zeros(len(H0_grid))
for i in range(len(H0_grid)):
CDF[i] = np.sum(posterior[:i+1])
# normalize the CDF by the highest value
CDF /= np.max(CDF)
# we can plot the CDF to see what it looks like if we want
if make_plot:
plt.figure()
plt.plot(H0_grid, CDF, '.')
plt.ylabel(r'PDF (unnormalized)')
plt.xlabel(r'$H_0$')
plt.show()
# generate samples from a uniform distribution
samples = np.random.uniform(size=N_samples)
# map samples to the CDF
sample_idxs = np.searchsorted(CDF, samples)
return H0_grid[sample_idxs]
H0_samples = inverse_transform_sample(H0_grid, posterior, N_samples=10000, make_plot=False)
# compute moments of the distribution -- mean and standard deviation
H0_mean = np.mean(H0_samples)
H0_sigma = np.std(H0_samples)
display(Math(fr'$H_0 = {H0_mean:0.3f} \pm {H0_sigma:0.3f} \; \rm{{km/s/Mpc}}'))
plt.figure()
plt.hist(H0_samples, bins=50, density=True, histtype='step', label='Full distribution')
plt.axvline(H0_mean, color='C1', ls='--', label='mean')
plt.axvspan(H0_mean-H0_sigma, H0_mean+H0_sigma, color='C1', alpha=0.2, label=r'$1\sigma$ range')
plt.ylabel(r'PDF')
plt.xlabel(r'$H_0$')
plt.legend()
plt.show()
Does your recovered value of $H_0$ and distribution look reasonable?
The last thing we should do is assess the model's fit to the data. An advantage of the Bayesian analysis is that we can use our final posterior samples distribution to map to a probability distribution over our model, $\mathcal{P}(H_0) \to \mathcal{P}(H(\vec{z}))$, and compare that to our original observations, $\vec{H}_{\rm obs}$. There are various types of quantitative checks that could be done here, but in this case a qualitative check will be enough.
plt.figure()
plt.errorbar(df["z"], df["H_z"], yerr=df["H_z_err"], fmt="o", markersize=4, label='Data')
plt.plot(df["z"], H_z_model(df["z"], H0_samples[0]), color='C1', lw=0.5, alpha=0.5, label='Model')
# plot a separate line for different draws from our posterior. 1000 samples will be enough
for H0 in H0_samples[1:1000]:
plt.plot(df["z"], H_z_model(df["z"], H0), color='C1', lw=0.5, alpha=0.1)
plt.xlabel("z")
plt.ylabel(r"$H(z)\ \left[km / s/Mpc\right]$")
plt.legend()
plt.tight_layout()
plt.show()
Ask yourself at this point:
If you are dissatisfied with the model's fit, think about what assumptions we made during the analysis, and how those could have impacted our results. We will attempt correct some of these assumptions in the next part.
Our last analysis was overly restrictive, because we were inferring just on $H_0$ itself, but not on any of the other cosmological parameters in $\vec{\Omega}$. In Bayesian inference parlance, we were using an overly restrictive prior, specifically delta function priors on the cosmological parameters:
\begin{align} p(\Omega_{m,0}) &= \delta(\Omega_{m,0} - 0.3), \\ p(\Omega_{\Lambda,0}) &= \delta(\Omega_{\Lambda,0} - 0.7), \\ p(\Omega_{r,0}) &= \delta(\Omega_{r,0} - 9 \times 10^{-5}), \\ p(\Omega_{k,0}) &= \delta(\Omega_{k,0} - 0). \end{align}However, in reality we might not know all of these cosmological parameters unless we infer them from the data. Let's relax our assumptions a little bit more and develop new priors based on that:
Based on this new information, the next most complicated thing to do is to sample $\Omega_{\Lambda,0}$ and treat $\Omega_{m,0}$ as a derived parameter via $\Omega_{m,0} = 1 - \Omega_{\Lambda,0} - \Omega_{r,0}$. So now we will sample a 2D parameter space including [$H_0$, $\Omega_{\Lambda,0}$]. Since we are still assuming $\Lambda{\rm CDM}$, we can reuse the likelihood. However, now we need a new prior and a new sampling technique.
First, let's write a new function for the log prior, $\log\pi(H_0, \Omega_{\Lambda,0} | \Lambda{\rm CDM})$. It is most common when performing multi-parameter inference to assume the parameters are uncorrelated with one another, i.e., $\pi(H_0, \Omega_{\Lambda,0}) = \pi(H_0)\pi(\Omega_{\Lambda,0})$. Based on the previous test, you should already have a good prior on $H_0$, so now we just need to figure out one for $\pi(\Omega_{\Lambda,0})$. Based on the above assumptions and your previous knowledge, consider what a good prior distribution would be for $\Omega_{\Lambda,0}$ and include that in your implementation. Note when implementing priors, there are often "good" and "bad" answers but rarely is there ever any one "right" answer.
def log_prior(H0, Omega_lambda):
"""Return the prior probability for a given value of H0 and Omega_lambda"""
# YOUR CODE HERE
Note our previous likelihood only accepted $H_0$ as an input, but not $\Omega_{\Lambda,0}$. Modify your previous likelihood here to allow both inputs.
Hint: It is not a bad idea to define a separate function $\Omega_{m,0} = 1 - \Omega_{\Lambda,0} - \Omega_{r,0} - \Omega_{k,0}$ to use throughout this section.
def log_likelihood(df, H0, Omega_lambda):
"""
Gaussian likelihood -- return the probability to obtain the data given the model parameters
This should be the same as your previous likelihood, just now accepting 2 parameter inputs
Don't forget that the value of Omega_matter depends on Omega_lambda!
"""
z = np.array(df['z'])
Hobs = np.array(df['H_z'])
sigma = np.array(df['H_z_err'])
# YOUR CODE HERE
Next we have the sampling method. We could use grid search again. However, this would be a bit slower now with complexity $\mathcal{O}(N_{\rm gridpoints}^2)$ instead of $\mathcal{O}(N_{\rm gridpoints})$ and it would be annoying (using the CDF is less straightforward, for example). If we increased the number of parameters, you can see the gridsearch method will fail very fast, because in principle we could have $n_d$ parameter in our model we have to infer at once, and then the complexity of a grid search scales as $\mathcal{O}(N_{\rm gridpoints}^{n_d})$. So now is a good opportunity to learn a more useful sampling method.
A classic method that scales well to up to dozens of parameters is Markov chain Monte Carlo (MCMC). The button below will give you an explanation about it, which you should read if you are unfamiliar with MCMC. However, in this tutorial we are primarily intereted in how to use MCMC to sample from a posterior distribution, rather than how to implement an MCMC algorithm ourselves. So, use the below explanation as reference (especially the jargon) but focus primarily on getting an intuition via the ensuing parts of the tutorial, and know that sooner or later you will run into problems with your choice of MCMC sampler, which will inevitably require you to setup your own MCMC algorithm :)
Here is the app from the lecture to see and play around with MCMC in action. For the conventional Metropolis-Hastings MCMC algorithm, set the algorithm to RandomWalkMH. Set the proposal distribution and watch as the algorithm samples form the posterior!
RandomWalkMH algorithm uses a simple Gaussian proposal distribution, i.e. it is proposing new samples as a Gaussian with some width $\sigma$ around the previous sample. Play around with the "Proposal $\sigma$" slider to change the width of the proposal distribution. How does it impact the efficiency of the sampler?emcee to sample the posterior¶Thankfully there are many python packages that exist that implement MCMC using a wide variety of jump proposals for us that work on many general problems, so we do not have to worry about it! Here we will use the common package emcee to do this, which is one of the easiest MCMC samplers to use out of the box. See the emcee tutorial here and associated docs to learn more.
Note emcee is an "ensemble sampler," meaning it natively samples using multiple chains in parallel with each other to exchange information. Moreover, running multiple chains and ensuring they each achieve the same rough distribution of samples under different starting points of the Markov chains is important to be more confident we have achieved the correct final distribution. Checking out the docs for emcee.EnsembleSampler, it is very simple to setup, we just need to specify three required inputs:
nwalkersndimlog_probability which is the log of our prior times our likelihood, given the input parameterslog_probability is setup for you using your previous log_prior and log_likelihood functions, just decide values nwalkers (your choice -- anywhere from 1 to 16 is a reasonable choice) and set ndim.
def log_probability(theta):
'''
Setting up the log prior + log likelihood for the emcee sampler.
Note the input theta is a 2D array containing a value of H0 and a value of Omega_lambda
'''
H0 = theta[0]
Omega_lambda = theta[1]
log_pi = log_prior(H0, Omega_lambda)
# don't bother evaluating the likelihood if the prior is negative
if not np.isfinite(log_pi):
return log_pi
return log_pi + log_likelihood(df, H0, Omega_lambda)
ndim = 2 # we are sampling in 2 parameters
nwalkers = 4 # feel free to change this
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability)
As for running the sampler, we finally just need to choose the initialization: What is a reasonable choice of initial parameters for $H_0$ and $\Omega_{\Lambda,0}$? You need a number of initial samples equal to nwalkers, i.e. the initial sample array should be a 2D array of shape (N_walkers, Nparameters). Some good choices might include:
Note the choice of initialization may impact how long it takes for the sampler to "burn in." See the emcee tutorial for some hints how to initialize the chain, or my version below if you are stuck (note my version is intentionally not the most efficient initialization)
initial_sample = # YOUR IMPLEMENTATION HERE
It is always a good idea to check that the log_probability function returns reasonable values for a few of your initial samples, so you rule out your code is not obviously broken before sampling.
for i in range(len(initial_sample)):
print(log_probability(initial_sample[i]))
Finally, we can run the sampler! We can choose how many samples to generate before we stop.
nsamples = 20000 # feel free to change this
# keep this line -- we should rebuild the sampler if we need to come back and restart this cell
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability)
# run the sampler
sampler.run_mcmc(initial_sample, nsamples, progress=True);
fig, axes = plt.subplots(ndim, figsize=(7, 4), sharex=True)
# resulting shape of the samples array is (Nsamples, Nwalkers, Nparams)
samples = sampler.get_chain()
labels = [r"$H_0$", r"$\Omega_{\Lambda,0}$"]
for i in range(ndim):
ax = axes[i]
ax.plot(samples[:, :, i], "k", alpha=0.3)
ax.set_ylabel(labels[i])
ax.yaxis.set_label_coords(-0.1, 0.5)
axes[-1].set_xlabel("step number")
axes[-1].set_xlim(0, len(samples))
fig.tight_layout()
fig.subplots_adjust(hspace=0)
plt.show()
This are called chain trace plots, and are useful to qualitatively determine if our samples are being drawn efficiently from the posterior. To ensure efficient sampling, we want to see that these chains look like "fuzzy caterpillars." It is likely that your trace plots do not look like fuzzy caterpillars at first. Some common outcomes:
You may see several of the chains disagree early on but quickly converge to a particular value. The phase before this convergence takes place is the "burn-in" phase, and we do not want to consider these samples as being part of the posterior distribution. You can set how many of these initial samples to throw away by supplying a value for the discard keyword argument to sampler.get_chain().
Another common issue for MCMC samplers I alluded to before is that not all chains may have independently converged to the same distribution, because one or more chains get "stuck" in some local maximum of the likelihood surface, which is not the true global maximum (this is why we need good jump proposals)! If you are seeing this, you may consider rerunning the sampler multiple times for different random initializations, changing the prior distribution, or changing the number of walkers, until results look better.
Let's make another plot to diagnose the 2nd issue. Specifically we can try overplotting posterior histograms of the different walkers to ensure the distributions look the same. We'll also look at the distribution of the log likelihood + log prior
discard = 1000 # set number of samples to burn here
fig, axes = plt.subplots(1, ndim+1, figsize=(7, 4))
samples = sampler.get_chain(discard=discard)
labels = [r"$H_0$", r"$\Omega_{\Lambda,0}$", r"$\log\pi + \log\mathcal{L}$"]
for i in range(ndim):
ax = axes[i]
for j in range(nwalkers):
ax.hist(samples[:, j, i], density=True, histtype='step', color=f"C{j}", bins=50)
ax.set_xlabel(labels[i])
ax.set_yticks([])
for j in range(nwalkers):
axes[-1].hist(sampler.lnprobability.T[discard:,j], density=True, histtype='step', color=f"C{j}", bins=50)
axes[-1].set_xlabel(labels[-1])
axes[-1].set_yticks([])
axes[0].set_ylabel('PDF')
fig.tight_layout()
fig.subplots_adjust(wspace=0)
plt.show()
Hopefully we see some nice posteriors that agree with each other! We should also obtain a distribution over $\Omega_{m,0}$, our derived parameter.
flat_samples = sampler.get_chain(discard=discard, thin=10, flat=True)
H0_samples = flat_samples[:,0]
Omega_lambda_samples = samples[:,1]
Omega_matter_samples = # GET OMEGA_M HERE
Another useful diagnostic is a corner plot. These are extremely visually appealing ways to see the covariances between our different model parameters.
labels = [r"$H_0$", r"$\Omega_{\Lambda,0}$", r"$\Omega_{m,0}$"]
fig = corner.corner(np.vstack([H0_samples, Omega_lambda_samples, Omega_matter_samples]).T,
labels=labels, levels=(0.68,0.95), show_titles=True);
Congratulations, you not only just measured the Universe's expansion rate, but also what percentage of the Universe is matter vs dark energy. Note the perfect covariance between $\Omega_{\Lambda,0}$ and $\Omega_{m,0}$ is by construction, but any remaining covariance with $H_0$ is an intrinsic feature of our model and/or the data.
Finally, we should also evaluate the fit to the data itself as we did before:
plt.figure()
plt.errorbar(df["z"], df["H_z"], yerr=df["H_z_err"], fmt="o", markersize=4, label='Data')
label = 'Model'
for H0, Omega_lambda, Omega_matter in zip(H0_samples[:1000],
Omega_lambda_samples[:1000],
Omega_matter_samples[:1000]):
Hz_sample = H_z_model(df["z"], H0, Omega_lambda=Omega_lambda, Omega_m=Omega_matter)
plt.plot(df["z"], Hz_sample, color='C1', lw=0.5, alpha=0.1, label=label)
label=None
plt.xlabel("z")
plt.ylabel(r"$H(z)\ \left[km / s/Mpc\right]$")
plt.legend()
plt.tight_layout()
plt.show()
Hopefully this is looking okay! This is the most surefire way to check the the model results are reasonable.
If everything went well sampling wise, we should have successfully fit our (simulated) cosmological data using a $\Lambda{\rm CDM}$ cosmology! That said, we may see our data is still not entirely consistent with the model. You may also be skeptical of the resulting model parameters we got out. We treated our $\Lambda{\rm CDM}$ model very well by allowing all three model parameters ($H_0$, $\Omega_{\Lambda,0}$, $\Omega_{m,0}$) to vary, but if there is still a discrepancy, perhaps $\Lambda{\rm CDM}$ is at fault?
Congratulations on making it this far! A final and important topic in Bayesian inference is to do hypothesis testing, i.e. to perform probabilistic inference on the models themselves, as opposed to just their parameters. This next test is inspired by the DESI DR2 results from last year.
Here we'd like to put $\Lambda{\rm CDM}$, which treats dark energy via a cosmological constant $\Lambda$, to the test against an alternative model, which assumes that dark energy may evolve over time. There are many physically motivated models for evolving dark energy, many of which can be simplified by changing the contribution of dark energy to the Hubble expansion from a constant $\Omega_{\Lambda,0} \to \Omega_{\Lambda,0}(1 + z)^{3(1 + w(a))}$, where we have introduced the smoothly-varying dark energy equation of state variable,
\begin{align} w(a) = w_0 + w_a(1 - a), \end{align}and $a = (1 + z)^{-1}$ is the cosmic scale factor, while $w_0$ and $w_a$ are two new parameters, where $w_0$ is the present day value of $w$ and $w_0 + w_a$ is the early-universe value of $w$. The form $\Omega_{\Lambda,0}(1 + z)^{3(1 + w(a))}$ tells us that $w_0 \ne -1$ implies that dark energy today is inconsistent with the notion of a cosmological constant, while $w_a \ne 0$ implies the existence of evolving dark energy.
We'd like to infer on the models themselves, and determine if and by how much one model is preferred by the data over the other. Labelling $\Lambda{\rm CDM}$ as $\mathcal{M}_\Lambda$ and the evolving dark energy model as $\mathcal{M}_w$, we'd like to determine the odds ratio
\begin{align} \frac{\mathcal{P}(\mathcal{M}_w | D)}{\mathcal{P}(\mathcal{M}_\Lambda | D)}. \end{align}We note this is proportional to the ratio of evidences given by each model, otherwise known as the Bayes factor,
\begin{align} \mathcal{B}^w_\Lambda \equiv \frac{\mathcal{Z}(D | \mathcal{M}_w)}{\mathcal{Z}(D | \mathcal{M}_\Lambda)}. \end{align}Recalling the definition of the model evidence as a marginal likelihood, solving this will amount to computing the integrals
\begin{align} \mathcal{B}^w_\Lambda = \frac{\int d\theta_\Lambda d\theta_w\mathcal{L}(D | \theta_\Lambda, \theta_w, \mathcal{M}_w)\pi(\theta_\Lambda, \theta_w | \mathcal{M}_w)}{\int d\theta_\Lambda\mathcal{L}(D | \theta_\Lambda \mathcal{M}_\Lambda)\pi(\theta_\Lambda | \mathcal{M}_\Lambda)}, \end{align}where $\theta_\Lambda$, $\theta_w$ label the model parameters of the $\Lambda{\rm CDM}$ and evolving dark energy models respectively.
There are several ways we can approach solving this equation, all of which are frequently used in the PTA community. We will briefly summarize and then pick one:
Since it is the easiest to implement, we are going to evaluate the preference for evolving dark energy just using method 3., the Savage-Dickey density ratio. Afterwards, you may attempt to implement the other methods yourself as an exercise (comparing these methods could also make for a nice IPTA mini-project)!
The next steps will not be so different from what we did before. Essentially, we need to implement a new model, likelihood, and prior for evolving dark energy, and then attempt the sampling problem once again. Here is the full model you'll have to implement.
\begin{align} H(z) = H_0 \left(\Omega_m (1+z)^3 + \Omega_r (1+z)^4 + \Omega_{\Lambda} (1+z)^{3 (1+ w_0 +w_a)} e^{-3 w_a \left(\frac{z}{1+z}\right)}\right)^{1/2}. \end{align}The next set of cells are open-ended for you to implement yourself (no solutions manual!) but the structure will be nearly the same as before (so you you should be copying a lot of code rather than writing a lot of code from scratch). Note that since this is a tougher problem to sample in 4D, sp you may want to consider using more samples and more walkers than you did before.
def H_z_evolving_model(z, H0, w0, wa, Omega_m, Omega_lambda, Omega_r=9e-5, Omega_k=0):
# YOUR CODE HERE
def log_prior(H0, Omega_lambda, w0, wa):
"""Return the prior probability for a given value of H0 and Omega_lambda"""
# YOUR CODE HERE
def log_likelihood(df, H0, Omega_lambda, w0, wa):
"""
Gaussian likelihood -- return the probability to obtain the data given the model parameters
This should be the same as your previous likelihood, now accepting 4 parameter inputs
Don't forget that the value of Omega_matter depends on Omega_lambda!
"""
# YOUR CODE HERE
def log_probability(theta):
'''
Setting up the log prior + log likelihood for the emcee sampler.
Note the input theta is a 4D array containing a values of H0, Omega_lambda, w0, wa
'''
# YOUR CODE HERE
# SETUP SAMPLER HERE
# SETUP INITIAL SAMPLE HERE
# RUN SAMPLER HERE
Note if you run into the error:
ValueError: Initial state has a large condition number. Make sure that your walkers are linearly independent for the best performance
you can set skip_initial_state_check=True to bypass.
# MAKE TRACE PLOTS
# MAKE CORNER PLOT
# COMPARE YOUR EVOLVING DARK ENERGY MODEL WITH YOUR DATA
Once you are confident that your solution is correct, comment on you new parameter distributions and your model's fit to the data.
Assuming you have a nice set of samples w0_samples and wa_samples, we can also replicate Figure 14 from the DESI paper below:
from matplotlib import lines
labels = [r"$w_0$", r"$w_a$"]
fig = corner.corner(np.vstack([w0_samples, wa_samples]).T,
labels=labels, levels=(0.68,0.95), show_titles=True,
color='C2', smooth=0.5, plot_density=False, fill_contours=True);
axes = np.array(fig.axes).reshape((2,2))
# resize the axes
axes[0,0].set_xlim([-1.2, 0])
axes[1,0].set_xlim([-1.2, 0])
axes[1,0].set_ylim([-1.6, 0.2])
axes[1,1].set_xlim([-1.6, 0.2])
corner.overplot_lines(fig, [-1,0], color="k")
handles = [lines.Line2D([],[],color='C2',label='Results')]
handles += [lines.Line2D([],[],color='k',label=r'Expectations of $\Lambda{\rm CDM}$')]
fig.legend(handles=handles)
plt.show()
Using the concept of the posterior probability and Savage-Dickey density ratio, can you place a rough bound on the odds ratio for evolving dark energy as given by this data?
Congratulations! You have reached the end of the tutorial. Check out the following resources, or consider the following challenges:
Go back and change the underlying data we were using by inflating the measurement uncertainties $\sigma$ and rerun your analyses to see how noisier data translates to higher uncertainties on our model parameters.
Go back and see the impact of changing the prior distribution. An informative prior would be, for example, a Gaussian centered on $H_0 = 68$ km/s/Mpc (roughly based on the CMB results).
Calculate the Bayesian evidence for both models in the last example by implementing a nested sampling algorithm (see below).
There are a lot more resources on Bayesian inference and specific methods out here.
To see many of these methods and more discussed in the context of pulsar timing, I recommend the GWB methods paper Johnson et al. 2024&q=%20%20firstauthor%3A%22Johnson%22%20%20year%3A(2024)%20pulsar&sort=date%20desc%2C%20bibcode%20desc&p=0).
Here is a great 5 part blog series on Bayesian statistics by Jake VanderPlas here. In particular it touches on comparisons with Frequentist statistics that were not really covered here. The last part on model selection is also very useful.
Some further methods you may explore include:
Parallel tempering can be used in tandem with MCMC as a computationally expensive but very powerful method to create jump proposals. It relies on sampling the real likelihood surface in parallel with "smoothed out" versions of the likelihood surface, and proposing jumps using the smoothed version. This is helpful to explore difficult posterior morphologies such as multimodal distributions.
PTMCMCSampler is a common software used in the PTA community (the PT part actually stands for "parallel tempering" not pulsar timing). It is built on an "adaptive Metropolis" algorithm, in which the MH jump proposals are tuned adaptively during sampling, but was designed to use MPI for a truly parallel imepentation of the parallel tempering algorithm.
Johnson et al. 2024&q=%20%20firstauthor%3A%22Johnson%22%20%20year%3A(2024)%20pulsar&sort=date%20desc%2C%20bibcode%20desc&p=0) includes further details on parallel tempering specifically.
If you are not already frustrated with Metropolis-Hasting MCMC, you will be. At this point, you may consider Hamiltonian Monte Carlo (HMC) as a sampling method, which works more efficiently than MH MCMC by using gradients of the likelihood with respect to the parameters at the current location for jump proposals. This method has recently gained more popularity due to widespread use of programming backends such as JAX that enable automatic gradient computations.
Honestly I have not found many great tutorials on HMC. If anyone knows of one please let me know!
NumPyro for an example python library with an implementation of HMC (as well as many other MCMC implementations)etudes package (by Gabe Freedman)I briefly mentioned hierarchical Bayesian inference in the context of Bayesian model selection. In hierarchical Bayesian inference, the idea is that your various model settings, such as the choice of model you are using, or the prior bounds on parameters you are inferring, may themselves be uncertain, and you include that uncertainty explicitly in the analysis with another layer of prior distributions.
For example, PTA data analysis typically employs a hierarchical likelihood, in which the
Firstly, if you'd like to implement the nested sampling solution to the previous problem, I recommend checking out the nautilus-sampler. This is installed on your VM so you can play around with it.
If you prefer JAX and GPUs there is also BlackJAX-NS
A review of nested sampling in physics
Probabilistic inference (including Bayesian methods) are so widespread that people have developed specialized libraries that change the way you approach writing the code to perform inference using PPLs -- Probabilistic Programming Languages (or libraries). These are especially useful for bookkeeping if you are performing hierarchical inference across many layers of a model and also want to keep track of derived parameters.
Some Useful Python PPLs:
I ran out of time to write more about the above topics, but they are also interesting. Come and ask us things!