Practical 8: Fitting a Hidden Markov Model

Anastasia Giachanou, Tina Shahedi

Machine Learning with Python - Utrecht Summer School

Welcome to the Hidden Markov Model practical!

This practical is built on Emmeke Aarts's practical workshop on "Extracting personalised latent dynamics using a multilevel hidden Markov model".

In this practical, we introduce a dataset and use it to fit a two-state hidden Markov model (HMM) (so not multilevel).

Explore the documentation provided in libraries like hmmlearn or pomegranate, which are used for building and analyzing hidden Markov models. The complete documentation for both can be accessed in bellow:

1. hmmlearn
2. pomegranate

For this practical, we will use several Python libraries. The hmmlearn library is a popular choice for building Hidden Markov Models in Python. We start by ensuring that it is installed in our environment.

Learning objectives:

By the end of this practical, you will be able to:

• Understand the core principles of Hidden Markov Models and how they apply to emotion data.
• Fit a two-state HMM using hmmlearn to model latent mood dynamics.
• Visualize emission distributions (probability of mood values per state) and interpret emotional states.
• Predict and analyze hidden state sequences across individuals.
• Compare emotional profiles and transitions across subjects using visualizations and summary statistics.

!pip install hmmlearn
!pip install plotly

Requirement already satisfied: hmmlearn in /usr/local/lib/python3.11/dist-packages (0.3.3)
Requirement already satisfied: numpy>=1.10 in /usr/local/lib/python3.11/dist-packages (from hmmlearn) (2.0.2)
Requirement already satisfied: scikit-learn!=0.22.0,>=0.16 in /usr/local/lib/python3.11/dist-packages (from hmmlearn) (1.6.1)
Requirement already satisfied: scipy>=0.19 in /usr/local/lib/python3.11/dist-packages (from hmmlearn) (1.15.3)
Requirement already satisfied: joblib>=1.2.0 in /usr/local/lib/python3.11/dist-packages (from scikit-learn!=0.22.0,>=0.16->hmmlearn) (1.5.1)
Requirement already satisfied: threadpoolctl>=3.1.0 in /usr/local/lib/python3.11/dist-packages (from scikit-learn!=0.22.0,>=0.16->hmmlearn) (3.6.0)
Requirement already satisfied: plotly in /usr/local/lib/python3.11/dist-packages (5.24.1)
Requirement already satisfied: tenacity>=6.2.0 in /usr/local/lib/python3.11/dist-packages (from plotly) (8.5.0)
Requirement already satisfied: packaging in /usr/local/lib/python3.11/dist-packages (from plotly) (25.0)

# Importing necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import hmmlearn
import plotly.graph_objects as go
from hmmlearn import hmm
from scipy.stats import norm

Data

We will be working with an open access dataset from Rowland and Wenzel (2020), which is part of a study involving 125 undergraduate students from the University of Mainz in Germany. These students completed a 40-day ambulatory assessment six times a day, reporting on their affective experiences such as happiness, excitement, relaxation, satisfaction, anger, anxiety, depression, and sadness. These affective states were quantified using a visual analog slider, ranging from 0 to 100.

Before the data collection, participants were randomly assigned to either a group receiving weekly mindfulness treatment during the study or a control group. We will be working with this dataset dataset. This dataset has been cleaned and is provided by Haslbeck, Ryan, & Dablander (2023) and can be found in their OSF repository.

1. Load the dataset ('emotion_data.csv') that can be found in the folder with the practicals' datasets. You can use the function pd.read_csv(). Give the dataset the name 'emotion_data'. Then, inspect the first few rows of the data.

emotion_data = pd.read_csv("emotion_data.csv")
emotion_data.head()

	subj_id	dayno	beep	group	happy	excited	relaxed	satisfied	angry	anxious	depressed	sad
0	1	1	1	2	79.0	64.0	30.0	74.0	19.0	7.0	8.0	58.0
1	1	1	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	1	1	3	2	20.0	24.0	37.0	36.0	62.0	41.0	73.0	67.0
3	1	1	4	2	25.0	24.0	33.0	29.0	42.0	68.0	77.0	81.0
4	1	1	5	2	37.0	24.0	40.0	17.0	35.0	36.0	12.0	76.0

