k-means and FCM


原理简介

K-means

  1. 指定需要划分的簇的个数,即K值(类的个数)
  2. 随机的选择K个数据对象作为初始的聚类中心
  3. 计算其余的各个数据对象到这K个聚类中心的距离,把数据对象划分到距离它最近的那个中心所在的簇中
  4. 调整新类并计算出新的聚类中心
  5. 循环第三和第四步,直到聚类中心不再发生变化,聚类结束

FCM

  1. 指定需要划分的类别个数以及柔性参数的值,即c和$\alpha$

  2. 初始化隶属度矩阵,其中样本数据属于所有类别的隶属度之和应该为1

    $$
    \sum_{j=1}^{c}u_{ij}=1,\forall{i}=1,…,m
    $$
    其中$i$对应样本,$j$对应聚类中心,$u_{ij}$即为隶属度,它是一个位于区间$[0,1]$之间的数,它表示样本$x_i$属于第$j$类的概率为$u_{ij}$

  3. 根据隶属度矩阵计算新的c个聚类中心,计算公式如下:

    $$
    z_j=\frac{\sum_{i=1}^{m}u_{ij}^\alpha{x_i}}{\sum_{i=1}^{m}u_{ij}^\alpha}
    $$
    其中$x_i$表示一个样本,$z_j$表示一个聚类中心,它是一个向量

  4. 聚类中心不再变化为终止条件,判断此次迭代是否满足终止条件

  5. 如果满足终止条件,那么对于每一个样本,将它归类于隶属度最大的那一类

  6. 如果不满足终止条件,就更新隶属度矩阵$U$,并回到第3步继续迭代,更新公式如下:

    $$
    u_{ij}=\frac{1}{\sum_{k=1}^{c}(\frac{d_{ij}}{d_{ki}})^{\frac{2}{\alpha-1}}}
    $$
    其中$d_{ij}$表示样本$x_i$到聚类中心点$z_j$的欧式距离

评价指标

由于聚类属于无监督学习,用OA、AA等指标来评价聚类效果可能得不到理想的结果。故在本次作业中,使用轮廓系数作为聚类的评价指标。

轮廓系数取值范围为$[-1,1]$,取值越接近1则说明聚类性能越好,相反,取值越接近-1则说明聚类性能越差,一个样本的轮廓系数计算公式如下:

$$
s_i=\frac{b_i-a_i}{\max{a_i,b_i}}
$$
其中:

  • $a_i$:样本$x_i$与其所在簇内其他样本的平均距离
  • $b_i$:样本$x_i$与其他簇内样本的平均距离

当$b_i$大于$a_i$时,说明样本$x_i$跟同类的样本之间的平均距离小于其与其他类别样本之间的平均距离,令$\Delta_1=b_i-a_i$,$\Delta_1$越大,$s_i$就越接近1,此时样本类间的紧密程度较高,类间的稀疏程度较低,聚类效果就较为理想

同理当$b_i$小于$a_i$时,说明样本$x_i$跟同类的样本之间的平均距离大于其与其他类别样本之间的平均距离,令$\Delta_2=a_i-b_i$,$\Delta_2$越大,$s_i$就越接近-1,此时样本类间的紧密程度较低,类间的紧密程度甚至高于类内紧密程度,这时聚类效果就比较糟糕

对于$m$个样本,样本整体的轮廓系数即为各个样本轮廓系数的平均值:

$$
S=\frac{1}{m}\sum_{i=1}^{m}s_i
$$

数据标准化

由于两数据集中特征的量纲都相同,所以只需要消除量纲,将每个特征的数值除以该特征中的最大值即可,这样就将所有特征数值映射到了区间$[0,1]$内

类数的选择

因为我们知道sonar数据集和iris数据集对应的类别数分别为2、3。所以在聚类的过程中将k值和c值取成对应的类别即可

编程思路

K-means的实现

导入数据部分

iris数据集

直接从sklearn包中导入即可:

from sklearn.datasets import load_iris

# 读取数据
data=load_iris()
# 标签
target=data.target.T
# 样本数据
data=data.data
sonar数据集

