BST 570-1: Probabilistic Machine Learning, Fall 2026
Submission cutoff: Oct 6th 11:59 PM. Create a new file named a1.qmd, submit it on a branch of your course repository, and open a pull request into main. Make sure your code and its output appear in the rendered document. Check that quarto render a1.qmd succeeds; failure to render will result in a 10% grade deduction.
Ground rules on AI are in the syllabus: discuss freely, but do not submit generated code verbatim, and write the explanations yourself. You will be asked to explain your solution in class.
RNA sequencing produces, for each sample \(n\) and gene \(j\), a count \(y_{n,j}\): the number of reads that mapped to gene \(j\) in sample \(n\). A sample is one cell in single-cell data or one tissue sample in bulk data. The simplest model treats the reads as arriving independently at a gene-specific rate,
\[ y_{n,j} \mid \lambda_j \sim \text{Poisson}(\lambda_j), \quad n = 1, \ldots, N, \qquad\text{independently given } \lambda_j, \]
and places a Gamma prior on the rate,
\[ \lambda_j \sim \text{Gamma}(a, b), \qquad p(\lambda) = \frac{b^a}{\Gamma(a)}\, \lambda^{a-1} e^{-b\lambda}, \quad \lambda > 0, \]
with \(b\) a rate, so the prior mean is \(\frac{a}{b}\). Note that scipy.stats.gamma takes a scale, which is \(1/b\).
Genes are modelled independently, so until part (e) we work with a single gene and suppress the subscript \(j\). Most parts of this problem ask for a derivation or a written answer and for code; submit both, with the code’s output visible in the rendered document. Parts (c) and (f) are writing only. The data are simulated so that everyone’s numbers agree:
import numpy as np
import scipy.stats as ss
from scipy.special import gammaln, logsumexp
import matplotlib.pyplot as plt
rng = np.random.default_rng(570)
N, lam_true = 20, 42.3
y = rng.poisson(lam_true, size=N)
print(f"N = {N}, sum = {y.sum()}, mean = {y.mean():.3f}")
print(y)N = 20, sum = 854, mean = 42.700
[47 40 31 42 38 47 52 43 60 38 35 45 43 39 49 53 42 40 30 40]The appendix of the week 2 slides did this derivation on one slide. Do it yourself here, without looking, and then compare.
Write down the likelihood of all \(N\) counts as a function of \(\lambda\), multiply by the prior, and drop every factor that does not involve \(\lambda\). Recognize the kernel and state the posterior with its parameters.
Then do it a second time using the exponential-family recipe: write the Poisson as \(h(y)\exp(\eta\, T(y) - A(\eta))\), identify \(\eta\), \(T\) and \(A\), and read the posterior update off \(\tau \to \tau + \sum_n T(y_n)\), \(\nu \to \nu + N\). Confirm that the two routes agree.
Throughout this problem, use a weak prior with mean \(50\) reads that is worth one tenth of a sample:
a0, b0 = 5.0, 0.1# TODO: posterior parameters a1, b1 from a0, b0 and the data
# TODO: plot the prior and posterior densities on one axis, with a vertical line at lam_trueWrite: the denominator of Bayes’ theorem never appeared in your derivation. What is it called, and why were you able to skip it?
Show that the posterior mean can be written as
\[ \mathbb{E}[\lambda \mid y_{1:N}] = w \,\frac{a}{b} \;+\; (1 - w)\, \bar y, \qquad \bar y = \frac{1}{N}\sum_n y_n, \]
and find \(w\) in terms of \(b\) and \(N\).
# TODO: compute w and check the formula against a1 / b1 numerically
# TODO: repeat with the prior Gamma(500, 10), which has the same mean 50 but a different weightWrite:
Week 1 separated inference from prediction. The posterior answers the inference question. Prediction averages the likelihood of a new observation over it:
\[ p(y_{N+1} \mid y_{1:N}) = \int_0^\infty p(y_{N+1} \mid \lambda)\, p(\lambda \mid y_{1:N})\, d\lambda . \]
Compute this integral in closed form and simplify. The result is a standard two-parameter discrete distribution. Which one? Check your pmf against the documentation of scipy.stats if you are unsure.
Write: without using the name of the distribution, find the mean and variance of \(y_{N+1}\) from the laws of total expectation and total variance, conditioning on \(\lambda\). Show that the mean is the posterior mean \(\mu = a_1 / b_1\) and the variance is \(\mu + \mu^2 / a_1\). Which of the two terms is the posterior variance of \(\lambda\), and what happens to it as \(N \to \infty\)? Confirm both against the known formulas for the distribution you identified.
A second run. Suppose the first \(8\) samples were sequenced in one run and the remaining \(12\) in a second. Update the prior with the first batch, then use the resulting posterior as the prior for the second batch.
# TODO: (a_A, b_A) after the first 8 samples; (a_B, b_B) after adding the remaining 12
# TODO: compare (a_B, b_B) with (a1, b1) from part (a)Write:
Unequal depth. Real samples are not sequenced to the same depth. If sample \(n\) has a known size factor \(s_n\) (its depth relative to a typical sample), the model becomes
\[ y_n \mid \lambda \sim \text{Poisson}(s_n \lambda). \]
Redo the derivation in (a) and show that the posterior is \(\text{Gamma}\big(a + \sum_n y_n,\; b + \sum_n s_n\big)\).
rng_s = np.random.default_rng(7)
s = np.exp(rng_s.normal(0, 0.4, size=N))
s = s / np.exp(np.log(s).mean()) # normalize so the geometric mean is 1
y_s = rng_s.poisson(s * lam_true)
print(np.round(s, 2))
print(y_s)[1.14 1.28 1.02 0.79 0.95 0.76 1.16 1.94 0.93 0.89 1.38 1.31 1.18 0.78
1.12 1.5 0.66 0.95 0.53 0.68]
[42 56 42 36 40 36 46 92 35 38 62 43 50 31 33 59 34 42 25 28]# TODO: posterior for lambda from y_s using the size factors, and the posterior you would
# get by ignoring them and treating y_s as Poisson(lambda). Compare the two means.Write: in the exponential-family update, which of the two quantities, \(\tau\) or \(\nu\), changed and which did not? The “sample size” in the update is no longer \(N\). What is it, and what does that say about what a sample contributes to the posterior? In which direction does ignoring the size factors bias the estimate here, and why?
Now the samples come from two groups, \(z_n \in \{0, 1\}\), for example control and treated, or two cell types. The label \(z_n\) is known. For each gene,
\[ y_{n,j} \mid z_n, \lambda_j^0, \lambda_j^1 \sim \text{Poisson}\big(\lambda_j^{z_n}\big), \]
and we want to test, gene by gene, the hypothesis \(\lambda_j^0 = \lambda_j^1\). Two hundred genes, ten samples per group. The first forty genes are truly differentially expressed, twenty up by a factor of two and twenty down by a factor of two; the rest are not. Use the truth only to evaluate at the end.
rng_e = np.random.default_rng(572)
J, n0, n1 = 200, 10, 10
z = np.array([0] * n0 + [1] * n1)
lam0 = 10 ** rng_e.uniform(0.7, 2.3, size=J) # baseline rates from 5 to 200 reads
fold = np.ones(J); fold[:20] = 2.0; fold[20:40] = 0.5
lam1 = lam0 * fold
Y = rng_e.poisson(np.where(z[:, None] == 0, lam0[None, :], lam1[None, :])) # (N, J)
truly_de = fold != 1.0
print(Y.shape, " first sample, first 8 genes:", Y[0, :8])(20, 200) first sample, first 8 genes: [ 40 15 142 168 66 10 7 139]Two ways to build the test, both from what you already derived. Use the same prior as in part (a), \(\text{Gamma}(a_0, b_0) = \text{Gamma}(5, 0.1)\), for every rate.
(i) Through the posterior. Under the two-rate model, the posteriors of \(\lambda_j^0\) and \(\lambda_j^1\) are independent Gammas. Say why in one sentence. Draw \(K = 20{,}000\) pairs from them and estimate \(P(\lambda_j^1 > \lambda_j^0 \mid y)\) for every gene at once (an array of shape (K, J) does it). Report the two-sided version \(p_j = 2\min\big(P,\, 1 - P\big)\).
(ii) Through the marginal likelihood. The marginal likelihood \(p(y_{1:N})\) is the denominator we never had to compute in part (a). Here it does some work. Show that
\[ p(y_{1:N}) = \frac{\Gamma(a + s)}{\Gamma(a)} \, \frac{b^a}{(b + N)^{a + s}} \, \prod_{n=1}^{N} \frac{1}{y_n!}, \qquad s = \sum_n y_n , \]
and implement it in log space with gammaln. Check it on the single gene from part (a) by Monte Carlo: average the likelihood over \(K = 100{,}000\) draws from the prior, in log space with logsumexp, as in Lab 2, Task 5.
Then compare two models for each gene. \(M_0\): one rate for all \(20\) samples. \(M_1\): one rate per group, each with its own \(\text{Gamma}(a_0, b_0)\) prior, so that the marginal likelihood is a product of two factors of the form above. Report the log Bayes factor
\[ \log \text{BF}_j = \log p(y_{\cdot, j} \mid M_1) - \log p(y_{\cdot, j} \mid M_0) \]
for every gene, and, for gene \(j = 100\) only, also the maximized log likelihood of each model with \(\hat\lambda = \bar y\) within each block.
(iii) Prior sensitivity. Rerun (i) and (ii) with the prior \(\text{Gamma}(50, 1)\): the same mean of \(50\) reads, but now worth one full sample instead of a tenth. Report how the two-sided probabilities and the log Bayes factors change, and how many of the top \(40\) genes by each statistic are still truly differentially expressed.
# TODO: (i) p_two (J,) from posterior draws
# TODO: (ii) log_ml(y, a, b) with gammaln; Monte Carlo check on the part (a) gene; log_bf (J,)
# TODO: for each statistic, the 40 top-ranked genes and how many of them are truly DE
# TODO: how many of the 160 null genes have log_bf > 0, and how many have p_two < 0.05
# TODO: scatter of log_bf against log2 of the ratio of group sample means, coloured by truly_de
# TODO: (iii) repeat (i) and (ii) with a0, b0 = 50, 1 and recompute the counts aboveWrite:
In practice, RNA-seq counts from biological replicates are not modelled as Poisson. The standard tools for differential expression, of which the two most widely used are edgeR and DESeq2, model each count directly as
\[ y_n \sim \text{NegBin}\big(\text{mean } \mu,\ \text{variance } \mu + \phi\mu^2\big), \]
with a dispersion \(\phi\) per gene that is estimated from the data. Typical estimated dispersions are between \(0.01\) and \(0.5\). Where does this model come from? Each sample \(n\) is given its own rate \(\lambda_n\), and the rates are drawn from a common Gamma distribution that describes how expression varies from one mouse, or one patient, to the next:
\[ \lambda_n \sim \text{Gamma}(\alpha, \beta), \qquad y_n \mid \lambda_n \sim \text{Poisson}(\lambda_n), \qquad n = 1, \ldots, N, \text{ independently}. \]
This is a hierarchical model: the parameters \((\alpha, \beta)\) are shared by the group, the rates \(\lambda_n\) are not.
Derive the marginal distribution of one count, \(p(y_n \mid \alpha, \beta) = \int_0^\infty p(y_n \mid \lambda_n)\, p(\lambda_n \mid \alpha, \beta)\, d\lambda_n\). You have done this integral already in part (c); only the parameters differ. Then re-express the result in terms of its mean \(\mu\) and its variance \(\mu + \phi\mu^2\): find \(\mu\) and \(\phi\) as functions of \((\alpha, \beta)\). Finally, write down the likelihood of \((\alpha, \beta)\) given all \(N\) counts, and note that \(\lambda_1, \ldots, \lambda_N\) appear nowhere in it.
Write: in part (c) you integrated a Poisson against \(\text{Gamma}(a_1, b_1)\), the posterior of the single rate \(\lambda\) shared by every sample, and got a negative binomial with dispersion \(1/a_1\), about \(0.001\) on our data. Here you integrated a Poisson against \(\text{Gamma}(\alpha, \beta)\) over the sample-specific rate \(\lambda_n\) and got a negative binomial with dispersion \(\phi\). The two integrals have the same form. Say what is random in each case: what does \(\text{Gamma}(a_1, b_1)\) describe about \(\lambda\), and what does \(\text{Gamma}(\alpha, \beta)\) describe about \(\lambda_n\)? Then use your expression for \(\phi\) to explain why \(1/a_1\) shrinks to zero as \(N\) grows while \(\phi\) does not. Which of the two is the right model for \(20\) technical replicates of one library, and which for \(20\) patients?
A graduate admissions committee reads an application. It sees a recommendation letter and an SAT score. It does not see the grade the letter is based on, how hard the course was, or how capable the student is, and the last of these is what it wants to know. This is the student network of Koller and Friedman (2009, Probabilistic Graphical Models, Fig. 3.4).
Shaded nodes are what the committee observes. The grade has three values; everything else is binary. No code is expected in this problem. Everything is worked by hand from the graph and the tables below, and your submission for it is written answers and arithmetic.
| \(P(D = \text{easy})\) | \(P(D = \text{hard})\) |
|---|---|
| 0.6 | 0.4 |
| \(P(I = \text{low})\) | \(P(I = \text{high})\) |
|---|---|
| 0.7 | 0.3 |
| \(I\) | \(D\) | \(P(G = \text{A} \mid I, D)\) | \(P(G = \text{B} \mid I, D)\) | \(P(G = \text{C} \mid I, D)\) |
|---|---|---|---|---|
| low | easy | 0.30 | 0.40 | 0.30 |
| low | hard | 0.05 | 0.25 | 0.70 |
| high | easy | 0.90 | 0.08 | 0.02 |
| high | hard | 0.50 | 0.30 | 0.20 |
| \(I\) | \(P(S = \text{high} \mid I)\) |
|---|---|
| low | 0.05 |
| high | 0.80 |
| \(G\) | \(P(L = \text{strong} \mid G)\) |
|---|---|
| A | 0.90 |
| B | 0.60 |
| C | 0.01 |
For each claim below, decide it using the Bayes ball rules from lecture. No code. Give the verdict, and for each path between the two variables name the node where the ball is blocked, or the node and rule through which it gets past.
| claim | |
|---|---|
| 1 | \(D \perp I\) |
| 2 | \(D \perp I \mid G\) |
| 3 | \(D \perp I \mid L\) |
| 4 | \(S \perp L\) |
| 5 | \(S \perp L \mid I\) |
| 6 | \(D \perp S\) |
| 7 | \(D \perp S \mid L\) |
Write:
The committee never observes a student’s ability \(I\), and never sees the grade \(G\) behind a letter. But from hundreds of past applicants it has learned which courses are hard, so treat \(D\) as known.
Today’s applicant took a course known to be hard, scored high on the SAT, and has a weak letter: \(D = \text{hard}\), \(S = \text{high}\), \(L = \text{weak}\). Work by hand from the tables.
If you knew the grade. Show that \[ P(I \mid D, G, S, L) \;\propto\; P(I)\, P(G \mid I, D)\, P(S \mid I), \] say which variable dropped out and which independence from the graph lets it, and compute \(P(I = \text{high} \mid D = \text{hard}, G = g, S = \text{high})\) for \(g = \text{A}, \text{B}, \text{C}\).
Write: three variables survive on the right hand side. Describe each one’s relation to \(I\) in the graph. One of them, \(D\), shares no edge with \(I\). Why must it be there?
What the committee actually knows. It knows neither \(I\) nor \(G\). Compute
\[ P(I = \text{high} \mid D = \text{hard}, S = \text{high}, L = \text{weak}), \,\, \text{ and } \]
\[ P(G \mid D = \text{hard}, S = \text{high}, L = \text{weak}). \]
Write: the same weak letter is read differently depending on the course. Which structure in the graph produces this, and what is it called?