The dataset contains multiple affective states recorded over time for multiple participants. Each row in the dataset represents a unique measurement instance. Before we fit a hidden Markov model, we need to understand the structure of the data, including the variables and their relationships.

emotion_data.describe()

	subj_id	dayno	beep	group	happy	excited	relaxed	satisfied	angry	anxious	depressed	sad
count	30000.000000	30000.000000	30000.000000	30000.000000	21570.000000	21570.000000	21570.000000	21570.000000	21570.000000	21570.000000	21570.000000	21570.000000
mean	80.736000	20.500000	3.500000	1.488000	59.333380	40.937089	54.761196	54.989754	13.508067	13.200649	20.590264	15.969866
std	46.926961	11.543589	1.707854	0.499864	24.213283	28.544285	26.221243	26.250658	20.935787	21.287047	24.603760	22.706737
min	1.000000	1.000000	1.000000	1.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	37.000000	10.750000	2.000000	1.000000	39.000000	18.000000	33.000000	33.000000	0.000000	0.000000	0.000000	0.000000
50%	84.000000	20.500000	3.500000	1.000000	64.000000	35.000000	61.000000	61.000000	4.000000	3.000000	11.000000	5.000000
75%	121.000000	30.250000	5.000000	2.000000	77.000000	66.000000	74.000000	74.750000	19.000000	18.000000	31.000000	24.000000
max	164.000000	40.000000	6.000000	2.000000	100.000000	100.000000	100.000000	100.000000	100.000000	100.000000	100.000000	100.000000

we can see the summary of the features here:

No	Name	Label	Options
1	subjno	individual code	1-164
2	dayno	study days	1-40
3	beep	daily signals	1-6
4	group	condition allocation	1 = control, 2 = training
5	emo1_m	happy	0-100
6	emo2_m	excited	0-100
7	emo3_m	relaxed	0-100
8	emo4_m	satisfied	0-100
9	emo5_m	angry	0-100
10	emo6_m	anxious	0-100
11	emo7_m	depressed	0-100
12	emo8_m	sad	0-100

Next, we need to preprocess the data to make it suitable for fitting a hidden Markov model. First we are gonna check for missing values across the columns by applying the .isnull().sum() function to each column in the dataset:

# Count of NaN values in each column
NaN_counts = emotion_data.isnull().sum()
NaN_counts

	0
subj_id	0
dayno	0
beep	0
group	0
happy	8430
excited	8430
relaxed	8430
satisfied	8430
angry	8430
anxious	8430
depressed	8430
sad	8430

dtype: int64

For this case, we will fill the NaN value by mean.

HMMs model temporal sequences where the order and alignment of observations matter. Removing rows with NaNs would disrupt the time series continuity, making the sequence shorter and potentially meaningless. Also, the hmmlearn library expects a full matrix of observations. If any row has a missing value, it will raise an error. Mean imputation ensures the input stays complete.

Mean imputation is a simple, fast, and model-agnostic approach that preserves the overall distribution of the variable—though it's not always optimal.

# Remove rows with any NaN values
emotion_data_clean = emotion_data.fillna(emotion_data.mean())

Now, we will create a plot to visualize the hidden states over time for each subject in the dataset.

2. Visualize time sequences for individual subjects separately to analyze emotional responses over time. Start by transforming your dataset from wide to long format using melt() function from pandas, which prepares it for detailed analysis. Select only five emotions: 'happy', 'excited', 'relaxed', 'angry', 'depressed'. How does the dataframe look like now?

pandas.melt() transforms a DataFrame from wide format to long format.

• Wide format: Each emotion is in its own column.
• Long format: Emotions are stored in a single column (variable), and their values are stored in another (value).

Example:

Before melting (wide format):

id	time	happy	excited	angry
1	1	0.5	0.8	0.2
1	2	0.4	0.6	0.3

After melting (long format):

pd.melt(df, id_vars=['id', 'time'], value_vars=['happy', 'excited', 'angry'])

id	time	variable	value
1	1	happy	0.5
1	1	excited	0.8
1	1	angry	0.2
1	2	happy	0.4
1	2	excited	0.6
1	2	angry	0.3

