Introduction to Bayesian Statistics#
The key idea behind Bayesian statistics requires a complete change of perspective regarding the way we approach problems. In a (probably far-fetched) analogy with physics, it would be similar to how we perceive gravity in a Newtonian versus a General Relativistic universe. In this notebook we will discuss:
Bayes’ theorem
The main Bayesian components: prior, likelihood, posterior, and evidence
Using Bayesian statistics for parameter estimation
Marginalization over nuissance parameters (and the difference with conditioning)
Setup#
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
import requests
from IPython.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
url = "https://raw.githubusercontent.com/AstroStat-Academy/assets-public/main/colab/clone_and_cd_colab.py"
colab = requests.get(url).text
exec(colab)
url = "https://raw.githubusercontent.com/AstroStat-Academy/assets-public/main/styles/plot_style.py"
style = requests.get(url).text
exec(style)
The problem: Counting photons…#
An X-ray imaging observation hands us a list of events: for each photon the telescope detects, we record its position and energy. Since we are literally counting photons, we call them counts, and the count rate is just the number of counts per second. The count rate scales with the flux of the source in the telescope’s energy band.
Figure 1. Left - Some galaxy seen in optical wavelengths. Right - The same galaxy observed in X-rays. Each dot represents a photon that has been detected. The energy of the photon is known as well (not shown). The ellipse defines the extent of the source based on the image in the optical. |
Source-detection algorithms group the events into source regions — one list of photons per source. From there, a chain of modeling steps turns those photons into the physical quantity we actually care about:
The integrated energy of a source region estimates the energy collected by the telescope — once the detector response and absorption are accounted for.
With enough counts, we can fit the source’s spectrum.
Modeling the detector, the intergalactic absorption, and the spectrum gives the bolometric flux — the flux a perfect experiment would measure.
Finally, given the source’s distance, flux converts to luminosity. This an intrinsic property of the source, which is usually what we’re after.
How many counts are we going to get?#
Discuss with your teammate, then report.
[Spoiler] (click here to expand)
We can get an estimate of the photon counts by simply dividing.
COUNT_TO_LUMINOSITY = 1.0e38 # in erg/s
source_luminosity = 5.0e38 # in erg/s
expected_counts = source_luminosity / COUNT_TO_LUMINOSITY
print(f"Expected counts: {expected_counts:.3f}")
Discuss with your teammate, then report.
[Spoiler] (click here to expand)
This is a counting problem. We expect each measurement of a photon to be independent of the others, so… it’s Poisson!
Reminder: The Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time, if these events occur with a known constant mean rate and independently of the time since the last event.
Let’s plot a Poisson distribution.
outcomes = np.arange(0, int(expected_counts*3)+1, 1)
pmf = st.poisson(expected_counts).pmf(outcomes)
plt.figure()
plt.plot(outcomes, pmf, "ks")
plt.xlabel("Counts per measurement")
plt.ylabel("Probability mass function")
plt.show()
Discuss with your teammate, then report.
[Spoiler] (click here to expand)
Since it’s Poisson, the standard deviation is the square root of the expected value!
Let’s compute the uncertainty on the measured counts. We can either use a formula or a scipy function.
counts_uncertainty = expected_counts**0.5
print(f"Uncertainty: {counts_uncertainty:.3f}")
Estimating the luminosity of a source#
Discuss with your teammate, then report.
[Spoiler] (click here to expand)
The Poisson distribution above peaks both at 5 and 4 counts. If we got 4 counts, would \(5\times 10^{38}\,\rm erg\,s^{-1}\) still be a good answer for the expected luminosity of the source? There is something weird going on! We have one measurement and we ask a question about the source. Also, we are always going to get integer multiples of the count/luminosty factor this way. But nature doesn’t care about how many photons we measured, how far away we are from the source, etc.
Exercise 1: Select three offsets (can be negative, positive, or zero) around the mean value and inspect the resulting distributions. What do you notice?
outcomes = np.arange(0, int(expected_counts*3)+1, 1)
plt.figure()
# try different offsets around the expected mean value and plot the PMF
for offset, marker in zip([...,...,...], ["s", "o", "d"]):
mean = expected_counts + offset
plt.plot(outcomes, st.poisson(mean).pmf(outcomes), ":", marker=marker, alpha=0.8, label=f"Mean = {mean:.4g}")
plt.axvline(expected_counts, color="0.5", alpha=0.5, label="k=5")
plt.ylim(ymin=0)
plt.legend(loc="upper right")
plt.xlabel("Counts")
plt.ylabel("Probability mass function")
plt.show()
[Solution] (click here to expand)
for offset, marker in zip([-0.5, 0.0, 0.3], ["s", "o", "d"]):
We have multiple (actually infinite) hypotheses for the source luminosity#
Before, we had the source luminosity fixed, and we wondered what we would measure if we repeated the experiment many times. This is the frequentist approach. We also assigned some uncertainty to the data.
A frequentist assigns probabilities to data - parameters are fixed. The probability is a frequency of data outcomes.
But when asking about the nature of things (as is the case in all applied sciences), we are not interested in assigning uncertainty to the data. The data are what they are, and in some cases they may be fundamentally unique, e.g. because some experiments cannot be repeated. We need to assign probabilities to the quantity of interest. This is the Bayesian approach:
A Bayesian assigns probabilities to hypotheses. The probability is a degree of belief in a value of a parameter.
Figure 2. A schematic representation of (a) deductive logic or pure mathematics, and (b) inductive logic or plausible reasoning. |
Bayes Theorem: The Foundation of Bayesian Inference#
To navigate from what we can model, \(P(\textrm{data} | \textrm{hypothesis})\), to what we want to know, \(P(\textrm{hypothesis} | \textrm{data})\), we use Bayes Theorem that arises from the basic rules of probability.
The joint probability of two events (or propositions) A and B occurring, \(P(A, B)\), can be written in two ways:
Since \(P(A,B)\) is the same in both expressions, we can equate them:
Rearranging this gives us Bayes Theorem:
In the context of data and hypotheses, we write this as:
Here, \(I\) represents all our background information and assumptions.
We made everything conditional to \(I\), the background information: all relevant knowledge we have about the problem we are solving but are not part of the data or the specific hypothesis being tested.
For example, if we are testing whether a die is fair, the data are a sequence of outcomes (e.g., 1, 6, 3, 5, 3, 1), the hypothesis could be the probability of rolling a 6 is \(p_6=1/6\), and \(I\) = \(\big\{\) all dice have 6 sides and are cubes, a fair die has equiprobable sides, the laws of physics as we know them, \(\cdots \big\}\).
For brevity, we usually omit writing the \(I\) in the equations, but its always there… somewhere in the background!
All the terms of Bayes Theorem have specific names:
Prior Probability: \(P(\text{hypothesis} | I)\) is the degree of belief in the hypothesis before looking at the current data. It encapsulates our previous knowledge, physical constraints, or initial assumptions. For example it can just be a range.
Likelihood: \(P(\text{data | hypothesis}, I)\) is the probability of observing the data given that the hypothesis is true. This is what we often model based on our understanding of the data generation process.
Posterior Probability: \(P(\text{hypothesis} | \text{data}, I)\) is the (updated) degree of belief in the hypothesis after considering the data.
Evidence (or Marginal Likelihood): \(P(\text{data} | I)\) is the overall probability of observing the data, averaged over all possible hypotheses.
In Bayesian Analysis we assign degrees of belief to hypotheses (prior), which we “update” using experimental data (via the likelihood), to arrive at a new, refined degree of belief in the hypotheses (posterior).
A Deeper Look at the Components:#
Priors \(P(\text{hypothesis} | I)\)#
The prior represents our state of knowledge about the hypothesis before we consider the current dataset. It can come from previous experiments, theoretical predictions, or fundamental principles. The choice of prior can be influential, especially with limited data, and is a key aspect of Bayesian modeling.
Priors can be broadly categorized as informative or uninformative or, more accurately, weakly informative.
Informative Priors: These reflect specific, existing information about the parameter before observing the current data. This information might come from previous experiments, physical laws, or established theories. For example, if measuring the mass of a planet, an informative prior might restrict the mass to be positive and perhaps within a range plausible for planets.
Uninformative/Weakly Informative Priors: These are intended to let the data speak for themselves as much as possible, or to represent a state of relative ignorance. However, a truly “uninformative” prior is often elusive. The most common approach is to use uniform priors. Those assign equal probability to all possible values of a parameter within a given range. For a parameter \(\theta\) bounded by \(\theta_{min}\) and \(\theta_{max}\), \(P(\theta|I) \propto 1\).
fig, axs = plt.subplots(1, 1, figsize=(6, 4))
# Uniform Prior
x_uniform = np.linspace(-2, 12, 500)
pdf_uniform = st.uniform(loc=0, scale=10).pdf(x_uniform) # Uniform between 0 and 10
axs.plot(x_uniform, pdf_uniform, label='Uniform Prior P(θ|I) ~ U(0, 10)')
axs.set_title('Uniform Prior')
axs.set_xlabel('Parameter value θ')
axs.set_ylabel('Probability Density')
axs.legend()
plt.show()
The choice of prior should always be justified and its impact on the results assessed, especially if strong claims are made from weakly informative data.
Likelihood \(P(\text{data} | \text{hypothesis}, I)\)#
The likelihood is what connects the data to the parameters: it tells us how probable the observed data are for a given value of the parameters. The choice of likelihood states our assumptions about how the data were generated, what process produced them, and what kind of noise they carry. Its form follows from the nature of the data: continuous measurements, discrete counts, binary outcomes, and so on. Each of these processes corresponds to a different distribution.
Two cases cover most of what we meet in practice:
Gaussian (Normal) likelihood — for continuous data scattered symmetrically around a mean, e.g. a brightness measurement with normally-distributed noise. For \(N\) points \(x_i\) drawn from a Gaussian with mean \(\mu\) and standard deviation \(\sigma\):
Poisson likelihood — for count data: independent events at a constant mean rate, e.g. the number of X-ray photons \(k\) detected when the expected count is \(\lambda\) (this is our running example).
When none of the standard distributions is appropriate, the likelihood is built from a model of the measurement itself. This can be derived from the instrument’s physics, combined from several noise sources, or, when no closed form exists, evaluated numerically by simulating the data-generating process.
Two things are worth keeping in mind:
Although the likelihood is a probability distribution over the data, we read it as a function of the parameters at the fixed, observed data. In a true Bayesian way, this informs us about the hypothesis.
It is a modeling choice: if its assumptions are badly violated, the inference can mislead, which is why heavier-tailed alternatives (e.g. a Student’s-\(t\) in place of a Gaussian) are sometimes used to blunt the effect of outliers.
# Example Plots of Common Likelihood Functions
fig, axs = plt.subplots(1, 2, figsize=(12, 5))
axs = axs.ravel() # Flatten the 1x2 array of axes for easy iteration
# 1. Gaussian Likelihood
mu_g, sigma_g = 0, 1.0 # Example parameters: mean, std_dev
x_g = np.linspace(mu_g - 4*sigma_g, mu_g + 4*sigma_g, 200)
pdf_g = st.norm.pdf(x_g, loc=mu_g, scale=sigma_g)
axs[0].plot(x_g, pdf_g, label=f'Gaussian(μ={mu_g}, σ={sigma_g})')
axs[0].set_title('Gaussian Likelihood')
axs[0].set_xlabel('Data Value (x)')
axs[0].set_ylabel('Probability Density P(x|μ,σ)')
axs[0].legend()
axs[0].grid(True, alpha=0.5)
# 2. Poisson Likelihood
lambda_p = 5 # Example parameter: expected rate λ
k_p = np.arange(0, lambda_p * 3 + 1)
pmf_p = st.poisson.pmf(k_p, mu=lambda_p)
axs[1].plot(k_p, pmf_p, 'o-', label=f'Poisson(λ={lambda_p})')
axs[1].set_title('Poisson Likelihood')
axs[1].set_xlabel('Number of Events (k)')
axs[1].set_ylabel('Probability Mass P(k|λ)')
axs[1].legend()
axs[1].grid(True, alpha=0.5)
plt.tight_layout()
plt.show()
The Evidence \(P(\text{data} | I)\):#
The evidence is the probability of the data, averaged over all possible values of the hypothesis parameters, weighted by their prior probabilities. If \(\theta\) represents the parameters of our hypothesis, then:
In parameter estimation problems (where we are trying to find the best values of \(\theta\)), the evidence \(P(\text{data} | I)\) is a constant for a given dataset and choice of model. It ensures that the posterior probability distribution is properly normalized (i.e., integrates/sums to 1). While it often cancels out or is ignored in parameter estimation as we often care about \(P(\text{hypothesis} | \text{data}) \propto P(\text{data} | \text{hypothesis}) P(\text{hypothesis})\), the evidence is crucial for model comparison, a topic we will explore later. It allows us to compare how well different models (different forms of hypotheses or priors) can explain the data.
Back to our counting experiment: what does Bayes say about it…#
where in this case, since we care only for the value of the parameters, we can ignore the evidence term (which is just a normalization constant).
Assigning a prior#
We need to choose a prior \(P(\lambda | I)\) that describes our belief on the source luminosity and, more specifically, its expected counts, \(\lambda\).
Option 1: Uniform Prior over a finite range (Principle of Indifference)
We may have a constraint, e.g. for galaxies that contain an Active Galactic Nucleus \(\lambda\) is above 20, but we know from other measurements that this is not the case for our galaxy (this is background information \(I\)):
The constant \(C\) would be \(1/20\) in order to normalize the prior. However, for parameter inference we are actually interested in the shape of the posterior \(P(\lambda | k) \propto P(k | \lambda) P(\lambda)\), thus we can often ignore normalization constants.
Option 2: Uniform (Improper) Prior for \(\lambda \ge 0\) (Agnostic approach)
Maybe we have no other information and want to be agnostic about \(\lambda\) (other than it must be non-negative):
This is an improper prior because its integral over \([0, \infty)\) diverges. However, it can still lead to a proper posterior. Since it is just a constant, it does not affect the posterior shape.
Assigning a likelihood#
We can use the Poisson likelihood, just like in the frequentist approach:
What is the Posterior?#
In the case of a proper prior we get:
while for an improper one the posterior is:
Discuss with your teammate, then report.
[Spoiler] (click here to expand)
No! Bayesian inference is not just frequentism + priors. Even if the formula for the posterior \(P(\lambda|k)\) looks proportional to the likelihood \(P(k|\lambda)\) when using a flat prior, the interpretation is fundamentally different.
The posterior \(P(\lambda|k)\) is a probability distribution for the parameter \(\lambda\), given the observed data \(k\). The likelihood \(P(k|\lambda)\) is a probability distribution for the data \(k\), given a fixed parameter \(\lambda\).
The posterior \(P(\lambda | k) \propto \lambda^k e^{-\lambda}\) (for \(k\) fixed, \(\lambda \ge 0\)) is a probability density function of the continuous quantity \(\lambda\). In fact, it is a Gamma distribution: specifically, \(\text{Gamma}(\text{shape}=k+1, \text{scale}=1)\) or \(\text{Gamma}(\alpha=k+1, \beta=1)\) depending on parameterization (e.g., if PDF is \(f(x;\alpha,\beta) = \frac{\beta^\alpha x^{\alpha-1}e^{-\beta x}}{\Gamma(\alpha)}\)).
As a consequence, Maximum Likelihood Estimation (see previous lecture) is rooted in the frequentist idea of how often data like ours would be observed if a parameter had a certain value. Bayesian inference provides a probabilistic statement about the parameter itself, given the observed data and prior beliefs.
Exercise 2: Compute the posterior values for the trial \(\lambda\) values. What is the best-fitting value for \(\lambda\)?
measure_counts = 5.0
lambda_values = np.linspace(0, int(measure_counts*3)+1, 100)
posteriors = st.gamma(...).pdf(...)
plt.figure()
plt.plot(lambda_values, posteriors, "k-")
plt.xlabel(r"Expected counts ($\lambda$)")
plt.ylabel(r"Posterior density $P(\lambda \mid k)$")
plt.axvline(5.0, color="r", ls=":", label="Measured counts")
plt.legend(loc="upper right")
plt.show()
[Solution] (click here to expand)
Using a uniform prior means that the posterior is \(P(\lambda | k) \propto \lambda^k e^{-\lambda}\), which for k counts is exactly the Gamma(k+1, scale=1) probability density function over the rate \(\lambda\).
posteriors = st.gamma(a=measure_counts + 1, scale=1).pdf(lambda_values)
Parameter Marginalization: Focusing on What Matters#
Often, our statistical models involve multiple parameters. However, we might only be interested in a subset of these parameters, while the others are necessary for a complete model but not of direct interest. These latter parameters are often called nuisance parameters.
For example, when measuring the flux, \(\mu\), of a star, the noise level or background uncertainty, \(\sigma\), in our measurements might also be an unknown parameter that we need to account for in our model. While we need to estimate \(\sigma\) to properly assess \(\mu\), our primary goal is the value of \(\mu\), making \(\sigma\) a nuisance parameter in this context.
Bayesian inference provides a natural and powerful way to deal with nuisance parameters through marginalization. If we have a joint posterior probability distribution for all parameters, say \(P(\theta_1, \theta_2, ..., \theta_k | \text{data}, I)\), and we are only interested in a specific parameter (or subset of parameters) \(\theta_1\), we can obtain its marginal posterior probability distribution by integrating (or summing, for discrete parameters) the joint posterior over all other nuisance parameters (\(\theta_2, ..., \theta_k\)):
For a two-parameter case (e.g., \(\mu\) and \(\sigma\)):
The result, \(P(\theta_1 | \text{data}, I)\), is the posterior probability distribution for \(\theta_1\) alone, fully accounting for the uncertainty in the nuisance parameters. This means that our inferences about \(\theta_1\) (like its credible interval) have properly incorporated the fact that the other parameters were not known precisely.
This is a significant advantage of the Bayesian framework. It allows us to propagate uncertainty from all parts of our model into the final inference for the quantities we care about. In frequentist statistics, dealing with nuisance parameters can sometimes be more complex, often involving techniques like profiling the likelihood or using specific test statistics designed to be insensitive to them.
Figure 3. Difference between marginalization and conditioning. The starting point is the joint (two-dimensional in this case) probability distribution shown in the middle panel. Marginalization (top) performs the integral over a set (of usually nuissance) parameters. Conditioning (bottom) assumes a fixed value of a set of parameters, essentially ``cutting'' the posterior surface with a plane in the N-dimensional parameter space - in this two-dimensional case this plane corresponds to the horizontal dashed black line in the middle panel. The vertical dashed organe and blue lines are there simply to guide the eye. |
Example: Estimating Stellar Flux with Background as a Nuisance Parameter#
Let’s apply marginalization in a practical setting.
Suppose we observe a noisy measurement of the total brightness from a patch of sky. This total brightness includes contributions from:
The flux \(F\) of a star (our parameter of interest),
The background \(B\) from surrounding light or instrumental effects (a nuisance parameter), and
Gaussian noise with known standard deviation \(\sigma\).
We model the observed value \(x\) as:
which means that we should detect a value fluctuating around a mean (flux + background) due to the random noise.
We will place a broad prior on \(F\), reflecting that the star could be relatively bright or dim, and a more informative prior on \(B\), reflecting prior knowledge about typical background levels:
Although both \(F\) and \(B\) are unknown, we are only interested in the posterior for \(F\). To obtain it, we marginalize over the nuisance parameter \(B\):
Exercise 3: Define the priors, likelihood, joint posterior, and marginal posterior for \(F\).
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.integrate import quad
# Simulated data
true_flux = 10.0
true_background = 2.0
sigma = 1.0
# One observation (flux + background + noise)
np.random.seed(42)
obs = np.random.normal(loc=true_flux + true_background, scale=sigma)
print(f"Observed value: {obs:.2f}")
# Define priors
def prior_flux(F):
return ...
def prior_background(B):
return ...
# Likelihood
def likelihood(obs, F, B):
return ...
# Posterior (unnormalized, joint in F and B)
def joint_posterior(F, B):
return ...
# Marginal posterior over F (integrate out B)
def marginal_posterior_F(F):
integrand = lambda B: joint_posterior(F, B)
result, _ = quad(..., ..., ...)
return result
# Evaluate marginal posterior over a grid
F_vals = np.linspace(5, 15, 200)
posterior_vals = np.array([marginal_posterior_F(F) for F in F_vals])
# Normalize
posterior_vals /= np.trapz(posterior_vals, F_vals)
# Plot
plt.figure(figsize=(8, 4))
plt.plot(F_vals, posterior_vals, label='Marginal posterior of flux')
plt.axvline(true_flux, color='k', linestyle='--', label='True flux')
plt.xlabel('Flux')
plt.ylabel('Posterior density')
plt.legend()
plt.title('Marginalization over nuisance parameter (background)')
plt.show()
[Solution] (click here to expand)
# Define priors
def prior_flux(F):
return norm.pdf(F, loc=10, scale=5) # broad prior
def prior_background(B):
return norm.pdf(B, loc=2, scale=1) # informed prior
# Likelihood
def likelihood(obs, F, B):
return norm.pdf(obs, loc=F + B, scale=sigma)
# Posterior (unnormalized, joint in F and B)
def joint_posterior(F, B):
return likelihood(obs, F, B) * prior_flux(F) * prior_background(B)
# Marginal posterior over F (integrate out B)
def marginal_posterior_F(F):
integrand = lambda B: joint_posterior(F, B)
result, _ = quad(integrand, -np.inf, np.inf)
return result
# EOF