实验目的

熟悉和掌握k-means聚类算法
熟悉和掌握DBSCAN聚类算法

实验要求

采用Python、Matlab等高级语言进行编程，推荐优先选用Python语言
核心模型和算法需自主编程实现，不得直接调用Scikit-learn、PyTorch等成熟框架的第三方实现
代码可读性强：变量、函数、类等命名可读性强，包含必要的注释

实验原理

K-means算法是最简单的一种聚类算法。算法的目的是使各个样本与所在类均值的误差平方和达到最小（这也是评价K-means算法最后聚类效果的评价标准）
K-means聚类算法的一般步骤：
1. 初始化。输入基因表达矩阵作为对象集X，输入指定聚类类数N，并在X中随机选取N个对象作为初始聚类中心。设定迭代中止条件，比如最大循环次数或者聚类中心收敛误差容限。
2. 进行迭代。根据相似度准则将数据对象分配到最接近的聚类中心，从而形成一类。初始化隶属度矩阵。
3. 更新聚类中心。然后以每一类的平均向量作为新的聚类中心，重新分配数据对象。
4. 反复执行第二步和第三步直至满足中止条件。

举一个简单的例子来说明问题：

设有一组数据集x1=(2,1),x2=(1,3),x3=(6,7),x4=(4,7)

(1)选取聚类中心，该中心可以任意选取，也可以通过直方图进行选取，还可以通过取前2个值进行选取。我们选择两个聚类中心。

(2)计算每一个样本值到聚类中心的距离；并划分新的聚类中心；

评价标准：
$$
J(c, \mu)=\sum_{i=1}^m\left|x^{(i)}-\mu_{c^{(i)}}\right|^2
$$
假设有 $\mathrm{M}$ 个数据源， $\mathrm{C}$ 个聚类中心。 $\mu \mathrm{c}$ 为聚类中心。该公式的意思也就是将每个类中的数据与每个聚类中心做差的平方和, $\mathrm{J}$ 最小, 意味着分割的效果最好。

实验内容

K-means聚类算法

采用数据集 “data/clustering1.csv”进行k-means聚类算法实验。
数据集简介：本数据集来自UCI数据库Sales_Transactions_Dataset_Weekly数据，经过简单筛选后维度为666\times52，对应666种商品在52周的销量。

将每一种商品看作一个样本, 每一周看作一个属性, 对其进行聚类。
对数据进行 min-max 归一化至 $[0,1]$ 区间内:
随机挑选 $\mathrm{k}$ 个初始聚类中心, 使用欧氏距离度量各样本与聚类中心的距离,划分聚类。
更新聚类中心, 重新划分聚类。
重复第 4 步, 直到各类包含的样本不再改变或达到最大迭代次数, 输出结果。
对原样本数据计算均值和方差, 进行标准化, 使用 PCA 降维方法降至二维,打印散点图至二维平面中, 以不同的颜色区分不同类别的样本。
更换随机聚类中心, 重复上述第 3-5 步并打印散点图。
参数设置: k 分别取 2、3、5，最大迭代次数 maxiter 取 100

代码输出：样本二维图像，共6幅放在一起，三行两列，分别对应3个不同的k值，每个k值对应两次不同的初始聚类中心。
## DBSCAN聚类算法
采用数据集 “data/clustering2.csv”进行DBSCAN聚类算法实验。

数据集简介: 本数据来自 DIP 数据库, 为蛋白质相互作用网络数据, 网络 $G-$
$(E, V), \mathrm{E}$ 是网络中的节点的集合, 即蛋白质, $\mathrm{V}$ 是网络中的边的集合, 即蛋白质 -蛋白质相互作用。网络保存为邻接矩阵 $\mathrm{A}$ 的形式, 矩阵中元素 $A_{i j}$ 代表第 $\mathrm{i}$ 个蛋白质和第 $\mathrm{j}$ 个蛋白质是否存在相互作用, 若值为 1 , 则存在相互作用, 并有一条边连接这两个节点。矩阵中 $A_{i j}$ 和 $A_{j i}$ 对应同一条边。矩阵维度为 $1274 \times 1274$, 即共有 $1274$ 个蛋白质。存在相互作用的两个蛋白质之间的距离定为 $1$ , 无须计算。

