ISODATA算法 python实现


前言

ISODATA经常被用来与Kmeans算法进行对比,其本质也是按照欧式距离来对样本进行分类,不同的是ISODATA可以根据一个大概的指定类别数去确定最终的聚类数(两者可能不同),而Kmeans指定聚类数是多少后,最终的聚类就一定是多少。


一、ISODATA的流程

本质上只有分裂和合并两个步骤加更新中心三个步骤。了解这个算法,核心需要解决下面的三个问题:

Question 1. 什么时候分裂?

现有的聚类数太少就进行分裂。你一开始指定100个聚类,现在只有2个,那就进行分裂。(大的分裂方向,还有细节见下面流程图)

Question 2. 什么时候合并?

现有的聚类数太多就进行分裂。你一开始指定100个聚类,现在上一次刚好裂成200个,那就进行合并。(大的合并方向,还有细节见下面流程图)

Question 3. 现在有的中心数不上不下怎么办?

如果是奇数次迭代,那就尝试去分裂吧(虽然最后不一定分裂了)
如果是偶数次迭代,那就尝试去合并吧(虽然最后不一定合并了)

1.流程图(这里按迭代的奇偶来判断分裂或者合并)

在这里插入图片描述

注意:

在流程图中,“合并”步骤并不一定执行了合并,只有满足在所有的中心中,存在一些中心的距离太近(这个距离低于了设定的阈值)才会真正的执行合并的操作,其余不执行合并的操作。而在分裂中,只有现有的中心数太少或者满足”类内的距离太大而且样本数太少“进行分裂的操作。其中类内的距离太大则表示了这个聚类太过于松散,再加上类的数量太少的话,才进行分裂。

分裂的细节:如何分裂?

计算需要分裂的这个簇在各个维度上的方差,如果最大的方差超过了特定的阈值,就在这个最大方差的维度上分裂成两个,其他维度的值保持不变。

比如现在有一个中心 (1, 3) , 对于属于这个中心的所有样本,我们计算其在第一个维度 (数值1的维度) 的方差,再计算其在第二个维度 (数值3的维度) 的方差。假设维度1计算的方差结果为 0.3,维度2计算的方差为1.5,预先设定的阈值为0.5;所以我们要在第二个维度上把中心分成2个:(1, 3 + 1.5 * k ), (1, 3 - 1.5 *k) ,其中k又是控制分裂远近的一个超参数,在代码中取0.5。由此,我们得到了新分裂的两个中心,并把原来的中心去掉。

合并的细节: 如何合并?

合并使用加权平均的方法,两个权重是两个中心控制的两簇样本的数量百分比,加权求和即可。

二、使用步骤

1.代码实现

Tips: 注意需要用到sklearn的库来产生数据集:

# -*- encoding:utf-8 -*-
"""
@author:zsiming
@fileName:ISODATA.py
@Time:2022/1/9  12:33
"""
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.metrics import euclidean_distances


