-
Notifications
You must be signed in to change notification settings - Fork 60
/
gmm.py
141 lines (120 loc) · 4.61 KB
/
gmm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
# -*- coding: utf-8 -*-
# ----------------------------------------------------
# Copyright (c) 2017, Wray Zheng. All Rights Reserved.
# Distributed under the BSD License.
# ----------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
DEBUG = True
######################################################
# 调试输出函数
# 由全局变量 DEBUG 控制输出
######################################################
def debug(*args, **kwargs):
global DEBUG
if DEBUG:
print(*args, **kwargs)
######################################################
# 第 k 个模型的高斯分布密度函数
# 每 i 行表示第 i 个样本在各模型中的出现概率
# 返回一维列表
######################################################
def phi(Y, mu_k, cov_k):
norm = multivariate_normal(mean=mu_k, cov=cov_k)
return norm.pdf(Y)
######################################################
# E 步:计算每个模型对样本的响应度
# Y 为样本矩阵,每个样本一行,只有一个特征时为列向量
# mu 为均值多维数组,每行表示一个样本各个特征的均值
# cov 为协方差矩阵的数组,alpha 为模型响应度数组
######################################################
def getExpectation(Y, mu, cov, alpha):
# 样本数
N = Y.shape[0]
# 模型数
K = alpha.shape[0]
# 为避免使用单个高斯模型或样本,导致返回结果的类型不一致
# 因此要求样本数和模型个数必须大于1
assert N > 1, "There must be more than one sample!"
assert K > 1, "There must be more than one gaussian model!"
# 响应度矩阵,行对应样本,列对应响应度
gamma = np.mat(np.zeros((N, K)))
# 计算各模型中所有样本出现的概率,行对应样本,列对应模型
prob = np.zeros((N, K))
for k in range(K):
prob[:, k] = phi(Y, mu[k], cov[k])
prob = np.mat(prob)
# 计算每个模型对每个样本的响应度
for k in range(K):
gamma[:, k] = alpha[k] * prob[:, k]
for i in range(N):
gamma[i, :] /= np.sum(gamma[i, :])
return gamma
######################################################
# M 步:迭代模型参数
# Y 为样本矩阵,gamma 为响应度矩阵
######################################################
def maximize(Y, gamma):
# 样本数和特征数
N, D = Y.shape
# 模型数
K = gamma.shape[1]
#初始化参数值
mu = np.zeros((K, D))
cov = []
alpha = np.zeros(K)
# 更新每个模型的参数
for k in range(K):
# 第 k 个模型对所有样本的响应度之和
Nk = np.sum(gamma[:, k])
# 更新 mu
# 对每个特征求均值
mu[k, :] = np.sum(np.multiply(Y, gamma[:, k]), axis=0) / Nk
# 更新 cov
cov_k = (Y - mu[k]).T * np.multiply((Y - mu[k]), gamma[:, k]) / Nk
cov.append(cov_k)
# 更新 alpha
alpha[k] = Nk / N
cov = np.array(cov)
return mu, cov, alpha
######################################################
# 数据预处理
# 将所有数据都缩放到 0 和 1 之间
######################################################
def scale_data(Y):
# 对每一维特征分别进行缩放
for i in range(Y.shape[1]):
max_ = Y[:, i].max()
min_ = Y[:, i].min()
Y[:, i] = (Y[:, i] - min_) / (max_ - min_)
debug("Data scaled.")
return Y
######################################################
# 初始化模型参数
# shape 是表示样本规模的二元组,(样本数, 特征数)
# K 表示模型个数
######################################################
def init_params(shape, K):
N, D = shape
mu = np.random.rand(K, D)
cov = np.array([np.eye(D)] * K)
alpha = np.array([1.0 / K] * K)
debug("Parameters initialized.")
debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n")
return mu, cov, alpha
######################################################
# 高斯混合模型 EM 算法
# 给定样本矩阵 Y,计算模型参数
# K 为模型个数
# times 为迭代次数
######################################################
def GMM_EM(Y, K, times):
Y = scale_data(Y)
mu, cov, alpha = init_params(Y.shape, K)
for i in range(times):
gamma = getExpectation(Y, mu, cov, alpha)
mu, cov, alpha = maximize(Y, gamma)
debug("{sep} Result {sep}".format(sep="-" * 20))
debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n")
return mu, cov, alpha