本实验中, 邻域半径 $\varepsilon$ 设为 1 , 即密度直达的两个节点就是存在相互作用的两个蛋白质, 也就是说邻接矩阵中值为 1 的元素所在的行和列所对应的两个节点是密度直达的。
密度阈值 MinPts 设为 $10 、 15 、 20$ 。
依据密度聚类流程, 首先根据 MinPts 构造核心对象子集。
依次从中取出一个核心对象 (同时从子集中删去该节点), 将它的邻居加入待访问节点子集。依次从待访问节点子集中取出一个节点 $v$ (同时从子集中删去该节点), 若它是核心对象, 将它的邻居也放入待访问节点子集。重复这一过程, 直到待访问节点子集为空集, 此时所有访问过的结点构成一个聚类。
重复第 5 步, 直到核心对象子集为空集。
计算各个聚类的聚类系数：
$$coef=\frac{n_{edge}}{2\times n_{node}\times\left(n_{node-1}\right)}$$
保存所有聚类的聚类系数。

代码输出:不同密度阈值MinPts取值时对应的聚类数量和相应各个聚类的聚类系数

实验代码和结果

import pandas as pd
from matplotlib import pyplot as plt
import numpy as np
class KMeans:
    def distEclud(self, vecA, vecB):
        return np.linalg.norm(vecA - vecB)

    def initCentroids(self, dataSet, k):
        numSamples, dim = dataSet.shape
        centroids = np.zeros((k, dim))
        for i in range(k):
            index = int(np.random.uniform(0, numSamples))
            centroids[i, :] = dataSet[index, :]
        return centroids

    def kmeans(self, dataSet, k):
        numSamples = dataSet.shape[0]

        clusterAssment = np.mat(np.zeros((numSamples, 2)))
        clusterChanged = True

        centroids = self.initCentroids(dataSet, k)

        while clusterChanged:
            clusterChanged = False
            for i in range(numSamples):
                minDist = 100000.0
                minIndex = 0

                for j in range(k):
                    distance = self.distEclud(centroids[j, :], dataSet[i, :])
                    if distance < minDist:
                        minDist = distance
                        minIndex = j

                if clusterAssment[i, 0] != minIndex:
                    clusterChanged = True
                    clusterAssment[i, :] = minIndex, minDist ** 2

            for j in range(k):
                pointsInCluster = dataSet[np.nonzero(clusterAssment[:, 0].A == j)[0]]
                centroids[j, :] = np.mean(pointsInCluster, axis=0)

        return centroids, clusterAssment
def pca(data, k=2):
    mean_vals = np.mean(data, axis=0)
    mean_removed = data - mean_vals
    cov_mat = np.cov(mean_removed, rowvar=False)
    eig_vals, eig_vecs = np.linalg.eig(np.mat(cov_mat))
    eig_val_index = np.argsort(eig_vals)[::-1]
    eig_val_index = eig_val_index[:k]
    reduced_data = mean_removed * eig_vecs[:, eig_val_index]
    return reduced_data

data = pd.read_csv(r'E:\Google\data\clustering1.csv',sep=',',header=None)
data1=data.values
Max=data1.max()
Min=data1.min()
data1=(data1-Min)/(Max-Min)
kmeans_instance = KMeans()
centroids, cluster_assignment = kmeans_instance.kmeans(data1, k=5)
data=data.values
Mean=data.mean()
Std=data.std()
data=(data-Mean)/Std
rd=pca(data)
rd=np.column_stack((rd,cluster_assignment[:, 0]))
color=['r','g','b','c','m','y']
for i in range(len(rd)):
    plt.scatter(rd[i,0],rd[i,1],color=color[int(rd[i,2])])
plt.show()

k=2

k=3

k=5

