机器学习之Kmeans算法推导(EM算法)以及python实现

Kmeans算法

Kmeans将样本集合划分为$k$个类，将$n$个样本分到$k$个类中，每个样本到其所属类的中心距离最小。每个样本只能属于一个类，所以Kmeans算法是硬聚类。

参考资料李航博士的《统计学习方法》
这篇文章增加了EM算法推导以及Python实现

1. Kmeans算法的推导

模型：
给定$n$个样本集合$X=x_1,x_2,...,x_n$，每个样本的特征维度是$m$。将样本分为$k$个类$C_1,C_2,...,C_k$,
每个样本只能属于一个类。
首先采用的是欧式距离平方作为样本之间的距离，定义K means算法的损失函数为：

$$W(C)=\sum _{l=1}^k\sum _{C(i)=l}||x_i-{\overline{x}}_l|{|}^2$$

其中${\overline{x}}_l$代表第$l$个类的中心坐标。

下面运用EM算法的思路推导K means：
1.明确隐变量：
在K means算法中反映观测数据$x_i$来自与哪一个类是未知的，设隐变量${\gamma }_{il}$表示：

$${\gamma }_{il}=\{\begin{array}{ll}1,& if第i个观测值来自于第l个类\\ o,& else\end{array}$$

$$i=1,2,...,n;l=1,2,...,k$$

2.E步：
设完全参数的损失函数为：

$$W({\overline{x}}_l)=\sum _{l=1}^k\sum _{i=1}^n{\gamma }_{il}||x_i-{\overline{x}}_l|{|}^2$$

有隐变量的定义可得：

$${\hat{\gamma }}_{il}=[(\arg \underset{1\le l\le k}{\min }||x_i-{\overline{x}}_l|{|}^2)==l]$$

则：

$$Q({\overline{x}}_l)=\sum _{l=1}^k\sum _{i=1}^n{\hat{\gamma }}_{il}||x_i-{\overline{x}}_l|{|}^2$$

3.M步：
求解参数${\overline{x}}_l$
对$Q({\overline{x}}_l)$求导，并令其导数为0

$$\frac{\partial }{\partial {\overline{x}}_l}(\sum _{l=1}^k\sum _{i=1}^n{\hat{\gamma }}_{il}||x_i-{\overline{x}}_l|{|}^2)=0$$

得到：

$$2\sum _{i=1}^n{\hat{\gamma }}_{il}(x_i-{\overline{x}}_l)=0$$

故：

$${\overline{x}}_l=\frac{{\sum }_{i=1}^n{\hat{\gamma }}_{il}{x}_i}{{\sum }_{i=1}^n{\hat{\gamma }}_{il}}$$

Kmeans算法的总结：
E步：

$${\hat{\gamma }}_{il}=[(\arg \underset{1\le l\le k}{\min }||x_i-{\overline{x}}_l|{|}^2)==l]$$

M步：

$${\overline{x}}_l=\frac{{\sum }_{i=1}^n{\hat{\gamma }}_{il}{x}_i}{{\sum }_{i=1}^n{\hat{\gamma }}_{il}}$$

仔细一看这两个公式就会发现E步就是在求第$i$个样本是否属于第$l$个类，而M步则在求第$l$个类的中心点。

2. Kmeans算法的python实现

2.1 Python实现模型

import numpy as np
np.random.seed(0)


class MyKmeans(object):
    def __init__(self, k=3):
        """
        运用EM算法解决Kmeans问题
        :param k: 超参数，表示分类的数量

        涉及到的其它参数
        :param n: 表示样本的个数
        :param m: 表示单个样本的特征维度
        :param gamma: 此模型的隐变量，表示样本属于哪一个类， 维度(n, k)
        :param x_l: 表示类的中心坐标，维度(k, m)
        """
        self.k = k
        self.params = {
            'gamma': None,
            'x_l': None
        }

        self.n = None
        self.m = None

    def _init_params(self, x):
        """
        初始化参数
        :param x: 输入样本
        :return:
        """
        # 由于在第一次E步中会迭代gamma参数，所以gamma可以随意赋值
        gamma = np.zeros((self.n, self.k))
        # 避免初始化的类离观测数据太远，随机选择观测数据的k个点为类中心
        x_l = np.array([x[i] for i in np.random.randint(0, self.n, self.k)])
        self.params = {'gamma': gamma, 'x_l': x_l}

    def _E_step(self, x):
        """
        迭代隐变量
        :param x: 输入样本(n, m)
        :return:
        """
        x_l = self.params['x_l']
        for i in range(self.n):
            x_i = x[i]
            gamma_list = []
            for l in range(self.k):
                x_l_l = x_l[l]
                # 计算第i个样本和第l类中心的距离
                gamma_list.append(np.dot((x_i - x_l_l), (x_i - x_l_l).reshape(self.m, -1)))
            ret_list = np.zeros(self.k)
            index = np.argmin(gamma_list)
            ret_list[index] = 1
            self.params['gamma'][i] = ret_list

    def _M_step(self, x):
        """
        M步迭代计算类中心
        :param x:
        :return:
        """
        gamma = self.params['gamma']
        # 使用循环依次计算元素
        for l in range(self.k):
            gamma_l = gamma[:, l]
            x_l_list = []
            for i in range(self.n):
                if gamma_l[i]:
                    x_l_list.append(x[i])
            # 计算第l类的类中心坐标
            self.params['x_l'][l] = np.average(x_l_list, axis=0)

        # # 也可以使用矩阵直接计算
        # sum = np.dot(x.T, gamma).T
        # self.params['x_l'] = np.array([i / j for i, j in zip(sum, np.sum(gamma, axis=0))])

    def fit(self, x, max_iter=30):
        """
        模型入口
        :param x: 输入的样本（n, m）
        :param max_iter: 最大循环次数
        :return:
        """
        x = np.array(x)
        self.n, self.m = x.shape

        self._init_params(x)

        for i in range(max_iter):
            self._E_step(x)
            self._M_step(x)

2.2 测试模型

def get_samples(n_ex=100, n_classes=3, n_in=2, seed=0):
    # 生成100个样本，为了能够在二维平面上画出图线表示出来，每个样本的特征维度设置为2
    from sklearn.datasets.samples_generator import make_blobs
    from sklearn.model_selection import train_test_split
    y, _ = make_blobs(
        n_samples=n_ex, centers=n_classes, n_features=n_in, random_state=seed)
    return y


def run_my_model():
    from matplotlib import pyplot as plt
    my = MyKmeans()
    y = get_samples()
    my.fit(y)

    max_index = np.argmax(my.params['gamma'], axis=1)
    print(max_index)

    k1_list = []
    k2_list = []
    k3_list = []

    for y_i, index in zip(y, max_index):
        if index == 0:
            k1_list.append(y_i)
        elif index == 1:
            k2_list.append(y_i)
        else:
            k3_list.append(y_i)
    k1_list = np.array(k1_list)
    k2_list = np.array(k2_list)
    k3_list = np.array(k3_list)

    plt.scatter(k1_list[:, 0], k1_list[:, 1], c='red')
    plt.scatter(k2_list[:, 0], k2_list[:, 1], c='blue')
    plt.scatter(k3_list[:, 0], k3_list[:, 1], c='green')
    plt.show()

结果如图所示：