Assignment 2: A hidden Markov model for daily step counts

BST 570-1: Probabilistic Machine Learning, Fall 2026

Submission cutoff: Oct 26th Noon. Create a new file named a2.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 a2.qmd succeeds; failure to render will result in a 10% grade deduction.

How this assignment is organized. Each part has some background and then a numbered list of deliverables, which is everything you have to hand in for that part. Number your answers the same way. Include the code behind every number, table and plot; derivations and explanations are a few lines of your own algebra or prose.

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.

The data

A mobile device (or a wearable device) records the number of steps its owner takes each day. Over a year the counts drift between stretches of a few thousand steps a day and stretches of ten thousand or more, and the switches between stretches are not announced. A study has collected these daily counts from \(N = 40\) subjects. Most have a full year; a few have only four to eight weeks. On some days the device was not worn, and the count for that day is missing.

The data file data/steps.csv holds a \(40 \times 365\) array, one subject per row, no header. Missing days are empty fields and load as NaN. Column \(t\) is day \(t + 1\) of the year. Not everyone wore the device all year: subjects joined the study on different days, and some stopped wearing the device before the year ended. The file data/steps_subjects.csv gives, for each subject, the day of the year of their first and last measurement, first_day and last_day (1-based, inclusive). We call the slice of row \(i\) between those two days the record of subject \(i\): their sequence of daily counts, of length n_days, with NaN on the days inside it when the device was not worn. The entries outside the slice are empty and are not part of the record. For example, a subject with first_day 5 and last_day 50 has a \(46\)-day record, which in code is Y[i, 4:50]: the columns are 0-based, the days are 1-based. The first day of every record is observed.

import numpy as np
import pandas as pd
import scipy.stats as ss
from scipy.special import logsumexp, gammaln
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt

Y = pd.read_csv("data/steps.csv", header=None).to_numpy()
subjects = pd.read_csv("data/steps_subjects.csv")
n_days = subjects["n_days"].to_numpy()
start = subjects["first_day"].to_numpy() - 1          # 0-based column of the first measurement
stop = subjects["last_day"].to_numpy()                # slice end, exclusive
N, T = Y.shape
K = 3

records = [Y[i, start[i]:stop[i]] for i in range(N)]  # NaN inside a record = missing day
in_record = (np.arange(T)[None, :] >= start[:, None]) & (np.arange(T)[None, :] < stop[:, None])
print(f"N = {N}, T = {T};  full-year records: {(n_days == T).sum()}")
print(f"short records (days): {sorted(n_days[n_days < T])}")
print(f"missing days within records: {np.isnan(Y[in_record]).mean():.1%}")
print(subjects[n_days < T].head())
N = 40, T = 365;  full-year records: 28
short records (days): [33, 38, 41, 43, 44, 45, 51, 51, 53, 55, 55, 55]
missing days within records: 7.8%
    subject  first_day  last_day  n_days
0         1        184       221      38
2         3         93       135      43
6         7        274       326      53
9        10        207       261      55
14       15        310       353      44

Problem 1 ignores the order of the days and models each subject’s counts as a mixture: every day is a draw from one of \(K = 3\) activity levels. Problem 2 puts the days back in order and lets today’s level depend on yesterday’s (Markov chain). The emission model built in Problem 1 carries over unchanged.

Problem 1: A mixture model for daily step counts (30 points)