Why we do this?

We want to plot how each emotion changes over time. To do this with tools like plotly, the data is better to be in long format — so each row is one observation of an emotion at a specific time point.

# Creating a new column to represent the beep normalized over days
emotion_data_clean.loc[:, 'beep2'] = (emotion_data['dayno'] - 1) * 6 + emotion_data['beep']

# Reshaping the data for visualization: melt the DataFrame to long format
emo_long = emotion_data_clean.melt(id_vars=['subj_id', 'beep2'],
                                   value_vars=['happy', 'excited', 'relaxed', 'angry', 'depressed'],
                                   var_name='emotion', value_name='outcome')

emo_long.head()

	subj_id	beep2	emotion	outcome
0	1	1	happy	79.00000
1	1	2	happy	59.33338
2	1	3	happy	20.00000
3	1	4	happy	25.00000
4	1	5	happy	37.00000

3. Create a set of individual line plots to visualize emotional responses over time for a small group of subjects. Use seaborn.FacetGrid to generate separate plots for the first four subjects, making it easier to interpret each subject’s emotional patterns without clutter.

# Filter for specific subjects to avoid a crowded plot (e.g., first 4 subjects)
emo_filtered = emo_long[emo_long['subj_id'].isin([1, 2, 3, 4])]

# The next line initializes a FacetGrid from Seaborn:
# col='subj_id': creates one subplot per subject.
# col_wrap=1: stacks all plots vertically.
# height=3, aspect=3: controls the size and shape of each subplot for better readability.
g = sns.FacetGrid(emo_filtered, col='subj_id', col_wrap=1, height=3, aspect=3)

# The next line maps a line plot to each facet: x='beep2': time points (beeps), y='outcome': emotion score,
# hue='emotion': different emotions shown in different colors, estimator=None: plots raw values instead of averaging them.
g.map_dataframe(sns.lineplot, x='beep2', y='outcome', hue='emotion', estimator=None)

# Enhancing the plot
g.add_legend()
g.set_titles("Subject {col_name}")
g.set_axis_labels("Beep Time Point", "Emotional Outcome")
plt.show()

Setting Up the Hidden Markov Model

In this section, you will fit a 2-state hidden Markov model. For this, you need a DataFrame in which only the affective variables you want to use in the model are included, plus the subject ID as the first column. Please feel free to use any subset of the provided affective variables, but do make sure to select at least two affective variables.

4. Create a dataframe emotion_mHMM which has the subject ID variable in the first column, followed by only the affective variables you want to use in the hidden Markov model. We are using 'happy', 'excited', 'relaxed', 'angry', 'depressed'

# Create the emotion_mHMM dataset
emotion_mHMM = emotion_data_clean[['subj_id', 'happy', 'excited', 'relaxed', 'angry', 'depressed']]

# Inspect the first few rows
print(emotion_mHMM.head())

   subj_id     happy    excited    relaxed      angry  depressed
0        1  79.00000  64.000000  30.000000  19.000000   8.000000
1        1  59.33338  40.937089  54.761196  13.508067  20.590264
2        1  20.00000  24.000000  37.000000  62.000000  73.000000
3        1  25.00000  24.000000  33.000000  42.000000  77.000000
4        1  37.00000  24.000000  40.000000  35.000000  12.000000

5. Set up the first set of model input arguments: m =2 for the number of states, n_dep for the number of dependent variables (in our case this will be 5), and starting values for gamma (this is the initial state transition matrix and can be defined with an np.array) and the emission distribution (that refer to the starting values for the Gaussian emission distributions for each of the number of selected variables).

# Define the number of states and dependent variables
m = 2
n_dep = 5

# Define starting values for the transition matrix (gamma)
# The matrix is 2x2 because we have 2 hidden states.
# Each row represents the probability of moving from one state to another.

#start_gamma[0] = [0.7, 0.3] means:
#If you're in state 0, there's a 70% chance you'll stay in state 0, and a 30% chance you'll switch to state 1.
start_gamma = np.array([[0.7, 0.3], [0.3, 0.7]])