从本地的sonar.xlsx导入:

from openpyxl import load_workbook

# 导入数据
def get_data(path):
    """
    传入文件路径,返回样本数据以及标签
    """
    workbook=load_workbook(path)
    # print(workbook.sheetnames)
    sheet=workbook['Sheet1']
    # 存储样本特征集
    sonar=np.zeros([208,60])
    # 存储标签
    target=np.zeros([208,1])
    for i in range(208):
        for j in range(60):
            sonar[i,j]=sheet.cell(row=i+1,column=j+1).value
    # print(sonar.shape)
    # print(sonar[0,59])
    for i in range(208):
        if sheet.cell(row=i+1,column=61).value=='R':
            target[i]=0
        else:
            target[i]=1
    return sonar,target

随机选择聚类中心

产生k个随机索引即可

# 指定k-means中的k值
    k=2
    # 随机选择k个聚类中心
    flags= random.sample(range(0,data.shape[0]),k)
    # print(flags)
    center_points=np.zeros([k,data.shape[1]])
    j=0
    for i in range(k):
        center_points[[j],:]=data[flags[i],:]
        j+=1

算法核心的实现

编写一个类,名称为k-means,它的属性以及方法如下:

利用python实现的代码如下:

## 创建一个k-means类
class K_means():
    def __init__(self,data,center_points, k):
        self.center=center_points #聚类中心点的坐标
        self.k=k #聚类数
        self.data=data #样本
        # dis前k列表示每个样本到各聚类中心的距离,最后一列表示类别
        self.dis=np.zeros([self.data.shape[0],self.k+1])


    # distance方法:计算当前到k个聚类中心的距离,找出最小距离并将每个样本分类
    def distance(self):
        for i in range(len(self.data)):
            for j in range(self.k):
                ## 距离取为欧氏距离,即为二范数
                self.dis[i,j]=np.linalg.norm(self.\
                    center[j,:]-self.data[i],ord=2)
            self.dis[i,self.k]=np.where(self.dis[i,:]\
               ==np.min(self.dis[i,0:self.k]))[0][0]


    ## 根据self.dis更新聚类中心self.center
    def update(self):
        for i in range(self.k):
            #设置统计数组每次循环将同类的样本存储起来
            count=np.zeros([1,self.data.shape[1]])
            # 逐类进行统计并计算新的聚类中心
            for j in range(len(self.data)):
                if self.dis[j,self.k]==i:
                    # 一行一行进行拼接
                    count=np.row_stack((count,self.data[j,:]))
                else:
                    continue
            # 更新聚类中心
            # print(count.shape)
            self.center[[i],:]= np.mean(count,axis=0)
        # print(self.center)

迭代更新部分

循环调用k-means类里的方法即可达到迭代更新的目的,终止条件:当此次的聚类中心点和上次的聚类中心点相比变化小于设定阈值时停止迭代

# 实例化一个K_means类
iris= K_means(data, target,center_points, k)
# print(iris.dis)
    
# center用来和iris.center比较,相同则迭代终止
center=np.zeros_like(iris.center)
    
# 迭代部分
i=0
while np.linalg.norm(iris.center-center)>0.0000001:
    # 刚用class写,踩到的一个坑,这里要调用copy方法,
    # 直接等于表示只修改指针,不可取
    center=iris.center.copy()
    iris.distance()
    # print(center)
    iris.update()
    i+=1
    print("迭代第{}次".format(i))
    
print("算法收敛,迭代终止")

评价聚类效果

利用sklearn包中的函数计算轮廓系数

from sklearn import metrics

## 计算轮廓系数来评价聚类效果
print(metrics.silhouette_score(data,\
    iris.dis[:,-1], metric='euclidean'))

FCM的实现

导入数据部分

与K-means相同

随机初始化隶属度矩阵

![class:FCM](E:%5Cblog_pic%5Cclass%EF%BC%9AFCM.jpg)# 指定FCM中的c值
c = 3
# 随机初始化隶属度矩阵
U = np.zeros([data.shape[0], c])
for i in range(data.shape[0]):
    for j in range(c):
        if np.sum(U[i, :]) < 1 and j < c - 1:
            U[i, j] = random.uniform(0, 1 - np.sum(U[i, :]))
    # 最后一列单独赋值
    U[i, c - 1] = 1 - np.sum(U[i, :])

