一、作业内容：

在二维空间内随机创建500个数据对象构建数据集D，每个维度的取值范围为[0,10]，设置三组不同的的(Eps，Minpts)值，使用DBSCAN算法对数据集D进行聚类分析。
作业要求：输出原始数据对象分布和聚类分布。可参考下图，不同聚类可用不同符号标注，如+，○，□，△等。

二、DBSCAN算法

1.主要概念

Ε邻域：给定对象半径为Ε内的区域称为该对象的Ε邻域；
核心对象：如果给定对象Ε领域内的样本点数大于等于MinPts，则称该对象为核心对象；
直接密度可达：对于样本集合D，如果样本点q在p的Ε领域内，并且p为核心对象，那么对象q从对象p直接密度可达。
密度可达：对于样本集合D，给定一串样本点p1,p2….pn，p= p1,q= pn,假如对象pi从pi-1直接密度可达，那么对象q从对象p密度可达。
密度相连：存在样本集合D中的一点o，如果对象o到对象p和对象q都是密度可达的，那么p和q密度相联。

2.算法执行步骤

（1）检测数据库中尚未检查过的对象p，如果p未被处理(归为某个簇或者标记为噪声)，则检查其邻域，若包含的对象数不小于minPts，建立新簇C，将其中的所有点加入候选集N；
（2）对候选集N 中所有尚未被处理的对象q，检查其邻域，若至少包含minPts个对象，则将这些对象加入N；如果q 未归入任何一个簇，则将q 加入C；
（3）重复步骤2)，继续检查N 中未处理的对象，当前候选集N为空；
（4）重复步骤1)~3)，直到所有对象都归入了某个簇或标记为噪声。

3.伪代码描述如下：

输入：数据对象集合D，半径Eps，密度阈值MinPts
输出：聚类C
DBSCAN（D, Eps, MinPts）
Begin
init C=0; //初始化簇的个数为0
for each unvisited point p in D
mark p as visited; //将p标记为已访问
N = getNeighbours (p, Eps);
if sizeOf(N) < MinPts then
mark p as Noise; //如果满足sizeOf(N) < MinPts，则将p标记为噪声
else
C= next cluster; //建立新簇C
ExpandCluster (p, N, C, Eps, MinPts);
end if
end for
End

其中ExpandCluster算法伪码如下：

ExpandCluster(p, N, C, Eps, MinPts)
add p to cluster C; //首先将核心点加入C
for each point p’ in N
mark p' as visited;
N’ = getNeighbours (p’, Eps); //对N邻域内的所有点在进行半径检查
if sizeOf(N’) >= MinPts then
N = N+N’; //如果大于MinPts，就扩展N的数目
end if
if p’ is not member of any cluster
add p’ to cluster C; //将p' 加入簇C
end if
end for
End ExpandCluster

《数据分析方法》

作业一：DBSCAN算法的实现

一、程序运行结果

1.原始数据集数据分布如图1所示。

2.三组不同Eps，Minpts值下聚类结果。

（1）Eps=0.3，Minpts=10时，聚类结果如图2所示。

（2）Eps=0.3，Minpts=1时，聚类结果如图3所示。

（3）Eps=0.03，Minpts=10时，聚类结果如图4所示。

二、聚类结果分析
1.随m值增大时，
聚类结果较差

2.随n值增大时，
要求较大的内存支持，I/0消耗也很大

三、程序主要代码

from sklearn import datasets
import numpy as np
import random
import matplotlib.pyplot as plt
import time
import copy


def find_neighbor(j, x, eps):
    N = list()
    for i in range(x.shape[0]):
        temp = np.sqrt(np.sum(np.square(x[j] - x[i])))  # 计算欧式距离
        if temp <= eps:
            N.append(i)
    return set(N)