# Define starting values for the emission distribution
# Each np.array is a 2×2 matrix where row 0 is for hidden state 0, row 1 is for hidden state 1.
# First column = mean, second column = standard deviation (SD) of the Gaussian.
start_emiss = [
    np.array([[70, 15], [30, 15]]),  # happy
    np.array([[70, 15], [30, 15]]),  # excited
    np.array([[70, 15], [30, 15]]),  # relaxed
    np.array([[15, 10], [40, 15]]),  # angry
    np.array([[15, 10], [40, 15]])   # sad
]

Next, we convert the data into a format required by hmmlearn

hmmlearn expects the format: means_: shape (n_components, nfeatures) → (2 states × 5 variables) `covars`: same shape, variances

So we are going to extract from our data the following:

• Means per state → start_means
• Variances per state → start_covars (note: variance = std²)

# Convert start_emiss to the format required by hmmlearn
start_means = np.array([emiss[:, 0] for emiss in start_emiss]).T
start_covars = np.array([emiss[:, 1] for emiss in start_emiss]).T ** 2

We also drop the subj_id column from the emotion_mHMM DataFrame (HMMs learn from the sequence of observed features, not the identity.) to ensure the observations matrix consists only of emotional state data for Hidden Markov Model (HMM) analysis.

# Prepare the data for HMM excluding the subject ID
observations = emotion_mHMM.drop(columns=['subj_id']).values
# Ensure the shape of observations
print(f"Shape of observations: {observations.shape}")  # Should be (n_samples, n_dep)

Shape of observations: (30000, 5)

Fitting a 2-State Hidden Markov Model with Custom Initialization

6. Initialize a Gaussian Hidden Markov Model (hmm.GaussianHMM) with a n_component = m (m was already set up to 2) and n_iter=500. Set the model's initial state probabilities (np.array([0.5, 0.5])), transition matrix, means, and covariances using predefined variables (start_gamma, start_means, start_covars)

# Initialize the HMM model without initializing the parameters ('stmc')
model = hmm.GaussianHMM(n_components=m, n_iter=500, init_params='')

# Set initial parameters
model.startprob_ = np.array([0.5, 0.5])  # Example initial state probabilities
model.transmat_ = start_gamma
model.means_ = start_means
model.covars_ = start_covars

Next, we will proceed to fit a Hidden Markov Model to this preprocessed data.

7. Fit (guess how the function is called... yes, right it is called fit()) a 2-state hidden Markov model. So use the fit method of the HMM model to train it on the prepared observations data. Assign the fitted model to a variable named out_2st_emotion.

# Fit the HMM model
model.fit(observations)

# Save the fitted model
out_2st_emotion = model

print("Model fitting completed.")

Model fitting completed.

Inspecting General Model Output

8. Let’s evaluate the fitted Hidden Markov Model using two key metrics. First, count how many unique individuals are included in the dataset (use .nunique() on subj_id). Second, compute the overall average log-likelihood of the model across all data points (model.score()). Print the number of subjects and the average log-likelihood to summarize key metrics These metrics will help you assess who the model is fitted on and how well it explains the observed emotional data.

# Number of subjects
num_subjects = emotion_data_clean['subj_id'].nunique()

# compute the log-likelihood of the entire dataset under the current model — a higher (less negative) value means a better fit.
avg_log_likelihood = model.score(observations)
print(f"Number of subjects: {num_subjects}")
print(f"Average Log likelihood over all subjects: {avg_log_likelihood}")

Number of subjects: 125
Average Log likelihood over all subjects: -637693.4192899191

We can see that there are 25 subjects, our model was trained on emotional data from 125 individuals. Also, the log-likelihood is -637,693: This is the log probability of the observed data given the model. It’s negative because it’s a log of a small probability. Log-likelihood doesn’t tell us if a model is "good" in isolation. But we can use it to compare models with different numbers of states

9. Document the structure of your Hidden Markov Model. First, use the n_components attribute to find out how many hidden states were used. Next, use the .shape attribute of the observations array to determine how many affective (dependent) variables were included. Print both values to summarize the model’s configuration.