算法核心的实现

编写一个类,名称为FCM,它的属性以及方法如下:

利用python实现的代码如下:

class FCM():
    def __init__(self, c, alpha, data, U):
        # 聚类数目
        self.c = c
        # 柔性参数(衡量模糊度)
        self.alpha = alpha
        # 样本数据
        self.data = data
        # 隶属度矩阵
        self.U = U
        # 聚类中心
        self.center = np.zeros([self.c, self.data.shape[1]])
        # 代价函数
        self.J = 0
        # 存储距离
        self.dis = np.zeros([self.data.shape[0], self.c])
        # 聚类结果
        self.target = np.zeros([self.data.shape[0], 1])


    # 更新聚类中心
    def update_center(self):
        for i in range(self.c):
            self.center[[i], :] = np.sum(np.power(self.U[:,\
                [i]],self.alpha)*self.data, axis=0) / np.sum\
                    (np.power(self.U[:, [i]], self.alpha))

    # 更新隶属度矩阵
    def update_U(self):
        # 计算距离,并存入self.dis中
        for i in range(self.dis.shape[0]):
            for j in range(self.dis.shape[1]):
                self.dis[i, j] = np.linalg.norm\
                    (self.data[[i], :] - self.center[[j], :])
        # print(self.dis)
        # 根据公式计算新的隶属度矩阵
        for i in range(self.U.shape[0]):
            for j in range(self.U.shape[1]):
                self.U[i, j] = 1 / np.sum(np.power(self.dis\
                    [i, j] / self.dis[[i], :], 2-(self.alpha-1)))
        # print(self.U)

    # 计算代价函数
    def compute_J(self):
        self.J = 0
        for j in range(self.c):
            J = 0
            for i in range(self.data.shape[0]):
                J += np.power(self.U[i, j], self.alpha) * \
                    np.power(self.dis[i, j], 2)
            # print(J)
            self.J += J

    # 根据当前隶属度矩阵判断类别
    def compute_target(self):
        for i in range(self.data.shape[0]):
            self.target[i, 0] = np.argmax(self.U[[i], :])

迭代更新部分

循环调用FCM类里的方法即可达到迭代更新的目的,终止条件:当此次的聚类中心点和上次的聚类中心点相比变化小于设定阈值时停止迭代

iris = FCM(c, alpha, data, U)
# print(sonar.U)
# 保留更新前的中心点
center = np.ones_like(iris.center)
    
# 迭代部分
i = 1
while np.linalg.norm(iris.center-center)>0.0000001:
    center=iris.center.copy()
    # print(U)
    iris.update_center()
    # print(iris.center)
    iris.compute_J()
    # print(iris.center)
    iris.update_U()
    print("迭代第{}次".format(i))
    iris.compute_target()
    i+=1
        
print("算法收敛,迭代终止")

评价聚类效果

与K-means相同

from sklearn import metrics

# 计算轮廓系数来评价聚类效果
print(metrics.silhouette_score(data,\
    iris.target.reshape(-1), metric='euclidean'))

结果展示与分析

结果展示

K-means

iris数据集

轮廓系数为0.519,聚类效果还可以

sonar数据集

轮廓系数为0.198,效果一般

FCM

iris数据集

轮廓系数为0.527,与k-means效果相当

sonar数据集

轮廓系数为0.198,与k-means效果相当

结果分析

可以看出,在性能方面,k-means算法和FCM算法基本上差距不大;但在收敛速度上,k-means算法要略优于FCM算法,收敛速度慢,容易陷入局部最小是FCM算法的缺点,改进可以使用IFCM算法(直觉模糊c均值聚类),它相较于FCM算法增加了新的参数——非隶属度$\gamma$和不确定度$\pi$,能更加细腻的刻画出客观世界的模糊性质。

拓展——k-means算法用于图像分割

原理简介

