聚类算法-Cluster Algorithm

总结

• 聚类算法对比总结, Overview of clustering methods
• 

sklearn聚类算法

• 参考：
  • sklearn聚类算法详解
  • scikit-learn中的无监督聚类算法
• scikit-learn主要由分类、回归、聚类和降维四大部分组成，其中分类和回归属于有监督学习范畴，聚类属于无监督学习范畴，降维适用于有监督学习和无监督学习。scikit-learn的结构示意图如下所示：

• 

• scikit-learn中的聚类算法主要有：
  • K-Means(cluster.KMeans)
  • AP聚类（cluster.AffinityPropagation）
  • 均值漂移（cluster.MeanShift）
  • 层次聚类（cluster.AgglomerativeClustering）
  • DBSCAN（cluster.DBSCAN）
  • BRICH（cluster.Brich）
  • 谱聚类（cluster.Spectral.Clustering）
  • 高斯混合模型（GMM）∈期望最大化（EM）算法（mixture.GaussianMixture）

1 KMeans

1.1 算法描述

1. 随机选择k个中心
2. 遍历所有样本，把样本划分到距离最近的一个中心
3. 划分之后就有K个簇，计算每个簇的平均值作为新的质心
4. 重复步骤2，直到达到停止条件

• 停止条件：
  • 聚类中心不再发生变化；所有的距离最小；迭代次数达到设定值，
• 代价函数：误差平方和（SSE）

1.2 算法优缺点

优点：

• 算法容易理解，聚类效果不错
• 具有出色的速度
• 当簇近似高斯分布时，效果比较好

缺点：

• 需要自己确定K值,k值的选定是比较难确定
• 对初始中心点敏感
• 不适合发现非凸形状的簇或者大小差别较大的簇
• 特殊值/离群值对模型的影响比较大
• 从数据先验的角度来说，在 Kmeans 中,我们假设各个 cluster 的先验概率是一样的,但是各个 cluster 的数据量可能是不均匀的。举个例子,cluster A 中包含了10000个样本,cluster B 中只包含了100个。那么对于一个新的样本,在不考虑其与A cluster、 B cluster 相似度的情况,其属于 cluster A 的概率肯定是要大于 cluster B的。

1.3 效果评价

• 从簇内的稠密程度和簇间的离散程度来评估聚类的效果。
• 常见的方法有轮廓系数Silhouette Coefficient和Calinski-Harabasz Index

1.3.1轮廓系数

• 轮廓系数 （Silhouette Coefficient）结合了聚类的凝聚度（Cohesion）和分离度（Separation），用于评估聚类的效果。该值处于-1~1之间，值越大，表示聚类效果越好。具体计算方法如下：
  1. 对于第i个元素x_i，计算x_i与其同一个簇内的所有其他元素距离的平均值，记作a_i，用于量化簇内的凝聚度。
  2. 选取x_i外的一个簇b，计算x_i与b中所有点的平均距离，遍历所有其他簇，找到最近的这个平均距离,记作b_i，用于量化簇之间分离度。
  3. 对于元素x_i，轮廓系数s_i = (b_i – a_i)/max(a_i,b_i)
  4. 计算所有x的轮廓系数，求出平均值即为当前聚类的整体轮廓系数

先是计算每一个样本的轮廓系数，然后计算所有样本的轮廓系数，求平均值作为整体轮廓系数

从上面的公式，不难发现若s_i小于0，a_i > b_i, 说明x_i与其簇内元素的平均距离大于最近的其他簇，表示聚类效果不好。如果a_i趋于0，或者b_i足够大，那么s_i趋近与1，说明聚类效果比较好。

1.3.2 Calinski-Harabasz Index

这个不知道怎么翻译，估计是两个人名。

类别内部数据的协方差越小越好，类别之间的协方差越大越好，这样的Calinski-Harabasz分数会高。

在scikit-learn中， Calinski-Harabasz Index对应的方法是metrics.calinski_harabaz_score。

import numpy as np
from sklearn.cluster import KMeans
kmeans_model = KMeans(n_clusters=3, random_state=1).fit(X)
labels = kmeans_model.labels_metrics.calinski_harabaz_score(X, labels)

参考博客： KMeans

1.4 K值确定

• 结合业务分析，确定需要分类的个数，这种情况往往有业务上聚类的个数范围
• 手肘原则，选定不同的K值，计算每个k值时的代价函数。Kmeans聚类的效果评估方法是SSE，是计算所有点到相应簇中心的距离均值，当然，k值越大 SSE越小，我们就是要求出随着k值的变化SSE的变化规律，找到SSE减幅最小的k值，这时k应该是相对比较合理的值。

如图，在k=3之后，代价函数变化缓慢，选择聚类的个数为3

• 算法特点
  • Distances between points（点之间的距离）
• 优点：
  • 算法容易理解，聚类效果不错
  • 具有出色的速度：O(NKt)
  • 当簇近似高斯分布时，效果比较好
