【无监督学习】DBSCAN 聚类步骤及matlab实现

DBSCAN 聚类算法

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）是一种基于密度的聚类算法，它能够发现任意形状的簇，并且可以有效处理数据集中的噪声点。核心概念：

• Epsilon邻域（ε-neighborhood）：对于数据集中的每个点 P，其 ε - 邻域定义为距离点 P 不超过指定半径 ε 的所有点的集合；
• 核心点（Core Point）：如果一个点在其 ε - 邻域内至少包含 MinPts 个点（包括自身），则该点被视为核心点。这里的 MinPts 是一个用户指定的参数，表示最小点数阈值；
• 边界点（Border Point）：如果一个点不是核心点，但它落在某个核心点的 ε - 邻域内，则该点被称为边界点；
• 噪音点（Noise Point）：既不是核心点也不是边界点的点被认为是噪音点或异常值。
  

特点：

• 无需预先指定簇的数量：与 K-means 等算法不同，DBSCAN 不需要用户事先确定要生成的簇数量；
• 发现任意形状的簇：由于 DBSCAN 基于密度的概念来定义簇，因此它可以识别出非球形或其他复杂形状的簇；
• 处理噪声点：DBSCAN 能够有效地识别并标记数据集中的噪声点，而不是强制将它们分配到某个簇中；
• 对输入参数敏感：DBSCAN 的表现高度依赖于 ε 和 MinPts 的选择，选择不当可能导致不理想的聚类结果，例如将一个簇分割成多个小簇或将多个簇合并为一个大簇。

DBSCAN适用于多种领域和情况，特别是当面对的数据具有以下特征时：

• 存在噪声的数据集：如传感器网络数据、图像处理中的边缘检测等；
• 具有复杂形状的簇：比如地理信息系统（GIS）中的空间数据分析、市场营销中的客户细分等；
• 高维数据：虽然 DBSCAN 在理论上可以应用于任何维度的数据，但在实践中，随着维度的增加，寻找合适的 ε 变得更加困难。不过，在某些特定的高维应用中，如生物信息学中的基因表达数据分析，DBSCAN 仍然被证明是有效的。

1.参数选择

DBSCAN 是一种基于密度的聚类算法，其核心参数 邻域半径（Eps） 和 核心点邻域内最少样本数（MinPts） 直接决定了聚类的效果。

1. 使用 K - 距离图确定 Eps

   • K - 距离图（k-Distance Graph） 是选择 Eps 的核心工具，步骤如下：
     • 计算 K - 距离：
       • 对每个样本 $p_i$，计算其到所有其他样本的欧氏距离；
       • 将距离从小到大排序，取第 $k$ 个距离（即第 $k$ 近邻的距离），记为 $d(k)$；
       • k 的取值：通常设为 $MinPts-1$（例如，若 $MinPts=5$，则 $k=4$，排除自身后第 5 个点为第 4 距离）。
   • 绘制 K - 距离图：
     • 将所有样本的 $d(k)$ 升序排列并绘制曲线；
     • 拐点选择：曲线中斜率突然增大的点（突变点）对应的距离值即为 $Eps$ 的合理估计。

2. 设定 MinPts 的经验法则

   • 最小值：$MinPts\ge 数据维度+1$（如 3 维数据，$MinPts\ge 4$）。
   • 常用值：$MinPts=2\times 数据维度$（平衡灵敏度与抗噪性）。
   • 调整策略：
     • 数据噪声多或维度高时，适当增大 $MinPts$。
     • 若数据分布稀疏，可略微降低 $MinPts$。

matlab实现：
全局历史均值法：每个数据点与其之前所有数据的平均值比较，若当前值超过历史均值的两倍，则标记为突增点。

clc;
clear;
% 输入数据（示例）
data = [0.01, 0.004, 0.004, 0.004, 0.005, 0.004, 0.007, 0.009, 0.008, 0.012, ...
        0.010, 0.007, 0.008, 0.006, 0.008, 0.006, 0.006, 0.007, 0.009, 0.011, ...
        0.015, 0.026, 0.058, 0.381]; 

% 初始化突增点索引
spike_points = [];

% 遍历数据（从第二个点开始）
for i = 2:length(data)
    % 计算当前点之前所有数据的均值
    historical_mean = mean(data(1:i-1));
    
    % 判断当前点是否超过历史均值的两倍
    if data(i) > 2 * historical_mean
        spike_points = [spike_points, i];
    end
end

% 输出结果
if ~isempty(spike_points)
    fprintf('突增点索引：%s\n', num2str(spike_points));
else
    fprintf('未检测到突增点\n');
end

% 可视化
figure;
plot(data, 'o-', 'DisplayName', '原始数据'); hold on;
plot(spike_points, data(spike_points), 'rx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName', '突增点');
xlabel('数据点索引');
ylabel('数值');
title('全局历史均值法检测结果');
legend;
grid on;

2.算法步骤

