EM、k-means、GMM代码样例

EM算法样例代码：

import numpy as np

# 定义高斯分布函数
def gaussian(x, mean, cov):
    n = x.shape[0]
    exp_part = np.exp(-0.5 * (x - mean).T.dot(np.linalg.inv(cov)).dot(x - mean))
    coef = 1 / np.sqrt(((2 * np.pi) ** n) * np.linalg.det(cov))
    return coef * exp_part

# EM算法
def EM(X, K, max_iter):
    n, m = X.shape
    # 初始化参数
    pi = np.ones(K) / K
    mu = X[np.random.choice(n, K, replace=False)]
    sigma = [np.eye(m) for i in range(K)]
    # 迭代
    for iter in range(max_iter):
        # E步
        gamma = np.zeros((n, K))
        for i in range(n):
            for j in range(K):
                gamma[i, j] = pi[j] * gaussian(X[i], mu[j], sigma[j])
            gamma[i] /= np.sum(gamma[i])
        # M步
        N_k = np.sum(gamma, axis=0)
        for j in range(K):
            mu[j] = np.sum(gamma[:, j].reshape(-1, 1) * X, axis=0) / N_k[j]
            sigma[j] = (X - mu[j]).T.dot(gamma[:, j].reshape(-1, 1) * (X - mu[j])) / N_k[j]
            pi[j] = N_k[j] / n
    return pi, mu, sigma

k-means算法样例代码：

import numpy as np

# k-means算法
def kmeans(X, K, max_iter):
    n, m = X.shape
    # 随机初始化聚类中心
    centers = X[np.random.choice(n, K, replace=False)]
    # 迭代
    for iter in range(max_iter):
        # 计算每个样本到各个聚类中心的距离
        dists = np.zeros((n, K))
        for j in range(K):
            dists[:, j] = np.sum((X - centers[j]) ** 2, axis=1)
        # 将样本划分到最近的聚类中心
        labels = np.argmin(dists, axis=1)
        # 更新聚类中心
        for j in range(K):
            if np.sum(labels == j) > 0:
                centers[j] = np.mean(X[labels == j], axis=0)
    return centers, labels

GMM算法样例代码：

import numpy as np

# 定义高斯分布函数
def gaussian(x, mean, cov):
    n = x.shape[0]
    exp_part = np.exp(-0.5 * (x - mean).T.dot(np.linalg.inv(cov)).dot(x - mean))
    coef = 1 / np.sqrt(((2 * np.pi) ** n) * np.linalg.det(cov))
    return coef * exp_part

# GMM算法
def GMM(X, K, max_iter):
    n, m = X.shape
    # 初始化参数
    pi = np.ones(K) / K
    mu = X[np.random.choice(n, K, replace=False)]
    sigma = [np.eye(m) for i in range(K)]
    # 迭代
    for iter in range(max_iter):
        # E步
        gamma = np.zeros((n, K))
        for i in range(n):
            for j in range(K):
                gamma[i, j] = pi[j] * gaussian(X[i], mu[j], sigma[j])
            gamma[i] /= np.sum(gamma[i])
        # M步
        N_k = np.sum(gamma, axis=0)
        for j in range(K):
            mu[j] = np.sum(gamma[:, j].reshape(-1, 1) * X, axis=0) / N_k[j]
            sigma[j] = np.zeros((m, m))
            for i in range(n):
                sigma[j] += gamma[i, j] * np.outer(X[i] - mu[j], X[i] - mu[j])
            sigma[j] /= N_k[j]
            pi[j] = N_k[j] / n
    return pi, mu, sigma