• 缺点：
  • 需要人工预先确定初试K值，且该值和真是的数据分布未必吻合
  • 对初始中心点敏感
  • 不适合发现非凸形状的簇或者大小差别较大的簇
  • 特殊值/离散值（噪点）对模型的影响比较大
  • 算法只能收敛到局部最优，效果受初始值影响很大
  • 从数据先验的角度来说，在 Kmeans 中,我们假设各个 cluster 的先验概率是一样的,但是各个 cluster 的数据量可能是不均匀的。举个例子,cluster A 中包含了10000个样本,cluster B 中只包含了100个。那么对于一个新的样本,在不考虑其与A cluster、 B cluster 相似度的情况,其属于 cluster A 的概率肯定是要大于 cluster B的。
• 适用场景
  • 通用, 均匀的 cluster size（簇大小）, flat geometry（平面几何）, 不是太多的 clusters（簇）
  • 非常大的 n_samples、中等的 n_clusters 使用 MiniBatch code
  • 样本量<10K时使用k-means，>=10K时用MiniBatchKMeans
  • 不太适用于离散分类
• 代码

print(__doc__)

# Author: Phil Roth <[email protected]>
# License: BSD 3 clause

import numpy as np
import matplotlib.pyplot as plt

from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

plt.figure(figsize=(12, 12))

n_samples = 1500
random_state = 170
X, y = make_blobs(n_samples=n_samples, random_state=random_state)

# Incorrect number of clusters
y_pred = KMeans(n_clusters=2, random_state=random_state).fit_predict(X)

plt.subplot(221)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.title("Incorrect Number of Blobs")

# Anisotropicly distributed data
transformation = [[0.60834549, -0.63667341], [-0.40887718, 0.85253229]]
X_aniso = np.dot(X, transformation)
y_pred = KMeans(n_clusters=3, random_state=random_state).fit_predict(X_aniso)

plt.subplot(222)
plt.scatter(X_aniso[:, 0], X_aniso[:, 1], c=y_pred)
plt.title("Anisotropicly Distributed Blobs")

# Different variance
X_varied, y_varied = make_blobs(n_samples=n_samples,
                                cluster_std=[1.0, 2.5, 0.5],
                                random_state=random_state)
y_pred = KMeans(n_clusters=3, random_state=random_state).fit_predict(X_varied)

plt.subplot(223)
plt.scatter(X_varied[:, 0], X_varied[:, 1], c=y_pred)
plt.title("Unequal Variance")

# Unevenly sized blobs
X_filtered = np.vstack((X[y == 0][:500], X[y == 1][:100], X[y == 2][:10]))
y_pred = KMeans(n_clusters=3,
                random_state=random_state).fit_predict(X_filtered)

plt.subplot(224)
plt.scatter(X_filtered[:, 0], X_filtered[:, 1], c=y_pred)
plt.title("Unevenly Sized Blobs")

plt.show()

2 DBSCAN

