markov-models/notebook.py
2021-01-05 11:10:29 +00:00

433 lines
11 KiB
Python

# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
#IMPORTS AND COMMON VARIABLES
import matplotlib
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from matplotlib.pyplot import savefig
import numpy as np
from math import sqrt
from constants import *
from maths import gaussian
from markov import MarkovModel
from markovlog import LogMarkovModel
fig_dpi = 200
fig_export = False
x = np.linspace(-4, 8, 300) # x values for figures
x_label = "Observation Space"
y_label = "Probability Density"
# %% [markdown]
# State Probability Functions (1)
# ===================
# %%
state_1_y = [gaussian(i, state1.mean, state1.std_dev) for i in x]
state_2_y = [gaussian(i, state2.mean, state2.std_dev) for i in x]
plt.plot(x, state_1_y, c='r', label="State 1")
plt.plot(x, state_2_y, c='b', label="State 2")
plt.legend()
plt.title("State Probability Density Functions")
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.grid(linestyle="--")
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
fig.set_tight_layout(True)
if fig_export:
savefig("report/res/pdfs.png")
plt.show()
# %% [markdown]
# Output Probability Densities (2)
# ==========
# %%
for obs in observations:
print(f'{obs} -> State 1: {gaussian(obs, state1.mean, state1.std_dev)},',
f'State 2: {gaussian(obs, state2.mean, state2.std_dev)}')
# %%
state_1_y = [gaussian(i, state1.mean, state1.std_dev) for i in x]
state_2_y = [gaussian(i, state2.mean, state2.std_dev) for i in x]
plt.plot(x, state_1_y, c='r', label="State 1")
plt.plot(x, state_2_y, c='b', label="State 2")
plt.legend()
plt.title("State Probability Density Functions With Observations")
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.grid(linestyle="--", axis='y')
state1_pd = [gaussian(i, state1.mean, state1.std_dev) for i in observations]
state2_pd = [gaussian(i, state2.mean, state2.std_dev) for i in observations]
#############################################
# Observation Marks
#############################################
config = {
"s": 65,
"marker": 'x'
}
[plt.axvline(x=i, ls='--', lw=1.0, c=(0,0,0), alpha=0.4) for i in observations]
plt.scatter(observations, state1_pd, color=(0.5, 0, 0), **config)
plt.scatter(observations, state2_pd, color=(0, 0, 0.5), **config)
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
fig.set_tight_layout(True)
if fig_export:
savefig("report/res/pdfs-w-obs.png")
plt.show()
# %% [markdown]
# # Forward Procedure (3)
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition)
model.populate_forward()
print(model.forward)
forward = model.forward
model.calculate_p_obs_forward()
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
state_x = np.arange(1, 10)
from numpy import log as ln
plt.plot(state_x, [ln(i) for i in model.forward[0, :]], c='r', label="State 1")
plt.plot(state_x, [ln(i) for i in model.forward[1, :]], c='b', label="State 2")
plt.ylim(top=0)
plt.legend()
plt.title("Forward Log-Likelihoods Over Time")
plt.xlabel("Observation (t)")
plt.ylabel("Log Likelihood")
plt.grid(linestyle="--")
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
fig.set_tight_layout(True)
if fig_export:
savefig("report/res/forward-logline.png")
plt.show()
# %% [markdown]
# # Backward Procedure (4)
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition)
model.populate_backward()
print(model.backward)
backward = model.backward
model.calculate_p_obs_backward()
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
state_x = np.arange(1, 10)
from numpy import log as ln
plt.plot(state_x, [ln(i) for i in model.backward[0, :]], c='r', label="State 1")
plt.plot(state_x, [ln(i) for i in model.backward[1, :]], c='b', label="State 2")
plt.ylim(top=0)
plt.legend()
plt.title("Backward Log-Likelihoods Over Time")
plt.xlabel("Observation (t)")
plt.ylabel("Log Likelihood")
plt.grid(linestyle="--")
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
fig.set_tight_layout(True)
if fig_export:
savefig("report/res/backward-logline.png")
plt.show()
# %% [markdown]
# # Compare Forward/Backward Final
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition)
model.populate_forward()
model.populate_backward()
print("forward:", model.calculate_p_obs_forward())
print("backward:", model.calculate_p_obs_backward())
print("diff: ", model.p_obs_forward - model.p_obs_backward)
# %% [markdown]
# # Occupation Likelihoods (5)
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
occupation = model.occupation
print(model.occupation)
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
fig = plt.figure(figsize=(6,6), dpi=fig_dpi, tight_layout=True)
ax = fig.add_subplot(1, 1, 1, projection="3d", xmargin=0, ymargin=0)
y_width = 0.3
X = np.arange(1, 10) - 0.5
Y = np.arange(1, 3) - 0.5*y_width
X, Y = np.meshgrid(X, Y)
Z = np.zeros(model.forward.size)
dx = np.ones(model.forward.size)
dy = y_width * np.ones(model.forward.size)
colours = [*[(1.0, 0.1, 0.1) for i in range(9)], *[(0.2, 0.2, 1.0) for i in range(9)]]
ax.bar3d(X.flatten(), Y.flatten(), Z,
dx, dy, model.occupation.flatten(),
color=colours, shade=True)
ax.set_yticks([1, 2])
ax.set_zlim(top=1.0)
ax.set_title("Occupation Likelihoods Over Time")
ax.set_xlabel("Observation")
ax.set_ylabel("State")
ax.set_zlabel("Occupation Likelihood")
ax.view_init(35, -72)
if fig_export:
savefig("report/res/occupation-bars.png")
fig.show()
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
state_x = np.arange(1, 10)
plt.plot(state_x, model.occupation[0, :], c='r', label="State 1")
plt.plot(state_x, model.occupation[1, :], c='b', label="State 2")
plt.legend()
plt.title("Occupation Likelihoods Over Time")
plt.xlabel("Observation (t)")
plt.ylabel("Occupation Likelihood")
plt.grid(linestyle="--")
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
fig.set_tight_layout(True)
if fig_export:
savefig("report/res/occupation-line.png")
plt.show()
# %% [markdown]
# # Re-estimate Mean & Variance (6)
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
print("mean: ", [state1.mean, state2.mean])
print("variance: ", [state1.variance, state2.variance])
print()
print("mean: ", model.reestimated_mean())
print("variance: ", model.reestimated_variance())
# %% [markdown]
# New PDFs (7)
# ===================
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
new_mean = model.reestimated_mean()
new_var = model.reestimated_variance()
new_std_dev = [sqrt(x) for x in new_var]
state_1_y = [gaussian(i, new_mean[0], new_std_dev[0]) for i in x]
state_2_y = [gaussian(i, new_mean[1], new_std_dev[1]) for i in x]
plt.plot(x, state_1_y, c='r', label="State 1")
plt.plot(x, state_2_y, c='b', label="State 2")
plt.legend()
plt.title("Re-estimated Probability Density Functions")
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.grid(linestyle="--")
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
plt.show()
# %% [markdown]
# # Compare PDFs (7)
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
new_mean = model.reestimated_mean()
new_var = model.reestimated_variance()
new_std_dev = [sqrt(x) for x in new_var]
#######################################
# Original
#######################################
state_1_y = [gaussian(i, state1.mean, state1.std_dev) for i in x]
state_2_y = [gaussian(i, state2.mean, state2.std_dev) for i in x]
plt.plot(x, state_1_y, '--', c='r', label="State 1", linewidth=1.0)
plt.plot(x, state_2_y, '--', c='b', label="State 2", linewidth=1.0)
#######################################
# Re-Estimated
#######################################
state_1_new_y = [gaussian(i, new_mean[0], new_std_dev[0]) for i in x]
state_2_new_y = [gaussian(i, new_mean[1], new_std_dev[1]) for i in x]
plt.plot(x, state_1_new_y, c='r', label="New State 1")
plt.plot(x, state_2_new_y, c='b', label="New State 2")
plt.legend()
plt.title("Re-estimated Probability Density Functions")
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.grid(linestyle="--")
fig = matplotlib.pyplot.gcf()
fig.set_dpi(fig_dpi)
fig.set_tight_layout(True)
if fig_export:
savefig("report/res/re-est-pdfs.png")
plt.show()
# %% [markdown]
# # Multiple Iterations
# %%
iterations = 50
fig = plt.figure(dpi=fig_dpi, tight_layout=True)
ax = fig.add_subplot(1, 1, 1, xmargin=0, ymargin=0)
iter_mean = [state1.mean, state2.mean]
iter_var = [state1.variance, state2.variance]
iter_state_transitions = state_transition
ax.plot(x, [gaussian(i, iter_mean[0], sqrt(iter_var[0])) for i in x], '--', c='r', linewidth=1.0)
ax.plot(x, [gaussian(i, iter_mean[1], sqrt(iter_var[1])) for i in x], '--', c='b', linewidth=1.0)
label1=None
label2=None
for i in range(iterations):
iter_model = MarkovModel(states=[State(iter_mean[0], iter_var[0], state1.entry, state1.exit),
State(iter_mean[1], iter_var[1], state2.entry, state2.exit)],
observations=observations,
state_transitions=iter_state_transitions).populate()
# NEW PARAMETERS
iter_mean = iter_model.reestimated_mean()
iter_var = iter_model.reestimated_variance()
iter_state_transitions[1:3, 1:3] = iter_model.reestimated_state_transitions()
print(f"mean ({i}): ", iter_mean)
print(f"var ({i}): ", iter_var)
print(iter_model.reestimated_state_transitions())
print()
state_1_y = [gaussian(i, iter_mean[0], sqrt(iter_var[0])) for i in x]
state_2_y = [gaussian(i, iter_mean[1], sqrt(iter_var[1])) for i in x]
style = '--'
linewidth = 1.0
if i == iterations - 1:
style = '-'
linewidth = 2.0
label1='State 1'
label2='State 2'
ax.plot(x, state_1_y, style, c='r', label=label1, linewidth=linewidth)
ax.plot(x, state_2_y, style, c='b', label=label2, linewidth=linewidth)
ax.set_title("Probability Density Function Iterations")
ax.set_xlabel(x_label)
ax.set_ylabel(y_label)
ax.grid(linestyle="--")
ax.legend()
if fig_export:
savefig("report/res/iterated-pdfs.png")
fig.show()
# %% [markdown]
# # Baum-Welch State Transition Re-estimations
# %%
model = MarkovModel(states=[state1, state2],
observations=observations,
state_transitions=state_transition).populate()
print(a_matrix)
print(model.reestimated_state_transitions())
model.reestimated_state_transitions()
# %%