更新了简单代码实现
高斯分布、多元高斯分布（Q：混合和多元的区别）
为什么要使用混合高斯分布
混合高斯分布概率密度函数定义
实质全概率公式：自然界中抽一个瓜正好抽中样本瓜x的概率
该分布由k个混合成分（高斯分布）组成，是混合系数、是第i个高斯分布的参数。
其中初始值为各成分的先验概率。
EM（Expectation Maximization）思想和算法
公式推导（贝叶斯、拉格朗日乘数法、最大似然函数）
假设
求随机选择一个瓜，其属于第z类的概率是多少
已知m个瓜的特征，求瓜分类的规律
最大似然（对数）函数
怎样迭代：
——通过样本加权平均值估计
求，有
因为，代入得
得到结果
——通过样本加权平均值估计
求，得
其左边展开，简写为，得
其右边展开，简写为，得
两式代回，得
得到结果
——由样本属于该成分的平均后验概率确定
由该函数对求偏导等于0，得
两边同时乘以,得
两边同时对k个混合成分求和，得
适用范围、计算代价
Q1：为什么似然函数无法得到解析解
推导得到各参数的最终公式后，我们就可以很简单的来构建一个GMM模型的代码
不会的东西就去Sklearn的GMM源代码里去学嘛，再不行就Github嘛
下面我们先定义需要的函数
# X: 样本特征矩阵
# Mu: 高斯分布均值矩阵
# Var: 高斯分布协方差矩阵
# W：后验概率矩阵
# Pi：簇比重
def log_prob(X, mu, var):
norm = multivariate_normal(mean=mu, cov=var)
return norm.logpdf(X)
def logLHSE(X, Pi, Mu, Var):
n_points, n_clusters = len(X), len(Pi)
pdfs = np.zeros(((n_points, n_clusters)))
for i in range(n_clusters):
pdfs[:, i] = np.mat(np.log(Pi[i]) + log_prob(X, Mu[i], np.diag(Var[i]))).reshape([len(X)])
return np.sum(pdfs)
def update_W(X, Mu, Var, Pi):
#W[N,K]
n_points, n_clusters = len(X), len(Pi)
pdfs = np.zeros(((n_points, n_clusters)))
for i in range(n_clusters):
pdfs[:, i] = Pi[i] * multivariate_normal.pdf(X, Mu[i], np.diag(Var[i]))
W = pdfs / pdfs.sum(axis=1).reshape(-1, 1)
return W
def update_Pi(W):
#Pi[K]
Pi = W.sum(axis=0) / W.sum()
return Pi
def update_Mu(X, W):
#Mu[K,D]
n_clusters = W.shape[1]
n_components = X.shape[1]
Mu = np.zeros((n_clusters, n_components))
for i in range(n_clusters):
Mu[i,:] = np.dot(W[:,i].T,X)/W[:,i].sum()
return Mu
#计算var，防止非正定矩阵reg_covar
def update_Var(X, Mu, W):
#Var[K,D,D]
n_clusters = Mu.shape[1]
reg_covar = 1e-7
Var = np.zeros((Mu.shape[0],n_clusters, n_clusters))
for i in range(Mu.shape[0]):
diff = np.mat(X - Mu[i])
Var[i] = diff.T * np.multiply(diff,W[:,i][:,np.newaxis]) / W[:,i].sum()
Var.flat[::len(X[0])+1] += reg_covar
return Var
def plot_clusters(X,y):
plt.figure(figsize=(15, 8))
sns.scatterplot(x=X.values[:, 0],y=X.values[:, 1],hue=y,
palette=sns.color_palette('hls', len(np.unique(y))), s=50)
fontdict = dict(fontsize=15)
plt.title('Cluster of credit card customer', fontdict, pad=10)
plt.show()
这样就搞定啦，是不是很简单呢。接下来我们用kaggle数据集简单示范下
别问为什么不包装，我懒
if __name__ == '__main__':
#导入你需要的数据
raw_df = pd.read_csv('./CC_GENERAL.csv')
raw_df = raw_df.drop('CUST_ID', axis = 1)
raw_df.fillna(method = 'ffill', inplace = True)
#数据预处理，PCA降维为可视化做准备
n_components = 2
# 标准化
scaler = StandardScaler()
scaled_df = scaler.fit_transform(raw_df)
# 归一化
normalized_df = normalize(scaled_df)
# 转成Dataframe
X = pd.DataFrame(normalized_df)
# PCA降维
pca = PCA(n_components)
X_pca = pca.fit_transform(X)
X_pca = pd.DataFrame(X_pca)
X_pca.columns = ['P1', 'P2']
# Kmeans初始均值+均分簇比重
n_clusters = 7
n_points = len(X)
Mu = cluster.KMeans(n_clusters=n_clusters).fit(X).cluster_centers_
Var = [np.cov(X.T) + 1e-3 * np.eye(X.shape[1])]*n_clusters
W = np.ones((n_points, n_clusters)) / n_clusters
Pi = W.sum(axis=0) / W.sum()
logfuc = []
for i in range(21):
if(i%4==0):
y = np.argmax(W,axis=1)
plot_clusters(X_pca,y)
logfuc.append(logLHSE(X, Pi, Mu, Var))
print('log-likelihood:%.3f'%logfuc[-1])
W = update_W(X, Mu, Var, Pi)
Pi = update_Pi(W)
Mu = update_Mu(X, W)
Var = update_Var(X, Mu, W)
metrics.silhouette_score(X,y, metric='cosine')
事实上，效果很一般，这种高维还重叠严重的数据感觉不适合聚类
都说轮廓系数越接近1越好，DBI指数越低越好，BIC指数越低越好
实际测试发现这里随簇类数量增加，轮廓系数降低，DBI增加,BIC下降