lab-2 Report
Abubakar Passum Abdul Gaffar
2025-09-18
Task 1:
The ring isdivided into 10 sectors. At any given time point, the robot is in one of the sectors and decides with equal probability to stay in that sector or move to the next sector. You do not have direct observation of the robot. However, the robot is equipped with a tracking device that you can access. The device is not very accurate though: If the robot is in the sector i, then the device will report that the robot is in the sectors [i − 2;i + 2] with equal probability.
- Build a hidden Markov model (HMM) for the scenario described above. Note that the documentation for the initHMM function says that the emissionProbs matrix should be of size [number of states]x[number of states]. This is wrong. It should be of size [number of states]x[number of symbols]. The package works correctly, though. It is just the documentation that is wrong.
Solution
# Load required packages
library(HMM)
library(entropy)
# Define the states and symbols (sectors 1 to 10)
states =as.character(1:10)
symbols =as.character(1:10)
# Initialize the transition probability matrix
transProbs = matrix(0, nrow = 10, ncol = 10, dimnames = list(states, symbols))
for(i in 1:10){
# Probability of staying in the current sector
transProbs[i, i] =0.5
# Probability of moving to the next sector (with wrap-around)
#next_sector =ifelse(i == 10, 1, i + 1) #ifelse(test, yes, no) is a vectorized conditional function.
# More generalized is to use modulo
next_sector = (i %% 10) + 1
transProbs[i, next_sector] = 0.5
}
# Initialize the emission probability matrix
emissionProbs = matrix(0, nrow = 10, ncol = 10, dimnames = list(states, symbols))
for (i in 1:10) {
# For each state i, the device reports i-2, i-1, i, i+1, i+2 with equal probability (0.2)
# Handle wrap-around for the ring
indices = ((i - 2):(i + 2)) %% 10
indices[indices == 0] = 10 # Convert 0 to 10
indices[indices < 1] = indices[indices = 1] + 10 # Ensure positive indices
for (j in indices) {
emissionProbs[i, j] = emissionProbs[i, j] + 0.2
}
}
# Uniform start probabilities
startProbs =rep(1/10, 10)
names(startProbs) = states
# Initialize HMM
hmm =initHMM(States = states, Symbols = symbols,
startProbs = startProbs,
transProbs = transProbs,
emissionProbs = emissionProbs)
print("HMM initialized successfully.")## [1] "HMM initialized successfully."Results:
HMM object hmm created with 10 states and 10 symbols. transProbs ensures 50% chance to stay in the current sector and 50% to move to the next (e.g., from 10 to 1). emissionProbs assigns 0.2 probability to each of the 5 adjacent sectors (including current) with modular wrap-around.
Task 2
Simulate the HMM for 100 time steps.
Solution
# (2) Simulate the HMM for 100 time steps
set.seed(12345) # For reproducibility
simulation =simHMM(hmm, 100)
observations =simulation$observation
true_states =simulation$states
result <- data.frame(
Time = 1:100,
States = true_states,
Observations = observations
)
result## Time States Observations
## 1 1 9 7
## 2 2 9 10
## 3 3 9 8
## 4 4 9 10
## 5 5 10 2
## 6 6 1 3
## 7 7 2 10
## 8 8 2 3
## 9 9 2 4
## 10 10 2 4
## 11 11 3 5
## 12 12 3 4
## 13 13 4 2
## 14 14 4 3
## 15 15 4 2
## 16 16 4 6
## 17 17 4 6
## 18 18 4 5
## 19 19 4 4
## 20 20 5 3
## 21 21 6 5
## 22 22 6 5
## 23 23 7 8
## 24 24 8 9
## 25 25 9 10
## 26 26 10 9
## 27 27 10 9
## 28 28 10 10
## 29 29 1 2
## 30 30 2 10
## 31 31 2 2
## 32 32 3 5
## 33 33 3 3
## 34 34 4 2
## 35 35 4 6
## 36 36 4 6
## 37 37 5 4
## 38 38 5 7
## 39 39 5 7
## 40 40 6 6
## 41 41 7 5
## 42 42 7 9
## 43 43 8 7
## 44 44 9 10
## 45 45 10 10
## 46 46 1 3
## 47 47 2 3
## 48 48 3 1
## 49 49 3 3
## 50 50 4 3
## 51 51 5 6
## 52 52 5 5
## 53 53 6 4
## 54 54 6 7
## 55 55 7 7
## 56 56 7 9
## 57 57 8 9
## 58 58 8 10
## 59 59 8 6
## 60 60 8 9
## 61 61 9 10
## 62 62 10 2
## 63 63 10 9
## 64 64 10 9
## 65 65 10 8
## 66 66 1 1
## 67 67 1 3
## 68 68 2 2
## 69 69 2 3
## 70 70 2 4
## 71 71 2 3
## 72 72 2 2
## 73 73 3 5
## 74 74 3 4
## 75 75 3 4
## 76 76 4 2
## 77 77 5 4
## 78 78 5 6
## 79 79 5 4
## 80 80 6 6
## 81 81 7 8
## 82 82 8 10
## 83 83 8 8
## 84 84 8 7
## 85 85 8 6
## 86 86 8 6
## 87 87 9 7
## 88 88 9 8
## 89 89 9 9
## 90 90 10 10
## 91 91 10 1
## 92 92 10 9
## 93 93 1 2
## 94 94 1 2
## 95 95 1 3
## 96 96 1 9
## 97 97 1 2
## 98 98 1 10
## 99 99 2 4
## 100 100 3 1# Print the simulation results
#cat("True States (first 20 and last 20):\n")
#print(head(true_states, 20))
#print(tail(true_states, 20))
#cat("\nObservations (first 20 and last 20):\n")
#print(head(observations, 20))
#print(tail(observations, 20))States represent the robot’s true sector, while observations reflect the noisy tracking data, often differing from states due to the ±2 sector error range.
Task 3
Discard the hidden states from the sample obtained above. Use the remaining observations to compute the filtered and smoothed probability distributions for each of the 100 time points. Compute also the most probable path.
Solution
# (3) Compute filtered, smoothed distributions and Viterbi path
# Filtered probabilities (convert from log scale and normalize)
filtered_log_probs =forward(hmm, observations)
filtered_probs =apply(filtered_log_probs, 2, function(x) {
exp_x =exp(x - max(x)) # Numerical stability
prop.table(exp_x)
})
# Smoothed probabilities
smoothed_probs =posterior(hmm, observations)
# Viterbi path
viterbi_path =viterbi(hmm, observations)Filtered probabilities use the forward algorithm, updated with each new observation. Smoothed probabilities use Posterior function, considering all 100 observations for a global estimate. The Viterbi algorithm finds the most likely state sequence. Log-scale probabilities are exponentiated and normalized using exp and prop.table.
Task 4
Compute the accuracy of the filtered and smoothed probability distributions, and of the most probable path. That is, compute the percentage of the true hidden states that are guessed by each method.
Hint: Note that the function forward in the HMM package returns probabilities in log scale. You may need to use the functions exp and prop.table in order to obtain a normalized probability distribution. You may also want to use the functions apply and which.max to find out the most probable states. Finally, recall that you can compare two vectors A and B elementwise as A==B, and that the function table will count the number of times that the different elements in a vector occur in the vector.
Solution
# (4) Compute accuracy
filtered_estimates =apply(filtered_probs, 2, which.max)
smoothed_estimates =apply(smoothed_probs, 2, which.max)
true_states_num =as.numeric(true_states)
filtered_accuracy =mean(filtered_estimates == true_states_num)
smoothed_accuracy =mean(smoothed_estimates == true_states_num)
viterbi_accuracy =mean(viterbi_path == true_states)
cat("Accuracy measures:\n")## Accuracy measures:cat("Filtered:", filtered_accuracy, "%\n")## Filtered: 0.53 %cat("Smoothed:", smoothed_accuracy, "%\n")## Smoothed: 0.74 %cat("Viterbi:", viterbi_accuracy, "%\n")## Viterbi: 0.56 %# Optional: View probability distributions for first few time points
print("Filtered probabilities for first 5 time points:")## [1] "Filtered probabilities for first 5 time points:"print(filtered_probs[, 1:5])## index
## states 1 2 3 4 5
## 1 0.0 0.0 0.0000000 0.1666667 0.3846154
## 2 0.0 0.0 0.0000000 0.0000000 0.1153846
## 3 0.0 0.0 0.0000000 0.0000000 0.0000000
## 4 0.0 0.0 0.0000000 0.0000000 0.0000000
## 5 0.2 0.0 0.0000000 0.0000000 0.0000000
## 6 0.2 0.0 0.0000000 0.0000000 0.0000000
## 7 0.2 0.0 0.0000000 0.0000000 0.0000000
## 8 0.2 0.4 0.2222222 0.1111111 0.0000000
## 9 0.2 0.4 0.4444444 0.3333333 0.0000000
## 10 0.0 0.2 0.3333333 0.3888889 0.5000000print("Smoothed probabilities for first 5 time points:")## [1] "Smoothed probabilities for first 5 time points:"print(smoothed_probs[, 1:5])## index
## states 1 2 3 4 5
## 1 0.0000000 0.0000000 0.00000000 0.2003536 0.50540765
## 2 0.0000000 0.0000000 0.00000000 0.0000000 0.04873128
## 3 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
## 4 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
## 5 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
## 6 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
## 7 0.1534318 0.0000000 0.00000000 0.0000000 0.00000000
## 8 0.4241889 0.3068636 0.06859401 0.0000000 0.00000000
## 9 0.4223794 0.5415141 0.47653910 0.2057820 0.00000000
## 10 0.0000000 0.1516223 0.45486689 0.5938644 0.44586107Task 5
Repeat the previous exercise with different simulated samples. In general, the smoothed distributions should be more accurate than the filtered distributions. Why ? In general, the smoothed distributions should be more accurate than the most probable paths, too. Why?
Solution
# (5) Repeat with different samples
n_simulations =10
accuracies =matrix(0, nrow = n_simulations, ncol = 3)
colnames(accuracies) =c("Filtered", "Smoothed", "Viterbi")
for (i in 1:n_simulations) {
sim =simHMM(hmm, 100)
obs =sim$observation
true =sim$states
# Filtered
flog =forward(hmm, obs)
fprobs =apply(flog, 2, function(x) prop.table(exp(x - max(x))))
fest =apply(fprobs, 2, which.max)
# Smoothed
sprobs =posterior(hmm, obs)
sest =apply(sprobs, 2, which.max)
# Viterbi
vpath =viterbi(hmm, obs)
# Store accuracies
true_num =as.numeric(true)
accuracies[i, ] =c(
mean(fest == true_num),
mean(sest == true_num),
mean(vpath == true)
)
}
cat("Results and Explanation")## Results and Explanationcat("\nAverage accuracies over", n_simulations, "simulations:\n")##
## Average accuracies over 10 simulations:print(colMeans(accuracies))## Filtered Smoothed Viterbi
## 0.513 0.695 0.536Understanding Smoothed Distributions vs. Filtered Distributions and Viterbi
Why Smoothed Distributions are More Accurate than Filtered Distributions
Filtering only uses information up to the current moment, while smoothing also looks ahead at what happens later.This gives a clearer and more accurate picture.
By including future information, smoothing helps reduce uncertainty that comes from only looking at the past.
Why Smoothed Distributions are More Accurate than the Most Probable Path (Viterbi)
The Viterbi method gives you the single most likely sequence, but that doesn’t always mean each step in that path is the best estimate.Smoothing instead shows the likelihood of each state at every step.
In Viterbi, if one step is wrong, the error can spread through the
rest of the path.
Smoothing looks at each step separately with all the information, so
mistakes don’t pile up.
Task 6: Analyzing the Evolution of Uncertainty
Is it always true that the later in time (i.e., the more observations you have received) the better you know where the robot is ? Hint: You may want to compute the entropy of the filtered distributions with the function entropy.empirical of the package entropy.
Solution
# (6) Analyze entropy of filtered distributions
entropies =apply(filtered_probs, 2, function(p) entropy.empirical(p))
plot(entropies, type = "l", main = "Entropy of Filtered Distributions Over Time",
xlab = "Time", ylab = "Entropy")
Entropy plot: Shows fluctuations, not a clear decrease over time .
Observation: Entropy varies but does not consistently decrease.
# Check if later time steps have lower entropy on average
early_entropy =mean(entropies[1:50])
late_entropy =mean(entropies[51:100])
cat("Average entropy in first 50 steps:", early_entropy, "\n")## Average entropy in first 50 steps: 1.000262cat("Average entropy in last 50 steps:", late_entropy, "\n")## Average entropy in last 50 steps: 0.8717555Entropy measures uncertainty; lower entropy indicates better knowledge. It’s not always true that later observations improve location knowledge due to the noisy tracking (±2 sectors), causing persistent ambiguity regardless of time, especially with the circular structure.
Task 7: Predicting Time Step 101
Consider any of the samples above of length 100. Compute the probabilities of the hidden states for the time step 101.
Solution
# (7) Predict state at time step 101
# Get filtered probabilities for last time point
last_filtered =filtered_probs[, 100]
# Predict state probabilities for time 101
predicted_probs_101 =last_filtered %*% transProbs
predicted_state_101 =which.max(predicted_probs_101)
cat("\nPredicted state probabilities for time step 101:\n")##
## Predicted state probabilities for time step 101:print(predicted_probs_101)## 1 2 3 4 5 6 7 8 9 10
## [1,] 0 0.1875 0.5 0.3125 0 0 0 0 0 0cat("Most probable state at time 101:", predicted_state_101, "\n")## Most probable state at time 101: 3