• DBSCAN（Density-Based Spatial Clustering of Application with Noise）基于密度的空间聚类算法。
• 两个参数：
  • • Eps邻域半径(epsilon,小量，小的值）
  • • MinPts(minimum number of points required to form a cluster)定义核心点时的阈值。

3个点

• 核心点：对应稠密区域内部的点
• 边界点：对应稠密区域边缘的点
• 噪音点：对应稀疏区域中的点

上图红色为核心点，黄色为边界点，蓝色为噪音点

几个概念：

• 核心对象：对于任一样本点，如果其Eps邻域内至少包含MinPts个样本点，那么这个样本点就是核心对象。（一个点）
• 直接密度可达：如果一个样本点p处于一个核心对象q的Eps邻域内，则称该样本点p从对象q出发时是直接密度可达的。

• 密度相连：对于样本点p和q，如果存在核心对象m,使得p、p均由m直接密度可达，则称p和q密度相连。

• DBSCAN的聚类是一个不断生长的过程。先找到一个核心对象，从整个核心对象出发，找出它的直接密度可达的对象，再从这些对象出发，寻找它们直接密度可达的对象，一直重复这个过程，直至最后没有可寻找的对象了，那么一个簇的更新就完成了。也可以认为，簇是所有密度可达的点的集合。

• DBSCAN核心思想：从某个选定的核心点出发，不断向密度可达的区域扩张，从而得到一个包含核心点和边界点的最大化区域，区域中任意两点密度相连。

优点：

• 不需要指定cluster的数目
• 聚类的形状可以是任意的
• 能找出数据中的噪音，对噪音不敏感
• 算法应用参数少，只需要两个
• 聚类结果几乎不依赖于节点的遍历顺序

缺点：

• 如果样本集较大时，聚类收敛时间较长，此时可以对搜索最近邻时建立的KD树或者球树进行规模限制来改进。
• DBSCAN算法的聚类效果依赖于距离公式的选取，实际中常用的距离是欧几里得距离，由于‘维数灾难’，距离的度量标准已变得不再重要。（分类器的性能随着特征数量的增加而不断提升，但过了某一值后，性能不升反而下降，这种现象称为维数灾难。对于维度灾难的理解： 维度灾难的理解 ）
• 不适合数据集中密度差异很大的情形，因为这种情形，参数Eps,MinPts不好选择。（个人理解，如果是密度大的，你选一个小的邻域半径就可以把这些数据点聚类，但对于那些密度小的数据点，你设置的小的邻域半径，并不能把密度小的这些点给全部聚类。）

聚类形状可以是任意的，来个图直观感觉一下：

在sklearn中的应用

from sklearn.cluster import DBSCAN
DBSCAN(eps=0.5,  # 邻域半径
min_samples=5,    # 最小样本点数，MinPts
metric='euclidean',
metric_params=None,
algorithm='auto', # 'auto','ball_tree','kd_tree','brute',4个可选的参数 寻找最近邻点的算法，例如直接密度可达的点
leaf_size=30, # balltree,cdtree的参数
p=None, # 
n_jobs=1)

• 完整版代码

from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import DBSCAN
import numpy as np
import matplotlib.pyplot as plt
from itertools import cycle  ##python自带的迭代器模块
from sklearn.preprocessing import StandardScaler

##产生随机数据的中心
centers = [[1, 1], [-1, -1], [1, -1]]
##产生的数据个数
n_samples=750
##生产数据:此实验结果受cluster_std的影响，或者说受eps 和cluster_std差值影响
X, lables_true = make_blobs(n_samples=n_samples, centers= centers, cluster_std=0.4,
                  random_state =0)


##设置分层聚类函数
db = DBSCAN(eps=0.3, min_samples=10)
##训练数据
db.fit(X)
##初始化一个全是False的bool类型的数组
core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
'''
   这里是关键点(针对这行代码：xy = X[class_member_mask & ~core_samples_mask])：
   db.core_sample_indices_  表示的是某个点在寻找核心点集合的过程中暂时被标为噪声点的点(即周围点
   小于min_samples)，并不是最终的噪声点。在对核心点进行联通的过程中，这部分点会被进行重新归类(即标签
   并不会是表示噪声点的-1)，也可也这样理解，这些点不适合做核心点，但是会被包含在某个核心点的范围之内
'''
core_samples_mask[db.core_sample_indices_] = True

##每个数据的分类
lables = db.labels_

##分类个数：lables中包含-1，表示噪声点
n_clusters_ =len(np.unique(lables)) - (1 if -1 in lables else 0)

##绘图
unique_labels = set(lables)
'''
   1)np.linspace 返回[0,1]之间的len(unique_labels) 个数
   2)plt.cm 一个颜色映射模块
   3)生成的每个colors包含4个值，分别是rgba
   4)其实这行代码的意思就是生成4个可以和光谱对应的颜色值
'''
colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))

plt.figure(1)
plt.clf()


for k, col in zip(unique_labels, colors):
    ##-1表示噪声点,这里的k表示黑色
    if k == -1:
        col = 'k'

    ##生成一个True、False数组，lables == k 的设置成True
    class_member_mask = (lables == k)

    ##两个数组做&运算，找出即是核心点又等于分类k的值  markeredgecolor='k',
    xy = X[class_member_mask & core_samples_mask]
    plt.plot(xy[:, 0], xy[:, 1], 'o', c=col,markersize=14)
    '''
       1)~优先级最高，按位对core_samples_mask 求反，求出的是噪音点的位置
       2)& 于运算之后，求出虽然刚开始是噪音点的位置，但是重新归类却属于k的点
       3)对核心分类之后进行的扩展
    '''
    xy = X[class_member_mask & ~core_samples_mask]
    plt.plot(xy[:, 0], xy[:, 1], 'o', c=col,markersize=6)

plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()

3 OPTICS

• 是基于密度的聚类算法，OPTICS（Ordering Point To Idenfy the Cluster Structure）,不显式地生成数据聚类，只是对数据对象集合中的对象进行排序，得到一个有序的对象列表。
• 核心距离（core-distance） 给定参数eps，MinPts，使得某个样本点成为核心对象（核心点）的最小邻域半径，这个最小邻域半径为该样本点的核心距离。

在DBSCAN中，给定领域半径eps和MinPts可以确定一个核心对象，如果eps比较大，在核心对象的邻域半径eps内的点的总数就会大于你所设定的MinPts，所以核心距离就是一个核心点在满足MinPts时的一个最小邻域半径。

• 可达距离（reachability-distance）

rd(y,x)表示使得‘x为核心点’且‘y从x直接密度可达’同时成立的最小邻域半径。

参考资料

4 Spectral Clustering 谱聚类

1）概述

• Spectral Clustering(SC,即谱聚类)，是一种基于图论的聚类方法,它能够识别任意形状的样本空间且收敛于全局最有解，其基本思想是利用样本数据的相似矩阵进行特征分解后得到的特征向量进行聚类.它与样本特征无关而只与样本个数有关。
• 基本思路：将样本看作顶点,样本间的相似度看作带权的边,从而将聚类问题转为图分割问题:找到一种图分割的方法使得连接不同组的边的权重尽可能低(这意味着组间相似度要尽可能低),组内的边的权重尽可能高(这意味着组内相似度要尽可能高).