将一张RGB三通道的图像中的每一个像素点都视为一个样本,那么每个样本都对应有3个特征,分别为R、G、B的值,如果像素点属于同一区域,那么它们的特征值一定很相近,此时就可以利用聚类算法让特征值相近的像素点聚为一类,最后将不同类别的像素点用不同的颜色或者灰度值输出即可

样本展示

在此处取聚类的k值为3,意思是将图片分割为3类,分别是天空、山坡和水面。

编程思路

整体思路和前面的思路基本一致,只需要修改图像输入部分以及加上最后的结果图展示部分

图像输入

利用pillow包中的Imge库将图片中的每一个像素点的rgb值都提取出来,成为一个行向量,然后将所有行向量拼接起来,变成一个样本矩阵。

from PIL import Image

# 读取数据
img=Image.open('3.bmp','r')
lst=list(img.getdata())
# print(type(img.size))
data=np.array(lst)

结果图展示部分

根据聚类得到的标签值(一个列向量),先将它重组为一个和原图像尺寸一致的矩阵,再根据不同的类别给像素点赋上不同的RGB值,然后再根据RGB值输出图像即可

from PIL import Image

## 输出分割后的图像
def plot_pic(*args):
    k=args[0]
    rows,columns=args[1],args[2]
    targets=args[3]
    img=np.zeros([rows,columns])
    img=targets.reshape([columns,rows]).T
    im = Image.new("RGB",(rows,columns))#创建图片
    for i in range(rows):
        for j in range(columns):
            if img[i,j]==0:
                im.putpixel((i,j),(255,0,0))
            if img[i,j]==1:
                im.putpixel((i,j),(0,255,0))
            if img[i,j]==2:
                im.putpixel((i,j),(0,0,255))
    im.show()
    im.save('3-result.jpg')

结果展示

首先来看一下轮廓系数

轮廓系数为0.669,说明聚类的效果比较好,再来看看分割后的图像

可以看出,整体上图像分割的效果还不错,基本上将天空,山坡和水区分开来。

完整代码

k-means

iris

from sklearn.datasets import load_iris
import numpy as np
from sklearn import metrics
import random


## 创建一个k-means类
class K_means():
    def __init__(self,data, target,center_points, k):
        self.center=center_points #聚类中心点的坐标
        self.k=k #聚类数
        self.data=data #样本
        self.target=target #样本对应的标签
        # dis前k列表示每个样本到各聚类中心的距离,最后一列表示类别
        self.dis=np.zeros([self.data.shape[0],self.k+1])
    
    
    # distance方法:计算当前到k个聚类中心的距离,找出最小距离并将每个样本分类
    def distance(self): 
        for i in range(len(self.data)):
            for j in range(self.k):
                ## 距离取为欧氏距离,即为二范数
                self.dis[i,j]=np.linalg.norm(self.\
                    center[j,:]-self.data[i],ord=2)
            self.dis[i,self.k]=np.where(self.dis[i,:]\
                ==np.min(self.dis[i,0:self.k]))[0][0]
   
        
    ## 根据self.dis更新聚类中心self.center    
    def update(self):
        for i in range(self.k):
            #设置统计数组每次循环将同类的样本存储起来
            count=np.zeros([1,self.data.shape[1]]) 
            # 逐类进行统计并计算新的聚类中心
            for j in range(len(self.data)):
                if self.dis[j,self.k]==i:
                    # 一行一行进行拼接
                    count=np.row_stack((count,self.data[j,:])) 
                else:
                    continue
            # 更新聚类中心
            # print(count.shape)
            self.center[[i],:]= np.mean(count,axis=0)
        # print(self.center)
            