num_states = model.n_components
num_dependent_vars = observations.shape[1]
print(f"Number of states used: {num_states}")
print(f"Number of dependent variables used: {num_dependent_vars}")

Number of states used: 2
Number of dependent variables used: 5

10. Let’s explore what the Hidden Markov Model has learned. Print the state transition probability matrix using model.transmat_ to understand how likely it is to switch or stay in the same state. Store the state transition matrix from the fitted HMM model in a variable called gamma.

# model.transmat_: Shows the probabilities of transitioning from one state to another.
gamma = model.transmat_
print("\nState transition probability matrix:")
print(gamma)

State transition probability matrix:
[[0.85621468 0.14378532]
 [0.38117036 0.61882964]]

We can have the following conclusions:

• State 1 → State 1: 85.6% chance of staying in State 1
• State 1 → State 2: 14.4% chance of switching to State 2
• State 2 → State 2: 61.9% chance of staying in State 2
• State 2 → State 1: 38.1% chance of switching to State 1

Interpretation:

• State 1 is more stable and persistent (likely to remain in that mood state).
• State 2 is less stable and more likely to switch back to State 1.
  This suggests that positive mood states (State 1) tend to last longer than negative mood states (State 2).

11. Now, we will access the model’s learned parameters for each emotion in both hidden states. Use the model.means_ and model.covars_ attributes from the model to extract and save these parameters. Organize the emission statistics into a dictionary called emiss, so each emotion shows its mean and standard deviation in both states. Print the emiss to inspect the mean and standard deviation for each dependent variable in each state. Which are your conclusions? Which emotions have the most distinct means between State 1 and State 2?

# Access the means and standard deviations (from covariances)
# model.means_: Contains the mean of each affective variable in each hidden state → shape (2, 5)
# model.covars_: Contains variances, so we take the square root to get standard deviations
means = model.means_
covariances = np.sqrt(np.diagonal(model.covars_, axis1=1, axis2=2))  # Extract standard deviations

# Save emission distribution parameters
means = model.means_
covariances = np.sqrt(np.diagonal(model.covars_, axis1=1, axis2=2))  # Extract standard deviations

emiss= {emotion: {'State 1 Mean': means[0][i], 'State 1 SD': covariances[0][i],
                         'State 2 Mean': means[1][i], 'State 2 SD': covariances[1][i]}
               for i, emotion in enumerate(['happy', 'excited', 'relaxed', 'angry', 'depressed'])}

# Inspect the emission distribution parameters
print("\nEmission distribution (continuous) for each of the dependent variables:")
for emotion, stats in emiss.items():
    print(f"\n${emotion.capitalize()}")
    for state, values in stats.items():
        print(f"{state}: {values:.3f}")

Emission distribution (continuous) for each of the dependent variables:

$Happy
State 1 Mean: 66.487
State 1 SD: 14.231
State 2 Mean: 40.370
State 2 SD: 22.513

$Excited
State 1 Mean: 46.179
State 1 SD: 23.239
State 2 Mean: 27.040
State 2 SD: 21.001

$Relaxed
State 1 Mean: 61.398
State 1 SD: 17.772
State 2 Mean: 37.165
State 2 SD: 23.263

$Angry
State 1 Mean: 8.085
State 1 SD: 7.099
State 2 Mean: 27.886
State 2 SD: 27.061

$Depressed
State 1 Mean: 12.970
State 1 SD: 10.159
State 2 Mean: 40.793
State 2 SD: 27.446

---

Emission Distributions by State

Emotion	State 1 (Mean ± SD)	State 2 (Mean ± SD)
Happy	66.5 ± 14.2	40.4 ± 22.5
Excited	46.2 ± 23.2	27.0 ± 21.0
Relaxed	61.4 ± 17.8	37.2 ± 23.3
Angry	8.1 ± 7.1	27.9 ± 27.1
Depressed	13.0 ± 10.2	40.8 ± 27.4

---

Interpretation

• State 1 = Positive Mood:
  High levels of happy, excited, and relaxed; low levels of angry and depressed.

• State 2 = Moderate to Negative Mood:
  Happiness and excitement are reduced compared to State 1, while negative emotions like anger and depression are much higher.

