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