import pandas as pd
import numpy as np
def DBSCAN(Matrix, eps, min_pts):
    Cluster = np.full(Matrix.shape[0], -1, dtype=int)
    cluster_id = 0

    for i in range(Matrix.shape[0]):
        if Cluster[i] != -1:
            continue

        neighbors = Neighbors(Matrix, i, eps)
        if len(neighbors) < min_pts:
            Cluster[i] = 0
        else:
            Expand(Matrix, Cluster, i, neighbors, cluster_id, eps, min_pts)
            cluster_id += 1
    return Cluster

def Neighbors(Matrix, center, eps):
    return np.where(Matrix[center] == 1)[0]

def Expand(Matrix, Cluster, center, neighbors, cluster_id, eps, min_pts):
    Cluster[center] = cluster_id

    i = 0
    while i < len(neighbors):
        n = neighbors[i]

        if Cluster[n] == -1:
            Cluster[n] = cluster_id
            new_neighbors = Neighbors(Matrix, n, eps)

            if len(new_neighbors) >= min_pts:
                neighbors = np.concatenate((neighbors, new_neighbors))

        i += 1

def Coefficient(Matrix, Cluster, cluster_id):
    cluster_nodes = np.where(Cluster == cluster_id)[0]

    if len(cluster_nodes) < 2:
        return 0.0

    internal_edges = np.sum(Matrix[np.ix_(cluster_nodes, cluster_nodes)])

    cluster_coefficient = 2.0 * internal_edges / (len(cluster_nodes) * (len(cluster_nodes) - 1))

    return cluster_coefficient
A = pd.read_csv(r'E:\Google\data\clustering2.csv',sep=',',header=None)
A=A.to_numpy()
min_pts_values = [10, 15, 20]
for min_pts in min_pts_values:
    Cluster = DBSCAN(A, eps=1, min_pts=min_pts)
    print(f'MinPts={min_pts}\n')
    unique_Cluster= np.unique(Cluster)
    for cluster in unique_Cluster:
        cluster_members = np.where(Cluster == cluster)[0]
        coefficient = Coefficient(A, Cluster, cluster)
        print(f'Cluster {cluster}: {coefficient}\n')

小结或讨论

k均值简单并且可以用于各种数据类型。它相当有效，尽管常常多次运行。k均值的某些变种（包括二分K均值）甚至更有效，并且不太受初始化问题的影响。然而，k均值并不适合所有的数据类型。它不能出来非球形簇、不同尺寸和不同密度的簇。最后，k均值仅限于具有中心（质心）概念的数据。

DBSCAN(Density-Based Spatial Clustering of Applications with Noise)是一个比较有代表性的基于密度的聚类算法。与划分和层次聚类方法不同，它将簇定义为密度相连的点的最大集合，能够把具有足够高密度的区域划分为簇，并可在噪声的空间数据库中发现任意形状的聚类（笔者认为是因为他不是基于距离的，基于距离的发现的是球状簇）。

该算法利用基于密度的聚类的概念，即要求聚类空间中的一定区域内所包含对象（点或其他空间对象）的数目不小于某一给定阈值。DBSCAN算法的显著优点是聚类速度快且能够有效处理噪声点和发现任意形状的空间聚类。但是由于它直接对整个数据库进行操作且进行聚类时使用了一个全局性的表征密度的参数，因此也具有两个比较明显的弱点：

（1）当数据量增大时，要求较大的内存支持I/O消耗也很大；

（2）当空间聚类的密度不均匀、聚类间距差相差很大时，聚类质量较差（有些簇内距离较小，有些簇内距离很大，但是Eps是确定的，所以，大的点可能被误判断为离群点或者边界点，如果Eps太大，那么小距离的簇内，可能会包含一些离群点或者边界点，KNN的k也存在同样的问题）。

（1）与K-MEANS比较起来，不需要输入要划分的聚类个数；

（2）聚类簇的形状没有偏倚；

（3）可以在需要时输入过滤噪声的参数；

Lilin

聚类实验