if __name__ == "__main__":
    
    # 读取数据
    data=load_iris()
    # 标签
    target=data.target.T
    # 样本数据
    data=data.data
    # 指定k-means中的k值
    k=3
    # 随机选择k个聚类中心
    flags= random.sample(range(0,data.shape[0]),k)
    # print(flags)
    center_points=np.zeros([k,data.shape[1]])
    j=0
    for i in range(k):
        center_points[[j],:]=data[flags[i],:]
        j+=1
    
    # 实例化一个K_means类
    iris= K_means(data, target,center_points, k)
    # print(iris.dis)
    
    # center用来和iris.center比较,相同则迭代终止
    center=np.zeros_like(iris.center)
    
    # 迭代部分
    i=0
    while np.linalg.norm(iris.center-center)>0.0000001:
    # while i<50:
        # 刚用class写,踩到的一个坑,这里要调用copy方法,
        # 直接等于表示只修改指针,不可取
        center=iris.center.copy()
        iris.distance()
        # print(center)
        iris.update()
        i+=1
        print("迭代第{}次".format(i))
        
    print("算法收敛,迭代终止")
    # print(iris.dis)
    # 计算轮廓系数来评价聚类效果
    print(metrics.silhouette_score(data, iris.dis[:,-1], metric='euclidean'))

sonar

from openpyxl import load_workbook
from sklearn import metrics
import numpy as np
import random


## 创建一个k-means类,编写一些方法
class K_means():
    def __init__(self,data, target,center_points, k):
        self.center=center_points #聚类中心点的坐标
        self.k=k #聚类数
        self.data=data #样本
        self.target=target #样本对应的标签
        # dis前k列表示每个样本到各聚类中心的距离,最后一列表示类别
        self.dis=np.zeros([self.data.shape[0],self.k+1])
    
    
    #distance方法:计算当前到k个聚类中心的距离,找出最小距离并将每个样本分类
    def distance(self): 
        for i in range(len(self.data)):
            for j in range(self.k):
                ## 距离取为欧氏距离,即为二范数
                self.dis[i,j]=np.linalg.norm\
                    (self.center[j,:]-self.data[i],ord=2)
            self.dis[i,self.k]=np.where(self.dis[i,:]\
                ==np.min(self.dis[i,0:self.k]))[0][0]
   
        
    ## 根据self.dis更新聚类中心self.center    
    def update(self):
        for i in range(self.k):
            #设置统计数组每次循环将同类的样本存储起来
            count=np.zeros([1,self.data.shape[1]]) 
            # 逐类进行统计并计算新的聚类中心
            for j in range(len(self.data)):
                if self.dis[j,self.k]==i:
                    # 一行一行进行拼接
                    count=np.row_stack((count,self.data[j,:])) 
                else:
                    continue
            # 更新聚类中心
            # print(count.shape)
            self.center[[i],:]= np.mean(count,axis=0)
        # print(self.center)
                       
# 导入数据
def get_data(path):
    """
    传入文件路径,返回样本数据以及标签
    """
    workbook=load_workbook(path)
    # print(workbook.sheetnames)
    sheet=workbook['Sheet1']
    # 存储样本特征集
    sonar=np.zeros([208,60])
    # 存储标签
    target=np.zeros([208,1])
    for i in range(208):
        for j in range(60):
            sonar[i,j]=sheet.cell(row=i+1,column=j+1).value
    # print(sonar.shape)
    # print(sonar[0,59])
    for i in range(208):
        if sheet.cell(row=i+1,column=61).value=='R':
            target[i]=0
        else:
            target[i]=1
    # print(target[97])
    # sonar1=sonar[0:97,:]
    # sonar2=sonar[97:208,:]
    # target1=target[0:97,:]
    # target2=target[97:208,:]  
    return sonar,target

if __name__ == "__main__":

    #读取数据
    path='sonar.xlsx'
    data,target=get_data(path)
    # 指定k-means中的k值
    k=2
    # 随机选择k个聚类中心
    flags= random.sample(range(0,data.shape[0]),k)
    # print(flags)
    center_points=np.zeros([k,data.shape[1]])
    j=0
    for i in range(k):
        center_points[[j],:]=data[flags[i],:]
        j+=1
    
    # 实例化一个K_means类
    sonar= K_means(data, target,center_points, k)
    # print(sonar.dis)
    
    # center用来和sonar.center比较,相同则迭代终止
    center=np.zeros_like(sonar.center)
    
    # 迭代部分
    i=0
    while np.linalg.norm(sonar.center-center)>0.00000001:
        # 刚用class写,踩到的一个坑,这里要调用copy方法,
        # 直接等于表示只修改指针,不可取
        center=sonar.center.copy()
        sonar.distance()
        # print(center)
        sonar.update()
        i+=1
        print("迭代第{}次".format(i))
        
    print("算法收敛,迭代终止")
    # print(sonar.dis)
    # 计算轮廓系数来评价聚类效果
    print(metrics.silhouette_score(data,\
        sonar.dis[:,-1], metric='euclidean'))