2）图解过程

如上图所示，断开虚线，六个数据被聚成两类。

3）Spectral Clustering算法函数

a）核心函数：sklearn.cluster.SpectralClustering

因为是基于图论的算法，所以输入必须是对称矩阵。

b）主要参数(参数较多， 详细参数 )

• n_clusters：聚类的个数。（官方的解释：投影子空间的维度）
• affinity：核函数，默认是’rbf’，可选：”nearest_neighbors”，”precomputed”,”rbf”或sklearn.metrics.pairwise_kernels支持的其中一个 - 内核之一。
• gamma :affinity指定的核函数的内核系数，默认1.0

c）主要属性

• labels_ ：每个数据的分类标签

• 算法特点
  • Graph distance (e.g. nearest-neighbor graph)（图形距离（例如最近邻图））
• 适用场景
  • 几个簇，均匀的簇大小，非平面几何
  • 中等的 n_samples, 小的 n_clusters
• 代码

from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import spectral_clustering
import numpy as np
import matplotlib.pyplot as plt
from sklearn import metrics
from itertools import cycle  ##python自带的迭代器模块

##产生随机数据的中心
centers = [[1, 1], [-1, -1], [1, -1]]
##产生的数据个数
n_samples=3000
##生产数据
X, lables_true = make_blobs(n_samples=n_samples, centers= centers, cluster_std=0.6,
                  random_state =0)

##变换成矩阵，输入必须是对称矩阵
metrics_metrix = (-1 * metrics.pairwise.pairwise_distances(X)).astype(np.int32)
metrics_metrix += -1 * metrics_metrix.min()
##设置谱聚类函数
n_clusters_= 4
lables = spectral_clustering(metrics_metrix,n_clusters=n_clusters_)

##绘图
plt.figure(1)
plt.clf()

colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
for k, col in zip(range(n_clusters_), colors):
    ##根据lables中的值是否等于k，重新组成一个True、False的数组
    my_members = lables == k
    ##X[my_members, 0] 取出my_members对应位置为True的值的横坐标
    plt.plot(X[my_members, 0], X[my_members, 1], col + '.')

plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()

5 Hierarchical Clustering 层次聚类

1）概述

• Hierarchical Clustering(层次聚类)：就是按照某种方法进行层次分类，直到满足某种条件为止。
• 主要分成两类：
  • a）凝聚：从下到上。首先将每个对象作为一个簇，然后合并这些原子簇为越来越大的簇，直到所有的对象都在一个簇中，或者某个终结条件被满足。
  • b）分裂：从上到下。首先将所有对象置于同一个簇中，然后逐渐细分为越来越小的簇，直到每个对象自成一簇，或者达到了某个终止条件。（较少用）

2）算法步骤（凝聚）

• a）将每个对象归为一类, 共得到N类, 每类仅包含一个对象. 类与类之间的距离就是它们所包含的对象之间的距离.
• b）找到最接近的两个类并合并成一类, 于是总的类数少了一个.
• c）重新计算新的类与所有旧类之间的距离.
• d）重复第2步和第3步, 直到最后合并成一个类为止(此类包含了N个对象).

3）图解过程

4）Hierarchical Clustering算法函数

• a）sklearn.cluster.AgglomerativeClustering
• b）主要参数( 详细参数 )
  • n_clusters：聚类的个数
  • linkage：指定层次聚类判断相似度的方法，有以下三种：
  • ward：组间距离等于两类对象之间的最小距离。（即single-linkage聚类）
  • average：组间距离等于两组对象之间的平均距离。（average-linkage聚类）
  • complete：组间距离等于两组对象之间的最大距离。（complete-linkage聚类）
• c）主要属性
  • labels_： 每个数据的分类标签

参考资料 聚类算法

• 算法特点
  • Distances between points（点之间的距离）
• 适用场景
  • 很多的簇，可能连接限制
  • 大的 n_samples 和 n_clusters

代码

# Authors: Gael Varoquaux
# License: BSD 3 clause (C) INRIA 2014

print(__doc__)
from time import time

import numpy as np
from scipy import ndimage
from matplotlib import pyplot as plt

from sklearn import manifold, datasets

digits = datasets.load_digits(n_class=10)
X = digits.data
y = digits.target
n_samples, n_features = X.shape

np.random.seed(0)

def nudge_images(X, y):
    # Having a larger dataset shows more clearly the behavior of the
    # methods, but we multiply the size of the dataset only by 2, as the
    # cost of the hierarchical clustering methods are strongly
    # super-linear in n_samples
    shift = lambda x: ndimage.shift(x.reshape((8, 8)),
                                  .3 * np.random.normal(size=2),
                                  mode='constant',
                                  ).ravel()
    X = np.concatenate([X, np.apply_along_axis(shift, 1, X)])
    Y = np.concatenate([y, y], axis=0)
    return X, Y