class ISODATA():
    def __init__(self, designCenterNum, LeastSampleNum, StdThred, LeastCenterDist, iterationNum):
        #  指定预期的聚类数、每类的最小样本数、标准差阈值、最小中心距离、迭代次数
        self.K = designCenterNum
        self.thetaN = LeastSampleNum
        self.thetaS = StdThred
        self.thetaC = LeastCenterDist
        self.iteration = iterationNum

        # 初始化
        self.n_samples = 1500
        # 选一
        self.random_state1 = 200
        self.random_state2 = 160
        self.random_state3 = 170
        self.data, self.label = make_blobs(n_samples=self.n_samples, random_state=self.random_state3)

        self.center = self.data[0, :].reshape((1, -1))
        self.centerNum = 1
        self.centerMeanDist = 0

        # seaborn风格
        sns.set()

    def updateLabel(self):
        """
            更新中心
        """
        for i in range(self.centerNum):
            # 计算样本到中心的距离
            distance = euclidean_distances(self.data, self.center.reshape((self.centerNum, -1)))
            # 为样本重新分配标签
            self.label = np.argmin(distance, 1)
            # 找出同一类样本
            index = np.argwhere(self.label == i).squeeze()
            sameClassSample = self.data[index, :]
            # 更新中心
            self.center[i, :] = np.mean(sameClassSample, 0)

        # 计算所有类到各自中心的平均距离之和
        for i in range(self.centerNum):
            # 找出同一类样本
            index = np.argwhere(self.label == i).squeeze()
            sameClassSample = self.data[index, :]
            # 计算样本到中心的距离
            distance = np.mean(euclidean_distances(sameClassSample, self.center[i, :].reshape((1, -1))))
            # 更新中心
            self.centerMeanDist += distance
        self.centerMeanDist /= self.centerNum

    def divide(self):
        # 临时保存更新后的中心集合,否则在删除和添加的过程中顺序会乱
        newCenterSet = self.center
        # 计算每个类的样本在每个维度的标准差
        for i in range(self.centerNum):
            # 找出同一类样本
            index = np.argwhere(self.label == i).squeeze()
            sameClassSample = self.data[index, :]
            # 计算样本到中心每个维度的标准差
            stdEachDim = np.mean((sameClassSample - self.center[i, :])**2, axis=0)
            # 找出其中维度的最大标准差
            maxIndex = np.argmax(stdEachDim)
            maxStd = stdEachDim[maxIndex]
            # 计算样本到本类中心的距离
            distance = np.mean(euclidean_distances(sameClassSample, self.center[i, :].reshape((1, -1))))
            # 如果最大标准差超过了阈值
            if maxStd > self.thetaS:
                # 还需要该类的样本数大于于阈值很多 且 太分散才进行分裂
                if self.centerNum <= self.K//2 or \
                        sameClassSample.shape[0] > 2 * (self.thetaN+1) and distance >= self.centerMeanDist:
                    newCenterFirst = self.center[i, :].copy()
                    newCenterSecond = self.center[i, :].copy()

                    newCenterFirst[maxIndex] += 0.5 * maxStd
                    newCenterSecond[maxIndex] -= 0.5 * maxStd

                    # 删除原始中心
                    newCenterSet = np.delete(newCenterSet, i, axis=0)
                    # 添加新中心
                    newCenterSet = np.vstack((newCenterSet, newCenterFirst))
                    newCenterSet = np.vstack((newCenterSet, newCenterSecond))

            else:
                continue
        # 更新中心集合
        self.center = newCenterSet
        self.centerNum = self.center.shape[0]

    def combine(self):
        # 临时保存更新后的中心集合,否则在删除和添加的过程中顺序会乱
        delIndexList = []

        # 计算中心之间的距离
        centerDist = euclidean_distances(self.center, self.center)
        centerDist += (np.eye(self.centerNum)) * 10**10
        # 把中心距离小于阈值的中心对找出来
        while True:
            # 如果最小的中心距离都大于阈值的话,则不再进行合并
            minDist = np.min(centerDist)
            if minDist >= self.thetaC:
                break
            # 否则合并
            index = np.argmin(centerDist)
            row = index // self.centerNum
            col = index % self.centerNum
            # 找出合并的两个类别
            index = np.argwhere(self.label == row).squeeze()
            classNumFirst = len(index)
            index = np.argwhere(self.label == col).squeeze()
            classNumSecond = len(index)
            newCenter = self.center[row, :] * (classNumFirst / (classNumFirst+ classNumSecond)) + \
                        self.center[col, :] * (classNumSecond / (classNumFirst+ classNumSecond))
            # 记录被合并的中心
            delIndexList.append(row)
            delIndexList.append(col)
            # 增加合并后的中心
            self.center = np.vstack((self.center, newCenter))
            self.centerNum -= 1
            # 标记,以防下次选中
            centerDist[row, :] = float("inf")
            centerDist[col, :] = float("inf")
            centerDist[:, col] = float("inf")
            centerDist[:, row] = float("inf")

        # 更新中心
        self.center = np.delete(self.center, delIndexList, axis=0)
        self.centerNum = self.center.shape[0]

    def drawResult(self):
        ax = plt.gca()
        ax.clear()
        ax.scatter(self.data[:, 0], self.data[:, 1], c=self.label, cmap="cool")
        # ax.set_aspect(1)
        # 坐标信息
        ax.set_xlabel('x axis')
        ax.set_ylabel('y axis')
        plt.show()


    def train(self):
        # 初始化中心和label
        self.updateLabel()
        self.drawResult()

        # 到设定的次数自动退出
        for i in range(self.iteration):
            # 如果是偶数次迭代或者中心的数量太多,那么进行合并
            if self.centerNum < self.K //2:
                self.divide()
            elif (i > 0 and i % 2 == 0) or self.centerNum > 2 * self.K:
                self.combine()
            else:
                self.divide()
            # 更新中心
            self.updateLabel()
            self.drawResult()
            print("中心数量:{}".format(self.centerNum))