FCM

iris

from sklearn import metrics
from sklearn.datasets import load_iris
import numpy as np
import random


class FCM():
    def __init__(self, c, alpha, data, U):
        # 聚类数目
        self.c = c
        # 柔性参数(衡量模糊度)
        self.alpha = alpha
        # 样本数据
        self.data = data
        # 隶属度矩阵
        self.U = U
        # 聚类中心
        self.center = np.zeros([self.c, self.data.shape[1]])
        # 代价函数
        self.J = 0
        # 存储距离
        self.dis = np.zeros([self.data.shape[0], self.c])
        # 聚类结果
        self.target = np.zeros([self.data.shape[0], 1])


    # 更新聚类中心
    def update_center(self):
        for i in range(self.c):
            self.center[[i], :] = np.sum(np.power(self.U[:,\
                [i]],self.alpha)*self.data, axis=0) / np.sum\
                    (np.power(self.U[:, [i]], self.alpha))

    # 更新隶属度矩阵
    def update_U(self):
        # 计算距离,并存入self.dis中
        for i in range(self.dis.shape[0]):
            for j in range(self.dis.shape[1]):
                self.dis[i, j] = np.linalg.norm\
                    (self.data[[i], :] - self.center[[j], :])
        # print(self.dis)
        # 根据公式计算新的隶属度矩阵
        for i in range(self.U.shape[0]):
            for j in range(self.U.shape[1]):
                self.U[i, j] = 1 / np.sum(np.power(self.dis\
                    [i, j] / self.dis[[i], :], 2-(self.alpha-1)))
        # print(self.U)

    # 计算代价函数
    def compute_J(self):
        self.J = 0
        for j in range(self.c):
            J = 0
            for i in range(self.data.shape[0]):
                J += np.power(self.U[i, j], self.alpha) * \
                    np.power(self.dis[i, j], 2)
            # print(J)
            self.J += J

    # 根据当前隶属度矩阵判断类别
    def compute_target(self):
        for i in range(self.data.shape[0]):
            self.target[i, 0] = np.argmax(self.U[[i], :])

         
            
if __name__ == "__main__":
    
    # 读取数据
    data=load_iris()
    # 标签
    target=data.target.T
    # 样本数据
    data=data.data
    # 指定FCM中的c值
    c = 3
    # 设置柔性参数alpha
    alpha = 2
    # 随机生成隶属度矩阵U
    # 随机初始化隶属度矩阵
    U = np.zeros([data.shape[0], c])
    for i in range(data.shape[0]):
        for j in range(c):
            if np.sum(U[i, :]) < 1 and j < c - 1:
                U[i, j] = random.uniform(0, 1 - np.sum(U[i, :]))
        # 最后一列单独赋值
        U[i, c - 1] = 1 - np.sum(U[i, :])
    iris = FCM(c, alpha, data, U)
    # print(sonar.U)
    # 保留更新前的中心点
    center = np.ones_like(iris.center)
    
    # 迭代部分
    i = 1
    while np.linalg.norm(iris.center-center)>0.0000001:
        center=iris.center.copy()
        # print(U)
        iris.update_center()
        # print(iris.center)
        iris.compute_J()
        # print(iris.center)
        iris.update_U()
        print("迭代第{}次".format(i))
        iris.compute_target()
        i+=1
        
    print("算法收敛,迭代终止")
    # print(iris.dis)
    # 计算轮廓系数来评价聚类效果
    print(metrics.silhouette_score\
          (data, iris.target.reshape(-1), metric='euclidean'))

sonar