X, y = nudge_images(X, y)


#----------------------------------------------------------------------
# Visualize the clustering
def plot_clustering(X_red, labels, title=None):
    x_min, x_max = np.min(X_red, axis=0), np.max(X_red, axis=0)
    X_red = (X_red - x_min) / (x_max - x_min)

    plt.figure(figsize=(6, 4))
    for i in range(X_red.shape[0]):
        plt.text(X_red[i, 0], X_red[i, 1], str(y[i]),
                 color=plt.cm.nipy_spectral(labels[i] / 10.),
                 fontdict={'weight': 'bold', 'size': 9})

    plt.xticks([])
    plt.yticks([])
    if title is not None:
        plt.title(title, size=17)
    plt.axis('off')
    plt.tight_layout(rect=[0, 0.03, 1, 0.95])

#----------------------------------------------------------------------
# 2D embedding of the digits dataset
print("Computing embedding")
X_red = manifold.SpectralEmbedding(n_components=2).fit_transform(X)
print("Done.")

from sklearn.cluster import AgglomerativeClustering

for linkage in ('ward', 'average', 'complete', 'single'):
    clustering = AgglomerativeClustering(linkage=linkage, n_clusters=10)
    t0 = time()
    clustering.fit(X_red)
    print("%s :\t%.2fs" % (linkage, time() - t0))

    plot_clustering(X_red, clustering.labels_, "%s linkage" % linkage)

plt.show()

AP聚类

====

2.1 简介

• AP(Affinity Propagation)通常被翻译为近邻传播算法或者亲和力传播算法，是在2007年的Science杂志上提出的一种新的聚类算法。AP算法的基本思想是将全部数据点都当作潜在的聚类中心(称之为exemplar)，然后数据点两两之间连线构成一个网络(相似度矩阵)，再通过网络中各条边的消息(responsibility和availability)传递计算出各样本的聚类中心。

2.2 算法原理

概念：

• 吸引度（responsibility）矩阵R： 其中r(i,k)描述了数据对象k适合作为数据对象i的聚类中心的程度，表示的是从i到k的消息。
• 归属度（availability）矩阵A： 其中a(i,k)描述了数据对象i选择数据对象k作为其据聚类中心的适合程度，表示从k到i的消息。
• 相似度（similarity）矩阵S： 通常S（i，j）取i，j的欧氏距离的负值，当i=j时，通常取整个矩阵的最小值或者中位数(Scikit-learn中默认为中位数)，取得值越大则最终产生的类数量越多。

算法步骤

1. 算法初始，吸引度矩阵和归属度矩阵均初始化为0矩阵。
2. 更新吸引度矩阵
   • 
3. 更新归属度矩阵
   • 
4. 根据衰减系数λ对两个公式进行衰减
   • 
5. 重复步骤2/3/4直至矩阵稳定或者达到最大迭代次数，算法结束。最终取a+r最大的k作为聚类中心。

2.3 算法特点

• Graph distance (e.g. nearest-neighbor graph)（图形距离（例如，最近邻图））

• 优点：
  • 与其他聚类算法不同,AP聚类不需要指定K(经典的K-Means)或者是其他描述聚类个数的参数
  • 一个聚类中最具代表性的点在AP算法中叫做Examplar,与其他算法中的聚类中心不同,examplar是原始数据中确切存在的一个数据点,而不是由多个 数据点求平均而得到的聚类中心
  • 多次执行AP聚类算法,得到的结果完全一样的，即不需要进行随机选取初值步骤.
  • AP算法相对于Kmeans优势是不需要指定聚类数量,对初始值不敏感
  • 模型对数据的初始值不敏感。
  • 对初始相似度矩阵数据的对称性没有要求。
  • 相比与k-centers聚类方法，其结果的平方差误差较小。
• 缺点：
  • AP算法需要事先计算每对数据对象之间的相似度，如果数据对象太多的话，内存放不下，若存在数据库，频繁访问数据库也需要时间。
  • AP算法的时间复杂度较高，一次迭代大概O(N3)
  • 聚类的好坏受到参考度和阻尼系数的影响

2.4 适用场景

• 许多簇，不均匀的簇大小，非平面几何
• 不可扩展的 n_samples
• 特别适合高维、多类数据快速聚类
• 图像、文本、生物信息学、人脸识别、基因发现、搜索最优航线、 码书设计以及实物图像识别等领域

代码

# -*- coding:utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import sklearn.datasets as ds
import matplotlib.colors
from sklearn.cluster import AffinityPropagation
from sklearn.metrics import euclidean_distances

#聚类算法之AP算法:
#1--与其他聚类算法不同,AP聚类不需要指定K(金典的K-Means)或者是其他描述聚类个数的参数
#2--一个聚类中最具代表性的点在AP算法中叫做Examplar,与其他算法中的聚类中心不同,examplar
#是原始数据中确切存在的一个数据点,而不是由多个数据点求平均而得到的聚类中心
#3--多次执行AP聚类算法,得到的结果完全一样的，即不需要进行随机选取初值步骤.
#算法复杂度较高,为O(N*N*logN),而K-Means只是O(N*K)的复杂度，当N》3000时,需要算很久
#AP算法相对于Kmeans优势是不需要指定聚类数量,对初始值不敏感