if __name__ == "__main__":
    isoData = ISODATA(designCenterNum=5, LeastSampleNum=20, StdThred=0.1, LeastCenterDist=2, iterationNum=20)
    isoData.train()


2.迭代过程

1. 原始数据如下图所示,可以看见我在这儿比较明显的生成三个簇的数据(然后指定类别数为5):

2. 从一个中心分裂成为两个中心(用颜色区分不同的聚类):
在这里插入图片描述

3. 未到达指定类别数(2 < 5)继续分裂为4个中心:

在这里插入图片描述

4.中心贴得太近了,需要合并:
在这里插入图片描述

5. 更新中心的位置和分裂:
在这里插入图片描述

6.中心贴得太近了,合并
在这里插入图片描述

7.后面将不再变化。
在这里插入图片描述


3. 总结

个人觉得:

从参数的角度来看,相比于Kmeans,由一个超参数数变成了六个超参数,不能说是改进。只能说某些先验知识比较完善的情况下,可能适用于数据流形分布比较复杂的场景。

  • 19
    点赞
  • 62
    收藏
    觉得还不错? 一键收藏
  • 19
    评论
### 回答1: ISODATA算法(Iterative Self-Organizing Data Analysis Technique Algorithm)是一种聚类算法,可以用来将无标签数据划分到不同的群组中。在Python中,可以使用Scikit-learn库中的函数实现ISODATA算法。例如,可以使用sklearn.cluster.IsoClustering函数进行聚类操作。 ### 回答2: ISODATA算法是一种常用的聚类分析方法,被广泛应用于数据分析、图像处理、模式识别和数据挖掘等领域。Python是目前最流行的编程语言之一,拥有丰富的科学计算、数据处理和机器学习库,方便快速地实现ISODATA算法ISODATA算法的主要步骤包括初始化、计算距离矩阵、聚类、计算质心、合并和分裂等部分。以下为各部分的详细介绍: 1. 初始化:设定聚类数、最小样本个数、最大样本个数、最大迭代次数、误差容限和质心的初始值等参数。 2. 计算距离矩阵:通过样本点之间的距离计算得到距离矩阵,可以使用scipy.spatial.distance库中的pdist函数实现。 3. 聚类:将样本点归为最近的聚类中心所属簇,可以使用scikit-learn库中的KMeans函数实现。 4. 计算质心:重新计算每个簇的质心。 5. 合并和分裂:根据簇内样本数、质心距离和误差容限等条件来决定是否将簇合并或者簇内部分裂。 6. 判断终止条件:当迭代次数达到最大迭代次数、簇的个数小于等于最小样本个数或者簇的个数大于等于最大样本个数时,算法停止迭代。 下面是一个简单的ISODATA算法Python实现: ``` import numpy as np from scipy.spatial.distance import pdist from sklearn.cluster import KMeans def isodata(data, k, min_samples, max_samples, max_iter, tol, var): n_samples, n_features = data.shape centers = data[np.random.choice(n_samples, k)] labels = np.zeros(n_samples) idx = np.array(range(n_samples)) it = 0 while (len(set(labels)) > 1) and (it < max_iter) and (k <= max_samples): dist_mat = pdist(data) dist_mat = np.reshape(dist_mat, (n_samples, n_samples)) # cluster assignments for i in idx: center_dists = np.sum((centers - data[i])**2, axis=1) labels[i] = np.argmin(center_dists) # re-compute centers for i in range(k): cluster_i = data[labels == i] if len(cluster_i) > 0: centers[i] = np.mean(cluster_i, axis=0) # split/merge clusters k_orig = k for i in range(k_orig): cluster_i = data[labels == i] if len(cluster_i) <= min_samples: labels[labels == i] = -1 k -= 1 elif len(cluster_i) >= (var * n_samples): centers = np.vstack((centers, centers[i])) labels[labels == i] = k labels[labels == -1] = i k += 1 # check termination conditions if len(set(labels)) <= 1 or k >= max_samples: break it += 1 return labels ``` 通过上述代码,我们可以完成一个简单的ISODATA算法Python实现。其中,使用了numpy、scipy和scikit-learn等库来实现各部分功能,具体使用时需要针对具体应用场景进行参数调整和算法优化,以获得更好的聚类效果。 ### 回答3: ISODATA算法是一种基于聚类的图像分割算法,主要用于将一幅图像分割成多个区域,每个区域具有相同的颜色和亮度。这种算法通过不断地计算和调整区域内像素的均值和标准差来确定区域的数量和大小。 在Python实现ISODATA算法,首先需要确定参数,如初始类别数量、最小误差、最小类别大小等。接着,需要将图像上的像素根据其像素值聚类成初始类别,并计算每个类别的均值和标准差。然后,根据计算出来的均值和标准差,将相邻的类别合并,直到满足停止条件为止。最后,将每个类别内的像素值赋予相应的颜色,得到分割后的图像。 这里给出一个简单的实现过程: ``` import cv2 def ISODATA(image, K, min_error, min_size): # 将图像转为灰度图 gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) # 初始化初始类别数量 num_classes = K # 初始化每个类别的均值和标准差 means = [] stds = [] for i in range(num_classes): means.append(gray[i * gray.shape[0] // num_classes:(i+1) * gray.shape[0] // num_classes, :].mean()) stds.append(gray[i * gray.shape[0] // num_classes:(i+1) * gray.shape[0] // num_classes, :].std()) # 循环直到满足停止条件 while num_classes > 1: # 合并相邻的类别 merge_idx = None min_dist = float("inf") for i in range(num_classes - 1): dist = abs(means[i+1] - means[i]) if dist < min_dist: min_dist = dist merge_idx = i if min_dist < min_error: # 计算新的类别均值和标准差 new_means = means[:merge_idx] + [(means[merge_idx] + means[merge_idx+1])/2] + means[merge_idx+2:] new_stds = stds[:merge_idx] + [(stds[merge_idx] + stds[merge_idx+1])/2] + stds[merge_idx+2:] # 判断每个类别是否满足最小类别大小条件 mask = cv2.inRange(gray, merge_idx * 255 / num_classes, (merge_idx+2) * 255 / num_classes) if mask.sum() < min_size: break # 更新参数 means = new_means stds = new_stds num_classes -= 1 else: break # 将每个类别内的像素值赋予相应的颜色 output = cv2.cvtColor(gray, cv2.COLOR_GRAY2BGR) for i in range(num_classes): mask = cv2.inRange(gray, i * 255 / num_classes, (i+1) * 255 / num_classes) output[mask > 0] = (255*i/num_classes, 0, 0) return output ``` 以上代码仅为简要示例,实际应用中需要对参数进行调整,并根据实际需求进行算法的优化。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 19
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值