For one subject, each observed day \(t\) has a hidden label \(z_t \in \{1, 2, 3\}\) drawn independently from \(\text{Categorical}(\pi)\), and a count whose distribution depends on the label through a mean \(\lambda_k\). Below, \(\theta_k\) stands for whatever parameters component \(k\) has: \(\lambda_k\) for the Poisson, and \((\mu_k, \phi)\) for the negative binomial of part (b). Days are independent, so missing days are simply left out: work with the observed counts of each record. This is the mixture EM of Lab 3 with the Gaussian replaced by a count distribution. The responsibilities are \[ r_t(k) = p(z_t = k \mid y_t) = \frac{\pi_k\, p(y_t \mid \theta_k)}{\sum_{k'} \pi_{k'}\, p(y_t \mid \theta_{k'})}, \] the weight update is their average, and the rate update is the responsibility-weighted mean \[ \bar y_k = \frac{\sum_t r_t(k)\, y_t}{\sum_t r_t(k)} . \]

Fit one model per subject, for the \(28\) subjects with a full-year record only. Every fitted quantity is therefore subject-specific; when that matters we write a superscript, as in \(\bar y^{(i)}_k\) for subject \(i\)’s weighted mean of component \(k\), and “median over subjects” means the median of the \(28\) values \(\bar y^{(1)}_k, \dots, \bar y^{(28)}_k\). Initialize the rates at the \(20\)th, \(50\)th and \(80\)th percentiles of the subject’s observed counts and the weights at \(1/3\), and order the components by rate when you report them.

(a) Poisson mixture (10 points)

The count is Poisson: \(y_t \mid z_t = k \sim \text{Poisson}(\lambda_k)\).

To check the model after fitting, compute for each subject \(i\) and each component \(k\) the responsibility-weighted variance \[ s^{2\,(i)}_k = \frac{\sum_t r_t(k)\, (y_t - \bar y^{(i)}_k)^2}{\sum_t r_t(k)} , \] with \(r_t(k)\) and \(y_t\) those of subject \(i\).

Deliverables.

  1. (1 point) With \(\pi = (0.3, 0.2, 0.5)\) and \(\lambda = (2000, 6000, 11000)\), compute the log likelihood \(\sum_t \log \sum_k \pi_k\, p(y_t \mid \lambda_k)\) over the observed days. It should be \(-8206.60\) for subject 1, which is row 0 of Y, and \(-49533.36\) for subject 2.
  2. (3 points) Fit the Poisson mixture to each full-year subject and tabulate the fitted rates and weights for the first five of them.
  3. (3 points) Compute the ratio \(s^{2\,(i)}_k / \bar y^{(i)}_k\) for every subject and component, and make one plot of it against the component’s mean \(\bar y^{(i)}_k\), one point per subject and component, with the Poisson’s prediction, a ratio of \(1\), as a horizontal line.
  4. (3 points) Say what the ratio would be if the Poisson model were right, and what the plot says about how the variance grows with the mean. In Assignment 1 part (f), overdispersion arose because each unit had its own rate; say what the unit is here.

(b) Negative binomial mixture (20 points)

Replace the Poisson with the negative binomial in its mean and dispersion form, \[ p(y \mid \mu, \phi) = \frac{\Gamma(y + 1/\phi)}{\Gamma(1/\phi)\, y!} \left(\frac{1}{1 + \phi\mu}\right)^{1/\phi} \left(\frac{\phi\mu}{1 + \phi\mu}\right)^{y}, \qquad \mathbb{E}[y] = \mu, \quad \text{Var}[y] = \mu + \phi\mu^2 . \] This is scipy.stats.nbinom with n = 1/phi and p = 1/(1 + phi*mu). Each component has its own mean \(\mu_k\); the dispersion \(\phi\) is shared by the three components.

Background: maximizing a function of one parameter numerically. The derivative with respect to \(\phi\) contains \(\psi(y_t + 1/\phi) - \psi(1/\phi)\) for every day, where \(\psi\) is the digamma function, and no rearrangement isolates \(\phi\). This is the normal situation for most models, and the remedy is routine: hand the function to an optimizer. For one parameter, scipy.optimize.minimize_scalar is the tool. Three conventions matter.

  1. The optimizer minimizes. Pass the negative log likelihood.
  2. Optimize an unconstrained quantity. The dispersion must be positive. Rather than telling the optimizer about the constraint, optimize \(\eta = \log\phi\) and set \(\phi = e^{\eta}\) inside the function. This also puts the search on the right scale. A dispersion acts multiplicatively, so doubling \(\phi\) matters the same whether it goes from \(0.01\) to \(0.02\) or from \(0.1\) to \(0.2\); in \(\eta\) both are the same step, and the optimizer’s tolerance, which is absolute, then resolves every magnitude equally well.
  3. Say where to look. Give a bounded interval in \(\eta\) that is wide enough to contain any plausible answer and use method="bounded". The optimizer evaluates the function about a dozen times inside it.

Here is the whole recipe on one simulated sequence of counts with a known mean. The function being minimized is the ordinary negative binomial negative log likelihood; the responsibilities come later.

rng = np.random.default_rng(570)
mu_sim, phi_sim = 6000.0, 0.08
y_sim = rng.negative_binomial(1 / phi_sim, 1 / (1 + phi_sim * mu_sim), size=365)

def nb_logpmf(y, mu, phi):
    r = 1.0 / phi
    return ss.nbinom.logpmf(y, r, r / (r + mu))

def neg_loglik(eta):                       # eta = log phi
    return -nb_logpmf(y_sim, mu_sim, np.exp(eta)).sum()

res = minimize_scalar(neg_loglik, bounds=(np.log(1e-4), np.log(5.0)), method="bounded")
phi_hat = np.exp(res.x)
print(f"phi_hat = {phi_hat:.4f}   (truth {phi_sim}),  function evaluations: {res.nfev}")

phis = np.logspace(-3, 0.5, 200)
plt.figure(figsize=(6, 3.2))
plt.plot(phis, [-neg_loglik(np.log(p)) for p in phis])
plt.axvline(phi_hat, color="C3", ls="--", label=f"maximizer {phi_hat:.3f}")
plt.xscale("log"); plt.xlabel(r"$\phi$"); plt.ylabel("log likelihood"); plt.legend(); plt.tight_layout(); plt.show()
phi_hat = 0.0750   (truth 0.08),  function evaluations: 11

Always draw the curve of the log likelihood against the parameter, as above, the first time you optimize a new function. It tells you whether the function is smooth and has one maximum, whether your bounds contain it, and whether the optimizer found it. Here the profile is unimodal, which is typical for the negative binomial dispersion, and the maximizer is a little off the truth because a single year of counts is not much information about a variance.

In the M-step the function to maximize is the emission part of \(Q\), with \(\mu_k\) already updated: every observed day enters once per component, weighted by \(r_t(k)\) and evaluated at that component’s mean. Write it as a function of \(\eta\) and call minimize_scalar exactly as above. The monotonicity argument from lecture only needs the M-step to increase \(Q\), not to maximize it exactly, so EM with a numerical M-step is still EM, sometimes called generalized EM. The digamma terms are inside scipy.stats.nbinom.logpmf, so you never need them explicitly.

Deliverables.

  1. (6 points) Write down the \(Q\) function of EM for this mixture, differentiate it with respect to \(\mu_k\), and show that for any fixed \(\phi\) the maximizer in \(\mu_k\) is the weighted mean \(\bar y_k\). Then write down the one-dimensional function of \(\phi\) that the M-step has to maximize numerically.
  2. (1 point) Compute the negative binomial mixture log likelihood with the weights and rates of (a) and \(\phi = 0.08\). It should be \(-324.66\) for subject 1 and \(-3271.61\) for subject 2.
  3. (6 points) Fit the negative binomial mixture to each full-year subject. Tabulate \(\mu\), \(\pi\) and \(\phi\) for the first five subjects, give the median and range of \(\hat\phi\) across subjects, and recompute the ratio \(s^{2\,(i)}_k / \bar y^{(i)}_k\) of (a) with the negative binomial responsibilities. Make one plot of it against the model’s value \(1 + \hat\phi^{(i)}\hat\mu^{(i)}_k\), one point per subject and component, with the identity line.
  4. (3 points) For each full-year subject, label each observed day by its most probable component, \(\arg\max_k r_t(k)\). Count the number of times, summed over subjects, that this label changes from one observed day to the next. Plot one subject’s counts over \(60\) days with each day coloured by its label and the three fitted means as horizontal lines. Keep the fitted mixtures and these labels; Problem 2 adds the HMM’s labels to this picture.
  5. (4 points) Describe two things that could go wrong if you optimized over \(\phi\) directly with bounds \((0, 5)\) instead of over \(\log\phi\).

Problem 2: A hidden Markov model (70 points)

The components overlap: a day with \(7000\) steps could belong to the medium or the high component, and the mixture of Problem 1 has no way to decide, because it treats days as independent. What could decide is the neighbours: a \(7000\)-step day between two \(11000\)-step days is a quiet day in a high stretch. The hidden Markov model uses them.

The hidden label of day \(t\) is now an activity regime that lasts for weeks, such as an injury, a vacation, an ordinary routine or a training block, and today’s regime depends on yesterday’s.

  • Observations: \(y_t\) is the number of steps on day \(t\), for \(t = 1, \ldots, T_i\), with some \(y_t\) missing.
  • Hidden states: \(z_t \in \{1, 2, 3\}\) is the regime on day \(t\): low, medium and high.
  • Initial distribution: \(z_1 \sim \text{Categorical}(\pi)\).
  • State dynamics: the regimes follow a Markov chain with a \(3 \times 3\) transition matrix \(A\), \[ p(z_t = k \mid z_{t-1} = j) = A_{jk}, \qquad j, k \in \{1, 2, 3\}. \]
  • Emission model: the negative binomial of Problem 1 (b), with a mean \(\mu_k\) per state and a dispersion \(\phi\) shared by the states.

Notation as in lecture: \(\alpha_t(k) = p(y_{1:t}, z_t = k)\), \(\beta_t(k) = p(y_{t+1:T} \mid z_t = k)\), \(\gamma_t(k) = p(z_t = k \mid y_{1:T})\) and \(\xi_t(j, k) = p(z_{t-1} = j, z_t = k \mid y_{1:T})\). When a subject index is needed we write \(y^{(i)}_t\), \(\mu^{(i)}_k\) and \(A^{(i)}\), and \(T_i\) for the record length.

(a) Forward-backward with missing days (20 points)

On a day with no observation, the chain still moves but emits nothing. The forward and backward recursions go through unchanged if, on such a day, the emission probability \(B_t(k) = p(y_t \mid z_t = k)\) is replaced by \(1\) for every \(k\).

For the fits, initialize the rates at the \(20\)th, \(50\)th and \(80\)th percentiles of the observed counts and \(A\) at \(0.8\) on the diagonal, and order the states by rate.

Deliverables.

  1. (4 points) Show that replacing \(B_t(k)\) by \(1\) on a missing day gives the correct recursions, and say what \(\sum_k \alpha_{T}(k)\) is at the end of a record with missing days. Then say which of the M-step updates use the \(\gamma_t\) and \(\xi_t\) of a missing day and which must not, and why.

  2. (2 points) With \(\pi = (0.3, 0.2, 0.5)\), rates \((2000, 6000, 11000)\) and \[ A = \begin{pmatrix} 0.90 & 0.08 & 0.02 \\ 0.08 & 0.80 & 0.12 \\ 0.02 & 0.10 & 0.88 \end{pmatrix}, \] compute the log likelihood \(\log p(y_{1:T_i} \mid \theta)\) from the forward pass, under Poisson emissions and under negative binomial emissions with \(\phi = 0.08\), for subjects 1 and 2. The four values should be

    subject row of Y record Poisson negative binomial
    1 0 38 days \(-8203.66\) \(-312.75\)
    2 1 365 days \(-49469.78\) \(-3166.65\)

    Subject 1 tests that you took the right slice of the row; subject 2 tests the missing days inside it. The Poisson column only tests the recursion; every fit from here on uses the negative binomial.

  3. (5 points) Fit a negative binomial HMM with \(K = 3\) to each full-year subject, estimating \(\pi\), \(A\), \(\mu\) and \(\phi\), and tabulate \(\mu\), \(\phi\) and the diagonal of \(A\) for the first five subjects.

  4. (5 points) For each full-year subject, label each observed day by its most probable state, \(\arg\max_k \gamma_t(k)\), and put it beside the mixture label from Problem 1 (b) deliverable 4. Summed over subjects, give the fraction of observed days on which the two labels agree and the number of label changes between consecutive observed days under each model. Plot one subject over \(60\) days, counts on top and the two label sequences below.

  5. (4 points) Explain what information the HMM uses that the mixture throws away, and why that makes the HMM change label less often. Then use the entries of \(A\) to say how much log likelihood, in nats (the units of a natural log), the transition probabilities charge for leaving the high state for a single medium day and coming back.

(b) How many states? (15 points)

So far \(K = 3\) was handed to you. Nothing in the data announces it, and the likelihood cannot choose it on its own.

To judge a model on days it has not seen, make some of the full-year records short on purpose. Among the \(28\) full-year subjects, truncate every other one to its first \(28\) days and the rest to their first \(309\) days; the \(12\) genuinely short records stay as they are. From here on, every model is fitted on this training set and scored on the held-out days by the log predictive likelihood \[ \log p\big(y_{T_{\text{train}}+1:T} \mid y_{1:T_{\text{train}}}\big) = \log p(y_{1:T}) - \log p(y_{1:T_{\text{train}}}), \] both terms from a forward pass with the parameters fitted on the training days only. The \(12\) short records are part of the training set, so they help fit anything that is shared across subjects, but they have no held-out days and do not enter the scores. Use the same split for everyone:

full = np.where(n_days == T)[0]
cut28, cut309 = full[0::2], full[1::2]
n_train = n_days.copy()
n_train[cut28], n_train[cut309] = 28, 309
train = [Y[i, start[i]:start[i] + n_train[i]] for i in range(N)]   # the full-year subjects start at column 0
print(f"{len(cut28)} subjects cut to 28 days, {len(cut309)} cut to 309 days, {(n_days < T).sum()} genuinely short")
14 subjects cut to 28 days, 14 cut to 309 days, 12 genuinely short

For a \(K\)-state HMM, initialize the rates at the percentiles \(100\,(k - \tfrac12)/K\) of the observed counts and \(A\) at \(0.8\) on the diagonal. The Bayesian information criterion is \(\text{BIC} = -2 \log p(y_{\text{train}} \mid \hat\theta) + d_K \log n_{\text{obs}}\), with \(d_K\) the number of free parameters and \(n_{\text{obs}}\) the number of observed training days; compute it per subject, with that subject’s \(n_{\text{obs}}\), and sum over the group, as you do with the log likelihoods. For the mixture of Problem 1, days are independent, so the held-out log predictive likelihood is simply the log likelihood of the held-out observed days under the fitted mixture.

Deliverables.

  1. (8 points) For the \(309\)-day group: fit the negative binomial HMM on the training days with \(K = 2, 3, 4\), and the negative binomial mixture of Problem 1 with \(K = 3\). Make one table with a row per model and, summed over the group, the training log likelihood, BIC and the held-out log predictive likelihood. State \(d_K\) for each row.
  2. (4 points) Explain why the training log likelihood must rise with \(K\) for the HMM rows, say what BIC and the held-out likelihood each charge for, and which \(K\) each chooses.
  3. (3 points) The HMM and the mixture at \(K = 3\) have the same emission model; the only thing the HMM adds is the Markov chain, which lets yesterday’s regime inform today’s. The held-out column of your table says how much that is worth. Take the difference between the HMM row and the mixture row, in nats, and divide it by the number of held-out observed days in the group, so that you have the gain per day. Then say what this number means for a single new day: by what factor does knowing the chain multiply the probability the model assigns to that day’s count, and why it is not larger than that.

(c) A Dirichlet prior on the rows of \(A\) (15 points)

Each subject has so far had their own \(A^{(i)}\), estimated from their record alone. A \(28\)-day record has \(27\) transitions to estimate nine probabilities. This part borrows strength across subjects with a Dirichlet prior that pulls every subject toward one population matrix.

The lecture put a Dirichlet prior on each row of \(A\) and derived the maximum a posteriori (MAP) update, \[ A_{j, \cdot} \sim \text{Dirichlet}(\alpha_{j, \cdot}) \quad\Longrightarrow\quad A_{jk} \leftarrow \frac{\bar n_{jk} + \alpha_{jk} - 1}{\sum_{k'} \big(\bar n_{jk'} + \alpha_{jk'} - 1\big)}, \qquad \bar n_{jk} = \sum_{t \ge 2} \xi_t(j, k). \]

The population matrix. First fit a single transition matrix \(\bar A\) shared by all \(40\) subjects, with subject-specific emission means \(\mu^{(i)}_k\) and a single dispersion \(\phi\) shared by everyone. This is one EM in which the E-step runs forward-backward on each subject, the \(A\) update pools \(\bar n^{(i)}_{jk}\) across subjects, the \(\mu\) update is per subject, and the \(\phi\) and \(\pi\) updates pool across subjects. Write it as a function that takes the list of records and returns all parameters; part (d) builds on it.

Shrinkage toward the population. Then give each subject their own \(A^{(i)}\) with the prior \[ A^{(i)}_{j, \cdot} \sim \text{Dirichlet}\big(\kappa\, \bar A_{j, \cdot} + 1\big), \] so that the MAP update adds \(\kappa\, \bar A_{jk}\) pseudo-transitions to subject \(i\)’s expected counts; the \(+1\) makes the pseudo-counts exactly \(\kappa\bar A_{jk}\) and keeps every \(\alpha_{jk}\) above one. The concentration \(\kappa\) is the number of imaginary transitions of population behaviour each subject is assumed to have exhibited. At \(\kappa = 0\) the update is the per-subject maximum likelihood estimate. As \(\kappa \to \infty\) the pseudo-counts swamp the data and every subject’s matrix becomes \(\bar A\), so no new fit is needed for that case: the shared fit above is the \(\kappa = \infty\) row.

Deliverables.

  1. (3 points) The MAP update adds \(\alpha_{jk} - 1\) to the expected counts, while the posterior mean of the row given the counts would add \(\alpha_{jk}\). Derive both from the Dirichlet density and say what \(\alpha_{jk} = 1\) for all \(j, k\) does to the update.
  2. (3 points) Suppose that at some iteration \(A_{jk} = 0\) exactly. Show from the definition of \(\xi_t(j, k)\) that \(A_{jk}\) stays zero at every later iteration, and say which choice of \(\alpha\) makes this impossible.
  3. (3 points) Fit the shared model on the training set with \(K = 3\). Print \(\bar A\).
  4. (4 points) Refit with the shrinkage prior for \(\kappa \in \{0, 5, 20, 100, 1000\}\): the same pooled EM as the shared fit, except that each subject now has their own \(A^{(i)}\), updated by the MAP formula with \(\alpha = \kappa \bar A + 1\). Take the shared fit as \(\kappa = \infty\). Make one table with the held-out log predictive likelihood summed over the \(28\)-day group and, separately, over the \(309\)-day group, one row per \(\kappa\). Below it, print the estimated \(A^{(i)}\) of one subject in the \(28\)-day group at \(\kappa = 0\) and at \(\kappa = 20\).
  5. (2 points) Explain why \(\kappa\) matters so much for one group and so little for the other, and why the \(309\)-day group eventually gets worse as \(\kappa\) grows.

(d) A mixture over transition matrices (20 points)

The shrinkage prior assumes there is one way to behave and that subjects scatter around it. Suppose instead that there are \(M\) kinds of people. Each subject has a type \(c_i \in \{1, \ldots, M\}\) drawn from \(\text{Categorical}(w)\), and uses the type’s matrix: \(A^{(i)} = A_{c_i}\). The type is never observed. The model now has two layers of latent variables, a type per subject and a regime per day. There is no Dirichlet prior here: each \(A_m\) is estimated by maximum likelihood from the expected transition counts of all subjects, weighted by how much each subject belongs to type \(m\), so even a \(28\)-day record contributes to a matrix with hundreds of transitions behind it. The prior in part (c) was needed because a matrix estimated from one short record is not.

Shaded nodes are observed, small dots are parameters, and the two plates are disjoint: the matrices \(A_1, \ldots, A_M\) and the weights \(w\) belong to the population, everything inside the subject plate belongs to one record. Read the picture from the type down. Conditioned on \(c_i = m\), subject \(i\) is the ordinary HMM of part (a) with transition matrix \(A_m\), so forward-backward under \(A_m\) delivers everything below the type node. The only new quantity is the posterior over the type, and it is the mixture responsibility of Problem 1 with the component density replaced by the marginal likelihood of the whole record under \(A_m\), which is one forward pass. A mixture whose components are HMMs, in other words, and the same EM handles both layers. The structure of one iteration is

# E-step, for every subject i and every type m:
#     gamma[i, m], xi[i, m], ll[i, m] = forward_backward(logB[i], pi, A[m])
#     log r[i, m] = log w[m] + ll[i, m], normalized over m           (logsumexp)
#     gamma[i]    = sum over m of r[i, m] * gamma[i, m]
# M-step:
#     w[m]  = mean over i of r[i, m]
#     A[m]  = sum over i of r[i, m] * xi[i, m], each row normalized
#     mu[i], phi  from gamma[i], exactly as in part (a)
# observed data log likelihood of subject i: logsumexp over m of (log w[m] + ll[i, m])

Extend your population EM from part (c) to \(M\) types; \(M = 1\) is the shared fit you already have. The two types are two transition matrices for the regime chain \(z_t\), and the symmetric start \(A_1 = A_2\) is a fixed point of EM, so initialize them differently, with \(0.6\) and \(0.95\) on the diagonal and the rest of each row split equally: \[ A_1 = \begin{pmatrix} 0.6 & 0.2 & 0.2 \\ 0.2 & 0.6 & 0.2 \\ 0.2 & 0.2 & 0.6 \end{pmatrix}, \qquad A_2 = \begin{pmatrix} 0.95 & 0.025 & 0.025 \\ 0.025 & 0.95 & 0.025 \\ 0.025 & 0.025 & 0.95 \end{pmatrix}. \] That is, start with one type whose regimes last a few days and one whose regimes last weeks, and let EM sort the subjects between them.

Deliverables.

  1. (6 points) Write the complete data log likelihood of all \(N\) records, including the \(\log w_{c_i}\) terms. From it, derive the type responsibility \[ r_{im} = p(c_i = m \mid y^{(i)}) = \frac{w_m\, p(y^{(i)} \mid A_m, \mu^{(i)}, \phi)}{\sum_{m'} w_{m'}\, p(y^{(i)} \mid A_{m'}, \mu^{(i)}, \phi)} \] and the M-step updates for \(w_m\) and \(A_m\). Then write down the observed data log likelihood of one record under this model, which is also what the held-out score must use.
  2. (10 points) Fit \(M = 2\) on the training set. Print \(A_1\), \(A_2\) and \(w\). Give the held-out log predictive likelihood for the \(28\)-day and the \(309\)-day groups, to be read next to the \(\kappa\) table of part (c), and the mean over subjects of \(\max_m r_{im}\), separately for the \(309\)-day group, the \(28\)-day group and the \(12\) genuinely short records.
  3. (4 points) Say what kind of person \(\bar A\) from part (c) describes and whether anyone in the data behaves that way. Point to where in the two tables the mixture beats the best \(\kappa\) and where it does not, and say what the pattern of \(\max_m r_{im}\) across the three groups says about why.