How to infer MAP estimate of HMM with PyStruct
Maximum a posteriori estimation of label series by HMM using PyStruct.
Preparation
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from pystruct.inference import inference_dispatch
The data is as follows. Try integerization and smoothing of time series data.
data
d = '17.2 19.7 21.6 21.3 22.1 20.5 16.3 18.4 21.0 16.1 17.5 18.5 18.4 18.3 16.0 21.2 18.8 24.3 23.3 20.5 16.9 22.4 20.1 24.5 24.2 22.7 19.6 23.6 23.3 24.6 25.0 24.3 22.2 22.7 19.5 20.5 17.3 17.2 22.0 20.9 21.5 22.3 24.0 22.4 20.2 15.7 20.4 16.3 17.7 14.3 18.4 16.6 13.9 15.2 14.8 15.0 11.5 13.4 13.5 17.0 15.0 17.5 12.3 11.8 14.5 12.4 12.9 15.8 13.8 11.4 6.5 5.9 7.2 5.6 4.6 7.5 8.9 6.6 3.9 5.7 7.3 6.1 6.8 3.1 2.6 7.9 5.2 2.0 4.0 3.4 5.7 8.1 4.7 5.4 5.9 3.6 2.9 5.7 2.1 1.6 2.3 2.4 1.2 4.2 4.2 2.4 5.6 2.5 3.0 6.1 4.9 7.1 5.0 7.2 5.2 5.1 10.4 8.3 6.9 6.8 7.8 4.2 8.0 3.2 7.9 5.9 9.5 6.4 9.2 11.7 11.6 15.5 16.7'
d = np.array([ float(c) for c in d.split()])
Now, build the HMM model and execute inference. --The unary term is the ndarray of (number of nodes, number of classes) --The pairwise term is the ndarray of (number of classes, number of classes). If you want to change it for each edge, ndarray of (number of edges, number of classes, number of classes) (There is a place where the API reference manual is wrong. Note.).
Run!
nClasses = 20 #Number of discrete classes. 20
p = 10 #Cost when the classes of adjacent nodes are different. (0 if same)
#Data term. Absolute value of difference from label
unaries = np.array([ [ abs(i-j) for j in range(nClasses) ] for i in d ])
#Pairwise term. 0 for the same label, p for different
pairwise = np.array((np.ma.ones((nClasses, nClasses)) - np.eye(nClasses)) * p)
#Setting the edge that connects the nodes. Since it is an HMM, it is one-dimensional.
edges = np.array([ [i, i+1] for i in range(unaries.shape[0]-1) ])
#Inference execution. Since this function maximizes, we add a minus to unaryies and pairwise.
res = inference_dispatch(-unaries, -pairwise, edges, inference_method='max-product')
#plot.
plt.plot(d, label="data")
plt.plot(res, label="result")
plt.legend()
Recommended Posts