The HMM successfully identifies two meaningful emotional states that align with intuitive psychological mood categories.

Which emotions have the most distinct means between State 1 and State 2?

📊 Difference in Emotion Means Between States

Emotion	State 1 Mean	State 2 Mean	Difference
Happy	66.49	40.37	26.12
Excited	46.18	27.04	19.14
Relaxed	61.40	37.17	24.23
Angry	8.09	27.89	19.80
Depressed	12.97	40.79	27.82

Most distinct: Depressed, Happy, and Relaxed show the largest mean differences. These variables help most clearly separate the emotional profiles of the two states.

Visualizing the Obtained Output

Visualizing the transition probabilities can be very helpful, especially when dealing with a large number of states and/or dependent variables.

12. In this question, you will create a Sankey diagram to visually represent the transition probabilities between hidden emotional states in your HMM. The Sankey diagram will help you understand how frequently the model predicts staying in the same state vs switching states. Visually compare transitions from State 1 and State 2. Your diagram should include: 1. Nodes for each state (e.g., State 1, State 2 as both source and target), 2. Arrows (links) representing transitions from and to each state. 3. Arrow thickness scaled to transition probability, 4. Labels and colors to distinguish directions (e.g., blue for staying, red for switching)

We give you a structure to help you with that.

# Step 1: Define state labels and number of states
states = [...]  # e.g., ['State 1', 'State 2']
num_states = ...

# Step 2: Prepare labels for source and target nodes
labels = ...
source = ...  # from-state indices
target = ...  # to-state indices
values = ...  # flatten the transition matrix (e.g., gamma.flatten())

# Step 3: Create the Sankey diagram
fig = go.Figure(data=[go.Sankey(
    node=dict(
        label=...,
        color=...,
        # padding, thickness, and line style (optional)
    ),
    link=dict(
        source=...,
        target=...,
        value=...,
        color=...  # Optional: different colors for stay vs switch
    )
)])

# Step 4: Customize layout
fig.update_layout(title_text="...", font_size=..., width=..., height=...)
fig.show()

# Prepare the data for the Sankey diagram
states = ['State 2', 'State 1']
num_states = len(states)

labels = states + states  # From and To states will have the same labels
source = np.repeat(np.arange(num_states), num_states)
target = np.tile(np.arange(num_states), num_states)
values = gamma.flatten()

# Create the Sankey diagram
fig = go.Figure(data=[go.Sankey(
    node=dict(
        pad=15,
        thickness=20,
        line=dict(color="black", width=0.5),
        label=labels,
        color=['white', 'white','white', 'white'] #Custom colors for nodes
    ),
    link=dict(
        source=source,
        target=target + num_states,  # Target states indices should start after source states
        value=values,
        color=['lightblue','lightblue','tomato','tomato']
    ))])

# Set plot title and layout
fig.update_layout(title_text="Transition Probabilities", font_size=14, title_x=0.5, width=700, height=500)

fig.show()

13. Visualize the transition probabilities gamma saved in the gamma object. Use a heatmap to represent the probabilities of transitioning from one state to another.

# Visualize the transition probabilities
plt.figure(figsize=(6, 4))
sns.heatmap(gamma, annot=True, cmap='viridis', cbar=True, xticklabels=['To state 1', 'To state 2'], yticklabels=['From state 1', 'From state 2'])
plt.title('State Transition Probability Matrix')
plt.xlabel('To State')
plt.ylabel('From State')
plt.show()

14. Visualize the emission distributions to understand the mean values of different dependent variables across states. Prepare the data by creating a DataFrame with the state, mean emission value, and dependent variable name using the emiss dictionary. Use Seaborn's barplot function to create the plot

# Prepare the data for visualization
emotions = ['happy', 'excited', 'relaxed', 'angry', 'depressed']
emiss_melt = pd.DataFrame({
    'State': np.repeat([1, 2], len(emotions)),
    'Mean': [emiss[emotion]['State 1 Mean'] for emotion in emotions] + [emiss[emotion]['State 2 Mean'] for emotion in emotions],
    'Dep': np.tile(emotions, 2)
})