#AP算法应用场景：图像、文本、生物信息学、人脸识别、基因发现、搜索最优航线、 码书设计以及实物图像识别等领域

#算法详解: http://blog.csdn.net/helloeveryon/article/details/51259459

if __name__=='__main__':
    #scikit中的make_blobs方法常被用来生成聚类算法的测试数据，直观地说，make_blobs会根据用户指定的特征数量、
    # 中心点数量、范围等来生成几类数据，这些数据可用于测试聚类算法的效果。
    #函数原型：sklearn.datasets.make_blobs(n_samples=100, n_features=2,
    # centers=3, cluster_std=1.0, center_box=(-10.0, 10.0), shuffle=True, random_state=None)[source]
    #参数解析：
    # n_samples是待生成的样本的总数。
    #
    # n_features是每个样本的特征数。
    #
    # centers表示类别数。
    #
    # cluster_std表示每个类别的方差，例如我们希望生成2类数据，其中一类比另一类具有更大的方差，可以将cluster_std设置为[1.0, 3.0]。

    N=400
    centers = [[1, 2], [-1, -1], [1, -1], [-1, 1]]
    #生成聚类算法的测试数据
    data,y=ds.make_blobs(N,n_features=2,centers=centers,cluster_std=[0.5, 0.25, 0.7, 0.5],random_state=0)
    #计算向量之间的距离
    m=euclidean_distances(data,squared=True)
    #求中位数
    preference=-np.median(m)
    print 'Preference:',preference

    matplotlib.rcParams['font.sans-serif'] = [u'SimHei']
    matplotlib.rcParams['axes.unicode_minus'] = False
    plt.figure(figsize=(12,9),facecolor='w')
    for i,mul in enumerate(np.linspace(1,4,9)):#遍历等差数列
        print 'mul=',mul
        p=mul*preference
        model=AffinityPropagation(affinity='euclidean',preference=p)
        af=model.fit(data)
        center_indices=af.cluster_centers_indices_
        n_clusters=len(center_indices)
        print ('p=%.1f'%mul),p,'聚类簇的个数为:',n_clusters
        y_hat=af.labels_

        plt.subplot(3,3,i+1)
        plt.title(u'Preference：%.2f，簇个数：%d' % (p, n_clusters))
        clrs=[]
        for c in np.linspace(16711680, 255, n_clusters):
            clrs.append('#%06x' % c)
            for k, clr in enumerate(clrs):
                cur = (y_hat == k)
                plt.scatter(data[cur, 0], data[cur, 1], c=clr, edgecolors='none')
                center = data[center_indices[k]]
                for x in data[cur]:
                    plt.plot([x[0], center[0]], [x[1], center[1]], color=clr, zorder=1)
            plt.scatter(data[center_indices, 0], data[center_indices, 1], s=100, c=clrs, marker='*', edgecolors='k',
                        zorder=2)
            plt.grid(True)
        plt.tight_layout()
        plt.suptitle(u'AP聚类', fontsize=20)
        plt.subplots_adjust(top=0.92)
        plt.show()

6 Mean-shift 均值迁移

1)概述

• Mean-shift（即：均值迁移）的基本思想：在数据集中选定一个点，然后以这个点为圆心，r为半径，画一个圆(二维下是圆)，求出这个点到所有点的向量的平均值，而圆心与向量均值的和为新的圆心，然后迭代此过程，直到满足一点的条件结束。(Fukunage在1975年提出)
• 后来Yizong Cheng 在此基础上加入了 核函数 和 权重系数 ，使得Mean-shift 算法开始流行起来。目前它在聚类、图像平滑、分割、跟踪等方面有着广泛的应用。

2）图解过程

为了方便大家理解，借用下几张图来说明Mean-shift的基本过程。

由上图可以很容易看到，Mean-shift 算法的核心思想就是不断的寻找新的圆心坐标，直到密度最大的区域。

3）Mean-shift 算法函数

a）核心函数：sklearn.cluster.MeanShift(核函数：RBF核函数) 由上图可知，圆心(或种子)的确定和半径(或带宽)的选择，是影响算法效率的两个主要因素。所以在sklearn.cluster.MeanShift中重点说明了这两个参数的设定问题。

b）主要参数

• bandwidth ：半径(或带宽)，float型。如果没有给出，则使用sklearn.cluster.estimate_bandwidth计算出半径(带宽).（可选）
• seeds :圆心（或种子），数组类型，即初始化的圆心。（可选）
• bin_seeding ：布尔值。如果为真，初始内核位置不是所有点的位置，而是点的离散版本的位置，其中点被分类到其粗糙度对应于带宽的网格上。将此选项设置为True将加速算法，因为较少的种子将被初始化。默认值：False.如果种子参数(seeds)不为None则忽略。