1. 输入数据与参数

   • 输入：数据矩阵 $X\in R^{n\times m}$，$n$ 为样本数，$m$ 为特征数；
   • 参数：
     • Eps：邻域半径，决定密度范围；
     • MinPts：核心点邻域内最少样本数。

2. 计算距离矩阵

   • 目标：计算所有样本间的欧氏距离；
   • 公式：$$d_{ij}=\sqrt{\sum _{k=1}^m(x_{ik}-x_{jk}{)}^2}$$
   • 优化：预计算距离矩阵以加速邻域查询。

3. 标记所有点为未访问

   • 初始化：创建 visited 数组，所有元素设为 false。

4. 遍历未访问点

   • 循环：逐个检查未访问点，判断其是否为核心点。

5. 检查是否为核心点

   • 邻域查询：统计当前点 $p$ 的 $Eps$ 邻域内样本数；
   • 判断条件：若 $邻域样本数\ge MinPts$，标记为核心点，否则标记为噪声。

6. 扩展聚类

   • 广度优先搜索 (BFS)：
     • 将核心点邻域内所有点加入队列；
     • 遍历队列中的点：
       • 若未访问，标记为已访问；
       • 若为核心点，将其邻域点加入队列。
     • 将访问过的点归入当前聚类。

7. 标记噪声

   • 非核心点且未被任何聚类包含：标记为噪声点（簇标签为 -1）。

3.MATLAB 实现

%% DBSCAN 对门店聚类分析
clc; clear; close all;

%% 生成模拟地理销售数据
rng(42); % 固定随机种子保证可复现性
n = 1200; % 总样本数

% 环形区域（50%数据）
n_ring = round(n*0.5);
theta = 2*pi*rand(n_ring, 1); 
r = 50 + 5*randn(n_ring, 1);   % 半径：N(50,5) 增加方差
X_ring = [r.*cos(theta), r.*sin(theta)];
sales_ring = 40 + 10*randn(n_ring, 1); 

% 中心区域（约41.67%数据）
n_center = round(n*10/24);
X_center = 6*randn(n_center, 2); % 中心坐标：N(0,6) 扩大分布范围
sales_center = 80 + 15*randn(n_center, 1); 

% 噪声点（约8.33%数据）
n_noise = round(n - n_ring - n_center);
X_noise = 100*rand(n_noise, 2) - 50;     
sales_noise = 20 + 10*rand(n_noise, 1);  % 修正为U(20,30)

% 合并数据
X = [X_ring; X_center; X_noise];
sales = [sales_ring; sales_center; sales_noise];

%数据打乱重排
shuffled_idx = randperm(n);
X = X(shuffled_idx, :);
sales = sales(shuffled_idx);

%% 数据标准化（仅地理位置）
X_scaled = zscore(X); % 只对地理坐标标准化
sales_normalized = (sales - min(sales)) / (max(sales) - min(sales)); % 销售额归一化[0,1]

%% K-距离图确定Eps参数
k = 5; % MinPts-1（后续DBSCAN的MinPts设为k+1=6）
k_distances = zeros(size(X_scaled,1),1);

% 计算每个点的第k近邻距离
for i = 1:size(X_scaled,1)
    distances = sort(pdist2(X_scaled(i,:), X_scaled, 'euclidean'));
    k_distances(i) = distances(k+1); % 排除自身（第0个是自己）
end

% 排序并绘制K-距离图
sorted_k_dist = sort(k_distances);
x = linspace(1, n, n);
x_new = linspace(min(x), max(x), n/10);
sorted_k_dist = interp1(x, sorted_k_dist, x_new);
gradient = diff(sorted_k_dist,1); % 计算梯度变化

% 初始化突增点索引
spike_points = [];
% 遍历数据（从第二个点开始）
for i = 2:length(gradient)
    % 计算当前点之前所有数据的均值
    historical_mean = mean(abs(gradient(1:i-1)));
    
    % 判断当前点是否超过历史均值的五倍
    if abs(gradient(i)) > 5 * historical_mean
        spike_points = [spike_points, i];
    end
end

% 输出结果
if ~isempty(spike_points)
    fprintf('突增点索引：%s\n', num2str(spike_points));
else
    fprintf('未检测到突增点\n');
end

% 可视化
figure;
plot(gradient, 'o-', 'DisplayName', '原始数据'); hold on;
plot(spike_points, gradient(spike_points), 'rx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName', '突增点');
xlabel('数据点索引');
ylabel('数值');
title('全局历史均值法检测结果');
legend;
grid on;

%% DBSCAN参数设置（根据K-距离图结果）
Eps = sorted_k_dist(spike_points(1));    % 使用自动检测的Eps
MinPts = k + 1;    % 根据k值设置MinPts

%% 优化后的DBSCAN实现
n_points = size(X_scaled,1);
D = squareform(pdist(X_scaled)); % 使用squareform更高效