# Visualize the emission distributions using bar plot
plt.figure(figsize=(10, 6))
sns.barplot(x='Dep', y='Mean', hue='State', data=emiss_melt, palette='Set3')
plt.title('Emission Distributions by State and Dependent Variable')
plt.xlabel('Mood State')
plt.ylabel('Mean Emission')
plt.legend(title='State')
plt.show()

15. Compute the density probabilities for each emotion and state. First define the length of the grid and generate a sequence of mood values ranging from 1 to 100. Select only those emotions that are present in the emiss dataset to ensure the accuracy of the density plots. Create a DataFrame emiss_dens with columns for state, emotion (Dep), mood value (X), and probability (Prob).Finally by using the normal distribution parameters calculate the density probabilitiesfor each emotion-state pair

# Prepare the data for density plot
len_grid = 200
grid = np.linspace(1, 100, len_grid)

# Ensure we only use emotions present in
emotions = [emotion for emotion in ['happy', 'excited', 'relaxed', 'angry', 'anxious', 'depressed', 'sad'] if emotion in emiss]

emiss_dens = pd.DataFrame({
    'State': np.repeat(np.repeat([1, 2], len_grid), len(emotions)),
    'Dep': np.tile(np.repeat(emotions, len_grid), 2),
    'X': np.tile(grid, len(emotions) * 2),
    'Prob': np.nan
})

for emotion in emotions:
    for i in range(1, 3):
        mask = (emiss_dens['Dep'] == emotion) & (emiss_dens['State'] == i)
        emiss_dens.loc[mask, 'Prob'] = norm.pdf(grid, loc=emiss[emotion][f'State {i} Mean'], scale=emiss[emotion][f'State {i} SD'])

emiss_dens.head()

	State	Dep	X	Prob
0	1	happy	1.000000	7.067787e-07
1	1	happy	1.497487	8.296296e-07
2	1	happy	1.994975	9.726449e-07
3	1	happy	2.492462	1.138921e-06
4	1	happy	2.989950	1.331994e-06

16. Create subplots for each emotion using seaborn.FacetGrid with specified attributes. Plot density probabilities for mood values differentiated by state using sns.lineplot. Add vertical dashed lines at the mean mood values for each state. Set y-axis limits from 0 to 0.05 for consistency, and include a legend and axis labels for mood value and density.

# Visualize the emission distributions using density plot
g = sns.FacetGrid(emiss_dens, col="Dep", col_wrap=3, height=4, aspect=1.2, hue='State', palette='Set1')
g.map(sns.lineplot, 'X', 'Prob')

# Add vertical lines for the means
for ax, emotion in zip(g.axes.flat, emotions):
    for i, state in enumerate([1, 2]):
        mean_val = emiss[emotion][f'State {state} Mean']
        ax.axvline(mean_val, color=g._colors[i], linestyle='--')

# Set the y-axis limit from 0 to 0.05
for ax in g.axes.flat:
    ax.set_ylim(0, 0.05)

g.add_legend(title='State')
g.set_axis_labels('Mood Value', 'Density')
plt.show()

From this plot, we see that for all emotions, the red and blue lines are clearly separated, meaning each state has a distinct mood distribution. In happy, State 1 (red) has a higher mood value than State 2 (blue), suggesting more intense happiness in that state. We also observe that wxcited, happy, and relaxed show overlapping but distinct bell curves, indicating similar variance with different means. Angry and depressed in State 1 have sharp peaks near the lower end, suggesting intense negative emotion, while State 2 has broader, flatter distributions with higher average mood values.

Across all emotions, State 1 appears to reflect more extreme mood values (either high or low), while State 2 shows more moderate values.

Obtaining the Most Likely Sequence of Hidden States

17. Using the predict method from hmmlearn, the fitted model, and the observations emotion_mHMM, obtain the most likely hidden state sequence. Save the state sequences in the object emotion_states_2st, and inspect the object.

# Obtain the most likely hidden state sequence
emotion_states_2st = model.predict(observations)

# Create a DataFrame to store the results
emotion_states_2st_df = pd.DataFrame({
    'subj': emotion_data_clean['subj_id'],
    'state': emotion_states_2st
})

# Inspect the first few entries
print(emotion_states_2st_df.head())