from numpy.ma.core import argsort
from openpyxl import load_workbook
from sklearn import metrics
import numpy as np
import random


def get_data(path):
    """
    传入文件路径,返回样本数据以及标签
    """
    workbook = load_workbook(path)
    # print(workbook.sheetnames)
    sheet = workbook['Sheet1']
    # 存储样本特征集
    sonar = np.zeros([208, 60])
    # 存储标签
    target = np.zeros([208, 1])
    for i in range(208):
        for j in range(60):
            sonar[i, j] = sheet.cell(row=i + 1, column=j + 1).value
    # print(sonar.shape)
    # print(sonar[0,59])
    for i in range(208):
        if sheet.cell(row=i + 1, column=61).value == 'R':
            target[i] = 0
        else:
            target[i] = 1
    return sonar, target


# 编写FCM类
class FCM():
    def __init__(self, c, alpha, data, U):
        # 聚类数目
        self.c = c
        # 柔性参数(衡量模糊度)
        self.alpha = alpha
        # 样本数据
        self.data = data
        # 隶属度矩阵
        self.U = U
        # 聚类中心
        self.center = np.zeros([self.c, self.data.shape[1]])
        # 代价函数
        self.J = 0
        # 存储距离
        self.dis = np.zeros([self.data.shape[0], self.c])
        # 聚类结果
        self.target = np.zeros([self.data.shape[0], 1])


    # 更新聚类中心
    def update_center(self):
        for i in range(self.c):
            self.center[[i], :] = np.sum(np.power(\
                self.U[:,[i]],self.alpha)*self.data, \
                    axis=0) / np.sum(np.power(self.U\
                        [:, [i]], self.alpha))

    # 更新隶属度矩阵
    def update_U(self):
        # 计算距离,并存入self.dis中
        for i in range(self.dis.shape[0]):
            for j in range(self.dis.shape[1]):
                self.dis[i, j] = np.linalg.norm(\
                    self.data[[i], :] - self.center[[j], :])
        # print(self.dis)
        # 根据公式计算新的隶属度矩阵
        for i in range(self.U.shape[0]):
            for j in range(self.U.shape[1]):
                self.U[i, j] = 1 / np.\
                    sum(np.power(self.dis\
                        [i, j]/self.dis[[i], \
                            :], 2-(self.alpha-1)))
        # print(self.U)

    # 计算代价函数
    def compute_J(self):
        self.J = 0
        for j in range(self.c):
            J = 0
            for i in range(self.data.shape[0]):
                J += np.power(self.U[i, j], self.alpha)\
                    * np.power(self.dis[i, j], 2)
            # print(J)
            self.J += J

    # 根据当前隶属度矩阵判断类别
    def compute_target(self):
        for i in range(self.data.shape[0]):
            self.target[i, 0] = np.argmax(self.U[[i], :])


if __name__ == "__main__":
    # 读取数据
    path = r'sonar.xlsx'
    data, target = get_data(path)
    # print(data.shape)
    # 设置聚类数目
    c = 2
    # 设置柔性参数alpha
    alpha = 2
    # 随机生成隶属度矩阵U
    # 随机初始化隶属度矩阵
    U = np.zeros([data.shape[0], c])
    for i in range(data.shape[0]):
        for j in range(c):
            if np.sum(U[i, :]) < 1 and j < c - 1:
                U[i, j] = random.uniform(0, 1 - np.sum(U[i, :]))
        # 最后一列单独赋值
        U[i, c - 1] = 1 - np.sum(U[i, :])
    sonar = FCM(c, alpha, data, U)
    # print(sonar.U)
    # 保留更新前的中心点
    center = np.ones_like(sonar.center)

    i = 1
    while np.linalg.norm(sonar.center-center)>0.0000001:
    # while i<100:
        center=sonar.center.copy()
        # print(center)
        sonar.update_center()
        # print(sonar.center)
        sonar.compute_J()
        # print(sonar.center)
        sonar.update_U()
        print("迭代第{}次".format(i))
        sonar.compute_target()
        i+=1
    # print(sonar.target)
    print("迭代结束")
    # print(type(sonar.target))
    # print(center)
    print(metrics.silhouette_score(data, sonar.target.\
        reshape(-1), metric='euclidean'))