def DBSCAN(X, eps, min_Pts):
    k = -1
    neighbor_list = []  # 用来保存每个数据的邻域
    omega_list = []  # 核心对象集合
    gama = set([x for x in range(len(X))])  # 初始时将所有点标记为未访问
    cluster = [-1 for _ in range(len(X))]  # 聚类
    marker = ["s" for _ in range(len(X))]
    for i in range(len(X)):
        neighbor_list.append(find_neighbor(i, X, eps))
        if len(neighbor_list[-1]) >= min_Pts:
            omega_list.append(i)  # 将样本加入核心对象集合
    omega_list = set(omega_list)  # 转化为集合便于操作
    while len(omega_list) > 0:
        gama_old = copy.deepcopy(gama)
        j = random.choice(list(omega_list))  # 随机选取一个核心对象
        k = k + 1
        Q = list()
        Q.append(j)
        gama.remove(j)
        while len(Q) > 0:
            q = Q[0]
            Q.remove(q)
            if len(neighbor_list[q]) >= min_Pts:
                delta = neighbor_list[q] & gama
                deltalist = list(delta)
                for i in range(len(delta)):
                    Q.append(deltalist[i])
                    gama = gama - delta
        Ck = gama_old - gama
        Cklist = list(Ck)
        for i in range(len(Ck)):
            cluster[Cklist[i]] = k
            marker[Cklist[i]] = "o"
        omega_list = omega_list - Ck
    return cluster, marker


centers = [[5, 5], [7, 6], [7, 3]]
# 产生的数据个数
n_samples = 500
# 生产数据:此实验结果受cluster_std的影响，或者说受eps 和cluster_std差值影响
X, lables_true = datasets.make_blobs(n_samples=n_samples, centers=centers, cluster_std=0.4,
                                     random_state=0, center_box=(0, 10))
# n_samples(int/array):如果参数为int，代表总样本数；如果参数为array-like，数组中的每个数代表每一簇的样本数。
# n_features(int):样本点的维度。
# centers(int):样本中心数。如果样本数为int且centers=None，生成三个样本中心；如果样本数（n_samples）为数组，则centers 要么为None，要么为数组的长度。
# cluster_std(float/sequence of floats):样本中，簇的标准差。
# center_box(pair of floats (min, max)):每个簇的上下限。
# shuffle(boolean):是否将样本打乱。
# random_state(int/RandomState instance /None):指定随机数种子，每个种子生成的序列相同，与minecraft地图种子同理。

plt.figure()
plt.scatter(X[:, 0], X[:, 1])
plt.show()
X = np.array(X)
eps = 0.3
min_Pts = 10
begin = time.time()
C, k = DBSCAN(X, eps, min_Pts)
end = time.time()
plt.figure()
plt.scatter(X[:, 0], X[:, 1], c=C)
plt.show()

附：MATLAB代码

clc;
clear all;
close all;
x=0+10*rand(500,2);
y=0+10*rand(1,500);
data1=x;
% 显示数据
plot(data1(:,1),data1(:,2),'bo');
hold on;
grid on;
data=data1; 
N=3;%设置聚类数目
[m,n]=size(data);
pattern=zeros(m,n+1);
center=zeros(N,n);%初始化聚类中心
pattern(:,1:n)=data(:,:);
for x=1:N
    center(x,:)=data( randi(500,1),:);%第一次随机产生聚类中心
end
while 1
distence=zeros(1,N);
num=zeros(1,N);
new_center=zeros(N,n);
 
for x=1:m
    for y=1:N
    distence(y)=norm(data(x,:)-center(y,:));%计算到每个类的距离
    end
    [~, temp]=min(distence);%求最小的距离
    pattern(x,n+1)=temp;         
end
k=0;
for y=1:N
    for x=1:m
        if pattern(x,n+1)==y
           new_center(y,:)=new_center(y,:)+pattern(x,1:n);
           num(y)=num(y)+1;
        end
    end
    new_center(y,:)=new_center(y,:)/num(y);
    if norm(new_center(y,:)-center(y,:))<0.10
        k=k+1;
    end
end
if k==N
     break;
else
     center=new_center;
end
end
[m, n]=size(pattern);
 
%最后显示聚类后的数据
figure;
hold on;
for i=1:m
    if pattern(i,n)==1 
         plot(pattern(i,1),pattern(i,2),'r*');
         plot(center(1,1),center(1,2),'ko');
    elseif pattern(i,n)==2
         plot(pattern(i,1),pattern(i,2),'g+');
         plot(center(2,1),center(2,2),'ko');
    elseif pattern(i,n)==3
         plot(pattern(i,1),pattern(i,2),'bs');
         plot(center(3,1),center(3,2),'ko');
    elseif pattern(i,n)==4
         plot(pattern(i,1),pattern(i,2),'yo');
         plot(center(4,1),center(4,2),'ko');
    else
         plot(pattern(i,1),pattern(i,2),'m.');
         plot(center(4,1),center(4,2),'ko');
    end
end
grid on;

运行结果