c）主要属性

• cluster_centers_ : 数组类型。计算出的聚类中心的坐标。

• labels_ :数组类型。每个数据点的分类标签。、

• 算法特点
  • Distances between points（点之间的距离）
  • 圆心(或种子)的确定和半径(或带宽)的选择，是影响算法效率的两个主要因素。
  • 该算法不是高度可扩展的，因为在执行算法期间需要执行多个最近邻搜索。
  • 该算法保证收敛，但是当质心的变化较小时，算法将停止迭代。
  • 通过找到给定样本的最近质心来给新样本打上标签。
• 适用场景
  • 许多簇，不均匀的簇大小，非平面几何
  • 不可扩展的 n_samples
  • 适用于类别数量未知，且样本数量<10K
  • 目前它在聚类、图像平滑、分割、跟踪等方面有着广泛的应用。

代码

print(__doc__)

import numpy as np
from sklearn.cluster import MeanShift, estimate_bandwidth
from sklearn.datasets.samples_generator import make_blobs

# #############################################################################
# Generate sample data
centers = [[1, 1], [-1, -1], [1, -1]]
X, _ = make_blobs(n_samples=10000, centers=centers, cluster_std=0.6)

# #############################################################################
# Compute clustering with MeanShift

# The following bandwidth can be automatically detected using
bandwidth = estimate_bandwidth(X, quantile=0.2, n_samples=500)

ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X)
labels = ms.labels_
cluster_centers = ms.cluster_centers_

labels_unique = np.unique(labels)
n_clusters_ = len(labels_unique)

print("number of estimated clusters : %d" % n_clusters_)

# #############################################################################
# Plot result
import matplotlib.pyplot as plt
from itertools import cycle

plt.figure(1)
plt.clf()

colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
for k, col in zip(range(n_clusters_), colors):
    my_members = labels == k
    cluster_center = cluster_centers[k]
    plt.plot(X[my_members, 0], X[my_members, 1], col + '.')
    plt.plot(cluster_center[0], cluster_center[1], 'o', markerfacecolor=col,
             markeredgecolor='k', markersize=14)
plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()

7 BIRCH

========

• BIRCH

简介

• Birch(利用层次方法的平衡迭代规约和聚类)：就是通过聚类特征(CF)形成一个聚类特征树，root层的CF个数就是聚类个数。

算法原理

• 相关概念
  • 聚类特征(CF)：每一个CF是一个三元组,可以用（N，LS，SS）表示.其中N代表了这个CF中拥有的样本点的数量;LS代表了这个CF中拥有的样本点各特征维度的和向量,SS代表了这个CF中拥有的样本点各特征维度的平方和。
  • 如上图所示：N = 5
  • LS=(3+2+4+4+3,4+6+5+7+8)=(16,30)
  • SS =(32+22+42+42+32,42+62+52+72+82)=(54,190)
  • 
• 聚类过程
  • 对于上图中的CF Tree,限定了B=7,L=5， 也就是说内部节点最多有7个CF(CF90下的圆),而叶子节点最多有5个CF(CF90到CF94)。叶子节点是通过双向链表连通的。
  • 
• 特点
  • Euclidean distance between points（点之间的欧式距离）
• 适用场景
  • 大数据集，异常值去除，数据简化
  • 大的 n_clusters 和 n_samples
• 完整版代码

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import Birch

# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本2个特征，共4个簇，簇中心在[-1,-1], [0,0],[1,1], [2,2]
X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]], cluster_std=[0.4, 0.3, 0.4, 0.3],
                  random_state =9)

##设置birch函数
birch = Birch(n_clusters = None)
##训练数据
y_pred = birch.fit_predict(X)
##绘图
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

8 GaussianMixtureModel(混合高斯模型，GMM)

• 正太分布也叫高斯分布，正太分布的概率密度曲线也叫高斯分布概率曲线

1）介绍

• 聚类算法大多数通过相似度来判断，而相似度又大多采用 欧式距离 长短作为衡量依据。
• 而GMM采用了新的判断依据： 概率 ，即通过属于某一类的概率大小来判断最终的归属类别。
• GMM的基本思想就是： 任意形状的概率分布都可以用多个高斯分布函数去近似 ，也就是说GMM就是有多个单高斯密度分布（Gaussian）组成的，每个Gaussian叫一个”Component”，这些”Component”线性加成在一起就组成了 GMM 的概率密度函数，也就是下面的函数。

2）原理

• 高斯混合模型是由K个高斯分布（正态分布）函数组成，而该算法的目的就是找出各个高斯分布最佳的均值、方差、权重。
• 
  • 指定K的值，并初始随机选择各参数的值
  • E步骤。根据当前的参数，计算每个点由某个分模型生成的概率
  • M步骤。根据E步骤估计出的概率，来改进每个分模型的均值、方差和权重
  • 重复步骤2、3，直至收敛。

列出来公式只是方便理解下面的函数中为什么需要那些参数。 K：模型的个数，即Component的个数（聚类的个数） 为第k个高斯的权重