# Create a table to inspect the state distribution across subjects
state_distribution = pd.crosstab(emotion_states_2st_df['subj'], emotion_states_2st_df['state'])
print(state_distribution)

   subj  state
0     1      0
1     1      0
2     1      1
3     1      1
4     1      1
state    0    1
subj           
1      158   82
2      148   92
3      217   23
4      109  131
6      215   25
...    ...  ...
156    209   31
159    128  112
161    159   81
163     81  159
164    212   28

[125 rows x 2 columns]

We see a contingency table: each row is a subject, and the columns show how many of their observations fall into State 0 and State 1. For example, subject 1 has 158 observations in State 0 and 82 in State 1.

Plotting the Inferred State Sequence Over Time

In this section we will plot the inferred state sequence over time for a selection of the subjects within the data. To plot the state sequences over time, a new variable denoting the beep moment needs to be added to the matrix containing the inferred states.

18. Now, create a DataFrame that includes subj_id, state, and beep moment. Filter the DataFrame to include subjects 1 to 10, and map the states to categorical labels ('state 1' and 'state 2'). Use Seaborn's 'Set2' palette for the states and FacetGrid to create horizontal bar plots for each subject, adjusting the figure size and aspect ratio for better readability. Customize the plot with axis labels, titles, and a legend.

# Create a DataFrame to store the results including the beep moment
emotion_states_2st_df = pd.DataFrame({
    'subj': emotion_data_clean['subj_id'],
    'state':emotion_states_2st,
    'beep': emotion_data_clean['beep2']
})

# Filter for a selection of subjects (e.g., subjects 1 to 10)
selected_subjects = emotion_states_2st_df[emotion_states_2st_df['subj'].isin(range(1, 11))]

# Convert the state to a categorical variable with meaningful labels
selected_subjects.loc[:, 'state'] = selected_subjects['state'].map({0: 'state 1', 1: 'state 2'})

# Define the palette
palette = sns.color_palette('Set2')

# Plot the state sequences over time using horizontal bars
plt.figure(figsize=(10, 15))
g = sns.FacetGrid(selected_subjects, col='subj', col_wrap=1, height=1, aspect=10, sharey=False)
g.map_dataframe(sns.histplot, x='beep', hue='state', multiple='stack', bins=500, palette='Set2', discrete=True)

# Customize the plot
g.set_axis_labels("Beep", "")
g.set_titles("Subject {col_name}")
for ax in g.axes.flatten():
    ax.set_yticks([])
    ax.set_ylabel('')

# Add a legend
g.add_legend(title='State')

# Customize the legend
state_labels = ['state 1', 'state 2']
handles = [plt.Line2D([0], [0], color=color, lw=6) for color in palette[:2]]
plt.legend(handles, state_labels, title='State', loc='upper right', bbox_to_anchor=(1.15, 1))

plt.show()

/tmp/ipython-input-30-3112736098.py:12: FutureWarning:

Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '['state 1' 'state 1' 'state 2' ... 'state 1' 'state 1' 'state 1']' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.

<Figure size 1000x1500 with 0 Axes>

This plot visualizes the inferred hidden state sequences over time for several individual subjects (from a Hidden Markov Model with two states). It shows how each subject transitions between state 1 and state 2 over a series of timepoints or “beeps.”

• X-axis (Beep): Each tick represents a timepoint where the subject’s emotional state was inferred.
• Y-axis: Each row corresponds to one subject (e.g., Subject 1, Subject 2, etc.).
• Bars (Color-coded by state):
  • Green (state 1): Indicates timepoints where the model assigned the subject to State 1.
  • Red/Orange (state 2): Indicates timepoints in State 2.

Some subjects (e.g., Subject 3, Subject 10) spend most of their time in State 1. Others (e.g., Subject 8) spend a large portion in State 2. Subject 4 and Subject 7 frequently switch between states, which indicates more dynamic or unstable patterns.

Many subjects have blocks of the same color, showing that the model detects state persistence (i.e., once in a state, the subject tends to stay in it for a while before switching). Potential Behavioral/Emotional Interpretation: From the plot we can interpret how often and how long each subject experiences each state.

End of Practical!