# -*- coding: UTF-8 -*-
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.datasets.samples_generator import make_blobs
from matplotlib.patches import Ellipse # 画椭圆的类
sns.set()
def draw_ellipse(position, covariance, ax=None, **kwargs):
# 做椭圆需要协方差和位置信息
ax = ax or plt.gca()
# 协方差换成主轴
if covariance.shape == (2, 2):
# 分解奇异值
print('covariance:\n', covariance)
# https://blog.csdn.net/rainpasttime/article/details/79831533?utm_source=app
U, s, Vt = np.linalg.svd(covariance) # compute_uv的取值是为0或者1，默认值为1，表示计算u,s,v。为0的时候只计算s
print('s:\n', s)
print('Vt:\n', Vt)
# arctan2(y, x) 是该点(x, y) 对应原点形成的三角形,求反正切对应的弧度数(过原点和x轴的那个角)
# 控制椭圆罩住的角度
angle = np.degrees(np.arctan2(U[1, 0], U[0, 0])) # 弧度换为角度，atan2 的用法一样（用C库所以更快）
print('angle is：', angle)
print('U:\n', U)
print('U[1, 0], U[0, 0]):\n', U[1, 0], U[0, 0])
width, height = 2 * np.sqrt(s) # 奇异值分配出去，椭圆大小
else: # 若不是2*2矩阵
angle = 0
width, height = 2 * np.sqrt(covariance) # 就用自己
# 画椭圆
for nsig in range(1, 4): # 遍历椭圆内概率密度分层的数量
# 绘制椭圆，第二个参数是椭圆的长直径，第三个是短直径，angle是逆时针旋转角度
ax.add_patch(Ellipse(position, nsig * width, nsig * height, angle, **kwargs))
def plot_gmm(gmm, X, label=True, ax=None):
ax = ax or plt.gca()
labels = gmm.fit(X).predict(X)
if label:
ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
else:
ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)
ax.axis('equal')
w_factor = 0.2 / gmm.weights_.max()
print('\ngmm.means_:\n', gmm.means_, '\ngmm.covariances_:\n', gmm.covariances_, '\ngmm.weights_:\n', gmm.weights_)
for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_): # 画出椭圆
print('w * w_factor:', w * w_factor)
draw_ellipse(pos, covar, alpha=w * w_factor) # alpha通过**kwargs传入
X, y_true = make_blobs(n_samples=400, centers=4, cluster_std=0.60, random_state=0)
X = X[:, ::-1]
rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2)) # 叉乘投影 (400,2) X (2,2)，变成椭圆
# gmm = GaussianMixture(n_components=4).fit(X)
# gmm = GaussianMixture(n_components=4, covariance_type='full', random_state=42).fit(X)
gmm = GaussianMixture(n_components=4, random_state=42).fit(X_stretched)
labels = gmm.predict(X)
probs = gmm.predict_proba(X)
# print(probs[:5].round(3))
size = 50 * probs.max(1) ** 2
# plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', s=size)
plot_gmm(gmm, X_stretched)
plt.show()
KAGAYA's Blog 