p（x |k） 则为第k个高斯概率密度,其均值为μk,方差为σk 上述参数，除了K是直接给定之外，其他参数都是通过EM算法估算出来的。(有个参数是指定EM算法参数的)

3）GaussianMixtureModel 算法函数

a）from sklearn.mixture.GaussianMixture b）主要参数（ 详细参数 ）

• n_components ：高斯模型的个数，即聚类的目标个数
• covariance_type : 通过EM算法估算参数时使用的协方差类型，默认是”full”
• full：每个模型使用自己的一般协方差矩阵
• tied：所用模型共享一个一般协方差矩阵
• diag：每个模型使用自己的对角线协方差矩阵

• spherical：每个模型使用自己的单一方差

• 算法特点
  • Mahalanobis distances to centers（Mahalanobis 与中心的距离）
• 优点
  • 可以给出一个样本属于某类的概率是多少
  • 不仅用于聚类，还可用于概率密度的估计
  • 可以用于生成新的样本点
• 缺点
  • 需要指定K值
  • 使用EM算法来求解
  • 往往只能收敛于局部最优
• 适用场景
  • 平面几何，适用于密度估计
  • Not scalable（不可扩展）
• 代码

import matplotlib as mpl
import matplotlib.pyplot as plt

import numpy as np

from sklearn import datasets
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import StratifiedKFold

print(__doc__)

colors = ['navy', 'turquoise', 'darkorange']


def make_ellipses(gmm, ax):
    for n, color in enumerate(colors):
        if gmm.covariance_type == 'full':
            covariances = gmm.covariances_[n][:2, :2]
        elif gmm.covariance_type == 'tied':
            covariances = gmm.covariances_[:2, :2]
        elif gmm.covariance_type == 'diag':
            covariances = np.diag(gmm.covariances_[n][:2])
        elif gmm.covariance_type == 'spherical':
            covariances = np.eye(gmm.means_.shape[1]) * gmm.covariances_[n]
        v, w = np.linalg.eigh(covariances)
        u = w[0] / np.linalg.norm(w[0])
        angle = np.arctan2(u[1], u[0])
        angle = 180 * angle / np.pi  # convert to degrees
        v = 2. * np.sqrt(2.) * np.sqrt(v)
        ell = mpl.patches.Ellipse(gmm.means_[n, :2], v[0], v[1],
                                  180 + angle, color=color)
        ell.set_clip_box(ax.bbox)
        ell.set_alpha(0.5)
        ax.add_artist(ell)
        ax.set_aspect('equal', 'datalim')

iris = datasets.load_iris()

# Break up the dataset into non-overlapping training (75%) and testing
# (25%) sets.
skf = StratifiedKFold(n_splits=4)
# Only take the first fold.
train_index, test_index = next(iter(skf.split(iris.data, iris.target)))


X_train = iris.data[train_index]
y_train = iris.target[train_index]
X_test = iris.data[test_index]
y_test = iris.target[test_index]

n_classes = len(np.unique(y_train))

# Try GMMs using different types of covariances.
estimators = {cov_type: GaussianMixture(n_components=n_classes,
              covariance_type=cov_type, max_iter=20, random_state=0)
              for cov_type in ['spherical', 'diag', 'tied', 'full']}

n_estimators = len(estimators)

plt.figure(figsize=(3 * n_estimators // 2, 6))
plt.subplots_adjust(bottom=.01, top=0.95, hspace=.15, wspace=.05,
                    left=.01, right=.99)


for index, (name, estimator) in enumerate(estimators.items()):
    # Since we have class labels for the training data, we can
    # initialize the GMM parameters in a supervised manner.
    estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)
                                    for i in range(n_classes)])

    # Train the other parameters using the EM algorithm.
    estimator.fit(X_train)

    h = plt.subplot(2, n_estimators // 2, index + 1)
    make_ellipses(estimator, h)

    for n, color in enumerate(colors):
        data = iris.data[iris.target == n]
        plt.scatter(data[:, 0], data[:, 1], s=0.8, color=color,
                    label=iris.target_names[n])
    # Plot the test data with crosses
    for n, color in enumerate(colors):
        data = X_test[y_test == n]
        plt.scatter(data[:, 0], data[:, 1], marker='x', color=color)

    y_train_pred = estimator.predict(X_train)
    train_accuracy = np.mean(y_train_pred.ravel() == y_train.ravel()) * 100
    plt.text(0.05, 0.9, 'Train accuracy: %.1f' % train_accuracy,
             transform=h.transAxes)

    y_test_pred = estimator.predict(X_test)
    test_accuracy = np.mean(y_test_pred.ravel() == y_test.ravel()) * 100
    plt.text(0.05, 0.8, 'Test accuracy: %.1f' % test_accuracy,
             transform=h.transAxes)

    plt.xticks(())
    plt.yticks(())
    plt.title(name)

plt.legend(scatterpoints=1, loc='lower right', prop=dict(size=12))

plt.show()

资料

结束