图像分割

import numpy as np
from PIL import Image
from sklearn import metrics
import random
from matplotlib import pyplot as plt


## 创建一个k-means类
class K_means():
    def __init__(self,data,center_points, k):
        self.center=center_points #聚类中心点的坐标
        self.k=k #聚类数
        self.data=data #样本
        # self.target=target #样本对应的标签
        # dis前k列表示每个样本到各聚类中心的距离,最后一列表示类别
        self.dis=np.zeros([self.data.shape[0],self.k+1])
    
    
    #distance方法:计算当前到k个聚类中心的距离,找出最小距离并将每个样本分类
    def distance(self): 
        for i in range(len(self.data)):
            for j in range(self.k):
                ## 距离取为欧氏距离,即为二范数
                self.dis[i,j]=np.linalg.norm\
                    (self.center[j,:]-self.data[i],ord=2)
            self.dis[i,self.k]=np.where(self.dis[i,:]==\
                np.min(self.dis[i,0:self.k]))[0][0]
   
        
    ## 根据self.dis更新聚类中心self.center    
    def update(self):
        for i in range(self.k):
            #设置统计数组每次循环将同类的样本存储起来
            count=np.zeros([1,self.data.shape[1]]) 
            # 逐类进行统计并计算新的聚类中心
            for j in range(len(self.data)):
                if self.dis[j,self.k]==i:
                    # 一行一行进行拼接
                    count=np.row_stack((count,self.data[j,:])) 
                else:
                    continue
            # 更新聚类中心
            # print(count.shape)
            self.center[[i],:]= np.mean(count,axis=0)
        # print(self.center)
            
## 输出分割后的图像
def plot_pic(*args):
    k=args[0]
    rows,columns=args[1],args[2]
    targets=args[3]
    img=np.zeros([rows,columns])
    # for i in range(rows):
    #     for j in range(columns):
    #         if 
    img=targets.reshape([columns,rows]).T
    # img=img*255/(k-1)
    im = Image.new("RGB",(rows,columns))#创建图片
    for i in range(rows):
        for j in range(columns):
            if img[i,j]==0:
                im.putpixel((i,j),(255,0,0))
            if img[i,j]==1:
                im.putpixel((i,j),(0,255,0))
            if img[i,j]==2:
                im.putpixel((i,j),(0,0,255))
    im.show()
    im.save('3-result.jpg')
    

if __name__ == "__main__":
    
    # 读取数据
    img=Image.open('3.bmp','r')
    lst=list(img.getdata())
    # print(type(img.size))
    data=np.array(lst)
    # 指定k-means中的k值
    k=3
    # 随机选择k个聚类中心
    flags= random.sample(range(0,data.shape[0]),k)
    # print(flags)
    center_points=np.zeros([k,data.shape[1]])
    j=0
    for i in range(k):
        center_points[[j],:]=data[flags[i],:]
        j+=1
    
    # 实例化一个K_means类
    pic= K_means(data,center_points, k)
    # print(pic.dis)
    
    # center用来和pic.center比较,相同则迭代终止
    center=np.zeros_like(pic.center)
    
    # 迭代部分
    i=0
    while np.linalg.norm(pic.center-center)>0.001:
        # 刚用class写,踩到的一个坑,这里要调用copy方法,
        # 直接等于表示只修改指针,不可取
        center=pic.center.copy()
        pic.distance()
        # print(center)
        pic.update()
        i+=1
        print("迭代第{}次".format(i))
        
    print("算法收敛,迭代终止")
    # print(pic.dis[:,-1])
    # 计算轮廓系数来评价聚类效果
    print(metrics.silhouette_score(data,\
        pic.dis[:,-1], metric='euclidean'))
    plot_pic(k,img.size[0],img.size[1],pic.dis[:,-1])

文章作者: Reset Ran
版权声明: 本博客所有文章除特別声明外,均采用 CC BY 4.0 许可协议。转载请注明来源 Reset Ran !
评论
  目录