visited = false(n_points, 1);
clusters = zeros(n_points, 1);
cluster_id = 0;

for i = 1:n_points
    if ~visited(i)
        visited(i) = true;
        
        % 查找Eps邻域（预计算距离更高效）
        neighbors = find(D(i,:) <= Eps);
        
        if numel(neighbors) >= MinPts
            cluster_id = cluster_id + 1;
            clusters(i) = cluster_id;
            
            % 使用队列优化BFS
            queue = neighbors;
            head = 1;
            tail = numel(queue);
            
            while head <= tail
                p = queue(head);
                if ~visited(p)
                    visited(p) = true;
                    p_neighbors = find(D(p,:) <= Eps);
                    
                    if numel(p_neighbors) >= MinPts
                        % 去重合并新邻域
                        new_nodes = setdiff(p_neighbors, queue);
                        queue = [queue, new_nodes]; %#ok<AGROW>
                        tail = tail + numel(new_nodes);
                    end
                end
                
                if clusters(p) == 0
                    clusters(p) = cluster_id;
                end
                
                head = head + 1;
            end
        else
            clusters(i) = -1;
        end
    end
end

%% 可视化与分析（新增三维可视化）
figure;

% 三维地理-销售额分布
subplot(1,2,1);
scatter3(X(:,1), X(:,2), sales, 15, clusters, 'filled');
colormap(jet);
xlabel('经度');
ylabel('纬度');
zlabel('销售额');
title('地理-销售额三维分布');
view(-30, 30);

% 二维聚类结果
subplot(1,2,2);
valid_clusters = clusters(clusters ~= -1);
if ~isempty(valid_clusters)
    gscatter(X(:,1), X(:,2), clusters, 'brgmck', '.', 15);
else
    scatter(X(:,1), X(:,2), 15, 'b.');
end
hold on;
plot(X(clusters==-1,1), X(clusters==-1,2), 'k.', 'MarkerSize', 5);
title('DBSCAN聚类结果');
xlabel('经度');
ylabel('纬度');
legend_text = arrayfun(@(x)sprintf('簇%d',x), unique(clusters(clusters>0)), 'UniformOutput',false);
legend([ '噪声';legend_text], 'Location', 'best');

%% 增强统计分析
fprintf('=== 聚类质量评估 ===\n');
fprintf('参数配置: Eps=%.2f, MinPts=%d\n', Eps, MinPts);
fprintf('有效簇数量: %d\n', cluster_id);
fprintf('噪声点比例: %.1f%%\n', mean(clusters==-1)*100);

if cluster_id > 0
    % 簇大小分布
    cluster_sizes = arrayfun(@(x)sum(clusters==x), 1:cluster_id);
    fprintf('\n=== 簇规模分布 ===\n');
    fprintf('平均规模: %.1f ± %.1f\n', mean(cluster_sizes), std(cluster_sizes));
    
    % 地理密度分析
    cluster_density = cluster_sizes ./ (pi * Eps^2); % 单位面积密度
    fprintf('平均密度: %.1f 点/单位面积\n', mean(cluster_density));
    
    % 销售额分析
    fprintf('\n=== 销售额分析 ===\n');
    sales_table = table();
    for c = 1:cluster_id
        mask = clusters == c;
        sales_table.Cluster(c) = c;
        sales_table.Mean(c) = mean(sales(mask));
        sales_table.Median(c) = median(sales(mask));
        sales_table.Std(c) = std(sales(mask));
        sales_table.IQR(c) = iqr(sales(mask));
    end
    disp(sales_table);
end

%% 业务建议生成系统
fprintf('\n=== 自适应业务策略 ===\n');
if cluster_id == 0
    fprintf('⚠️ 未发现有效聚类，建议重新调整参数或检查数据质量\n');
else
    % 按销售额排序簇
    [~, sorted_idx] = sort(sales_table.Mean, 'descend');
    
    for i = 1:numel(sorted_idx)
        c = sorted_idx(i);
        fprintf('\n🏆 高价值集群 %d（Top %d）\n', c, i);
        fprintf('  规模: %d 门店\n', sum(clusters==c));
        fprintf('  销售额: %.1f ± %.1f 万元\n', sales_table.Mean(c), sales_table.Std(c));
        
        if i == 1
            fprintf('  策略: 设立区域总部，优先资源投放\n');
        else
            fprintf('  策略: 优化运营流程，开展促销活动\n');
        end
    end
    
    if any(clusters==-1)
        fprintf('\n🚧 异常门店 (%d 家)\n', sum(clusters==-1));
        fprintf('  策略: 开展专项审计，制定关停/改造计划\n');
    end
end

参考资料

[1] 基于密度的聚类 DBSCAN 解释与实例计算_哔哩哔哩_bilibili
[2] 快速学会聚类算法系列之DBSCAN（附matlab代码)哔哩哔哩_bilibili
[3] 5-DBSCAN工作流程_哔哩哔哩_bilibili