Assignment: Is the Spectral Line Really There?
Bayesian Inference
A story
A faint astronomical source is observed by a photon-counting detector. A theoretical model — say,
emission from a particular atomic transition in a distant cloud — predicts that the spectrum of
this source should contain a narrow emission line sitting on top of a smooth continuum background.
More precisely, the model predicts that the flux as a function of frequency should have the form
(x − x0 )2
fnet (x) = B + A exp − ,
|{z} 2w2
| {z }
background spectral line
where x denotes frequency, B is a constant continuum level, and the bump is a Gaussian centered
at x0 with width w and peak amplitude A. Independent astrophysical arguments fix the location
and width of the line at
x0 = 1.5, w = 0.1,
but the amplitude A and the continuum B are not predicted from first principles — they depend
on the physical state of the source and are exactly what we want to learn from the observation.
Here is the twist. The detector is a photon counter, not a flux meter. In each of M = 101 narrow
frequency bins x1 , x2 , . . . , xM (uniformly spaced on [1, 2]) it simply records an integer: the number of
photons that happened to arrive in that bin during the exposure. Photon arrivals are fundamentally
a discrete, random process, and for a source of constant mean flux the number of counts in a bin is
Poisson-distributed about its expected value. Concretely, the number of counts in bin k satisfies
nk ∼ Poisson(Dk ), Dk = n0 fnet (xk ), k = 1, 2, . . . , M,
independently across bins, where n0 = 5 is an instrumental constant (absorbing the exposure time,
effective collecting area, and so on; assume it is known from calibration). What the observer actually
sees is therefore not a smooth curve but a jagged sequence of integer counts — the underlying
fnet (x) is corrupted by irreducible Poisson (“shot”) noise, and must be recovered statistically.
This raises a genuinely interesting question. When you stare at the noisy counts, it may look like
there is a bump in the middle; it may also look like the bump is consistent with an ordinary upward
fluctuation of the background. How confidently, given only the data, can we say that the line is
there at all, and that its amplitude is in the range predicted by the theory? That is the question you
will answer in this assignment, using Bayesian inference and the Metropolis algorithm you have just
learnt in class.
1
The data
The observed photon counts n = (n1 , n2 , . . . , nM ) with M = 101, recorded in the frequency windows
xk uniformly spaced on [1, 2], are
4, 8, 14, 18, 7, 8, 8, 10, 8, 6, 12, 13, 9, 6, 7, 9, 8, 6, 9, 12, 12, 7, 9, 7,
9, 17, 7, 12, 12, 13, 13, 15, 10, 16, 18, 20, 24, 18, 14, 28, 23, 15, 20, 24, 26, 23, 20, 19,
28, 26, 29, 29, 19, 31, 17, 15, 20, 24, 12, 26, 18, 23, 20, 15, 20, 12, 12, 10, 6, 9, 14, 13,
14, 12, 8, 7, 9, 13, 13, 11, 8, 14, 8, 8, 8, 12, 9, 8, 10, 10, 5, 10, 7, 4, 13, 12,
10, 9, 9, 10, 10,
listed in the order k = 1, 2, . . . , M .
Setting up the inference
Before any sampling, you need to turn the story above into a posterior distribution on (A, B).
Because the counts in different bins are independent, the joint probability of the entire dataset
factorizes over k. Each factor is a Poisson probability mass function whose mean depends on (A, B)
through Dk (A, B) = n0 B + A exp(−(xk − x0 )2 /2w2 ) . Taking the log turns the product into a
sum, and any terms that do not depend on (A, B) are irrelevant for inference — they only shift the
log-posterior by a constant, and so drop out of the Metropolis acceptance ratio.
For the prior, we will take a flat (uniform) prior on the rectangle
(
1 if Amin ≤ A ≤ Amax and Bmin ≤ B ≤ Bmax ,
p(A, B) ∝
0 otherwise,
with
Amin = Bmin = 0.1, Amax = Bmax = 6.
This encodes the weakest kind of prior information: we are willing to entertain any line amplitude
and any background level in the range [0.1, 6], with no preference among them. (If the theory had
made a sharper prediction about A, we could have encoded it here, but we are deliberately letting
the data do the talking.)
The physical question
Suppose the theoretical model, when combined with independent estimates of the source’s distance
and luminosity, predicts a line amplitude in the range
2.5 ≤ A ≤ 3.5.
A rival model predicts a much weaker line, essentially indistinguishable from the background. Given
the noisy data above, how strongly should we believe the first model? The relevant quantity is the
posterior probability Z 3.5
Pr 2.5 ≤ A ≤ 3.5 n = p(A | n) dA,
2.5
where p(A | n) is the marginal posterior of A obtained by integrating p(A, B | n) over all values of B.
Computing this integral directly is unpleasant, but MCMC makes it almost trivial: if {(Ai , Bi )}N s
i=1
2
are Ns samples from the joint posterior, then
Ns
1 X
Pr 2.5 ≤ A ≤ 3.5 n ≈ 1[2.5, 3.5] (Ai ),
Ns
i=1
where 1[2.5,3.5] (·) is the indicator of the interval [2.5, 3.5]. Notice that the Bi never enter this sum;
simply ignoring them is what performs the marginalization over B.
What you are asked to do
1. Look at the data. Enter the 101 counts into your program and plot nk against xk . Describe in
one or two sentences what you see. Is there visibly a line, or is the bump debatable? (Keep your
honest first impression — it will be interesting to revisit it after the analysis.)
2. Derive the log-posterior. Starting from the Poisson model nk ∼ Poisson(Dk (A, B)) and the
independence of the bins, write down the likelihood p(n | A, B) as a product over k, and then
derive the log-likelihood log p(n | A, B) as a function of (A, B), dropping any additive constants
that do not depend on (A, B). Combine this with the flat prior to write down the log of the
(unnormalized) posterior log p(A, B | n), being careful about what happens on the boundary of
the prior rectangle. Explain briefly why it is enough, for the Metropolis algorithm, to know the
log-posterior only up to an additive constant.
3. Implement the target. Turn the expression from Question 2 into a function that takes a pair
(A, B) and returns log p(A, B | n), returning −∞ whenever (A, B) leaves the prior rectangle.
4. Run the Metropolis sampler. Using a Gaussian random-walk proposal with standard deviation
σ = 0.2 applied independently to each coordinate, generate a chain of length N = 100,000 starting
from (A(0) , B (0) ) = (1, 1.5), and discard the first 10,000 iterations as burn-in. Report the empirical
acceptance rate and comment on whether it is in a reasonable range.
5. Visualize the posterior.
(a) Produce a scatter plot of the post-burn-in samples in the (A, B) plane. Comment on the
shape of the joint posterior: is it roughly elliptical? Do you see any correlation between A
and B, and does the direction of that correlation make physical sense? (Hint: think about
which combinations of A and B produce roughly the same total number of counts.)
(b) Produce separate histograms of the marginal samples of A and of B. For each, report the
sample mean and sample standard deviation.
6. Answer the physical question. Using the post-burn-in samples, compute the Monte Carlo
estimate of
Pr 2.5 ≤ A ≤ 3.5 n .
State your numerical answer. In one short paragraph, interpret it: does the data favor the first
theoretical model, the rival model, or does it leave the question open? Compare this quantitative
conclusion with the first impression you recorded in Question 1.
7. Sensitivity checks. (Short answers suffice.)
3
(a) Rerun the sampler with proposal standard deviations σ = 0.01 and σ = 2.0, keeping
everything else the same. Report the acceptance rate in each case and describe qualitatively
what goes wrong (or right) with the exploration of the posterior.
(b) Rerun the sampler with a different random seed (or a different initial state inside the
prior rectangle) and check that your estimate of Pr(2.5 ≤ A ≤ 3.5 | n) does not change
appreciably. This is a minimal check of Monte Carlo reproducibility.
What to hand in
A single document (notebook export or PDF) containing your derivation for Question 2, your code,
the requested plots, the numerical answers, and brief written commentary wherever the question
asks for interpretation. Keep the writing concise; a paragraph per question is usually enough.