| 18 | |
| 19 | |
| 20 | class HMM: |
| 21 | def __init__(self, M): |
| 22 | self.M = M # number of hidden states |
| 23 | |
| 24 | def fit(self, X, max_iter=30): |
| 25 | np.random.seed(123) |
| 26 | # train the HMM model using the Baum-Welch algorithm |
| 27 | # a specific instance of the expectation-maximization algorithm |
| 28 | |
| 29 | # determine V, the vocabulary size |
| 30 | # assume observables are already integers from 0..V-1 |
| 31 | # X is a jagged array of observed sequences |
| 32 | V = max(max(x) for x in X) + 1 |
| 33 | N = len(X) |
| 34 | |
| 35 | self.pi = np.ones(self.M) / self.M # initial state distribution |
| 36 | self.A = random_normalized(self.M, self.M) # state transition matrix |
| 37 | self.B = random_normalized(self.M, V) # output distribution |
| 38 | |
| 39 | print("initial A:", self.A) |
| 40 | print("initial B:", self.B) |
| 41 | |
| 42 | costs = [] |
| 43 | for it in range(max_iter): |
| 44 | if it % 10 == 0: |
| 45 | print("it:", it) |
| 46 | # alpha1 = np.zeros((N, self.M)) |
| 47 | alphas = [] |
| 48 | betas = [] |
| 49 | scales = [] |
| 50 | logP = np.zeros(N) |
| 51 | for n in range(N): |
| 52 | x = X[n] |
| 53 | T = len(x) |
| 54 | scale = np.zeros(T) |
| 55 | # alpha1[n] = self.pi*self.B[:,x[0]] |
| 56 | alpha = np.zeros((T, self.M)) |
| 57 | alpha[0] = self.pi*self.B[:,x[0]] |
| 58 | scale[0] = alpha[0].sum() |
| 59 | alpha[0] /= scale[0] |
| 60 | for t in range(1, T): |
| 61 | alpha_t_prime = alpha[t-1].dot(self.A) * self.B[:, x[t]] |
| 62 | scale[t] = alpha_t_prime.sum() |
| 63 | alpha[t] = alpha_t_prime / scale[t] |
| 64 | logP[n] = np.log(scale).sum() |
| 65 | alphas.append(alpha) |
| 66 | scales.append(scale) |
| 67 | |
| 68 | beta = np.zeros((T, self.M)) |
| 69 | beta[-1] = 1 |
| 70 | for t in range(T - 2, -1, -1): |
| 71 | beta[t] = self.A.dot(self.B[:, x[t+1]] * beta[t+1]) / scale[t+1] |
| 72 | betas.append(beta) |
| 73 | |
| 74 | |
| 75 | cost = np.sum(logP) |
| 76 | costs.append(cost) |
| 77 | |