(self, X, max_iter=30, eps=1e0)
| 27 | self.K = K # number of Gaussians |
| 28 | |
| 29 | def fit(self, X, max_iter=30, eps=1e0): |
| 30 | # train the HMM model using the Baum-Welch algorithm |
| 31 | # a specific instance of the expectation-maximization algorithm |
| 32 | |
| 33 | N = len(X) |
| 34 | D = X[0].shape[1] # assume each x is organized (T, D) |
| 35 | |
| 36 | self.pi = np.ones(self.M) / self.M # initial state distribution |
| 37 | self.A = random_normalized(self.M, self.M) # state transition matrix |
| 38 | self.R = np.ones((self.M, self.K)) / self.K # mixture proportions |
| 39 | print("initial A:", self.A) |
| 40 | print("initial R:", self.R) |
| 41 | self.mu = np.zeros((self.M, self.K, D)) |
| 42 | for i in range(self.M): |
| 43 | for k in range(self.K): |
| 44 | random_idx = np.random.choice(N) |
| 45 | x = X[random_idx] |
| 46 | random_time_idx = np.random.choice(len(x)) |
| 47 | self.mu[i,k] = x[random_time_idx] |
| 48 | self.sigma = np.zeros((self.M, self.K, D, D)) |
| 49 | for j in range(self.M): |
| 50 | for k in range(self.K): |
| 51 | self.sigma[j,k] = np.eye(D) |
| 52 | |
| 53 | costs = [] |
| 54 | for it in range(max_iter): |
| 55 | if it % 1 == 0: |
| 56 | print("it:", it) |
| 57 | alphas = [] |
| 58 | betas = [] |
| 59 | gammas = [] |
| 60 | Bs = [] |
| 61 | # components = [] |
| 62 | P = np.zeros(N) |
| 63 | |
| 64 | for n in range(N): |
| 65 | x = X[n] |
| 66 | T = len(x) |
| 67 | |
| 68 | # calculate B so we can lookup when updating alpha and beta |
| 69 | B = np.zeros((self.M, T)) |
| 70 | component = np.zeros((self.M, self.K, T)) # we'll need these later |
| 71 | for j in range(self.M): |
| 72 | for t in range(T): |
| 73 | for k in range(self.K): |
| 74 | p = self.R[j,k] * mvn.pdf(x[t], self.mu[j,k], self.sigma[j,k]) |
| 75 | component[j,k,t] = p |
| 76 | B[j,t] += p |
| 77 | Bs.append(B) |
| 78 | |
| 79 | alpha = np.zeros((T, self.M)) |
| 80 | alpha[0] = self.pi*B[:,0] |
| 81 | for t in range(1, T): |
| 82 | alpha[t] = alpha[t-1].dot(self.A) * B[:,t] |
| 83 | P[n] = alpha[-1].sum() |
| 84 | assert(P[n] <= 1) |
| 85 | alphas.append(alpha) |
| 86 |
no test coverage detected