K-Means算法的C语言实现

K-means均值聚类算法

K-means 有一个著名的解释:牧师—村民模型:

有四个牧师去郊区布道,一开始牧师们随意选了几个布道点,并且把这几个布道点的情况公告给了郊区所有的村民,于是每个村民到离自己家最近的布道点去听课。 听课之后,大家觉得距离太远了,于是每个牧师统计了一下自己的课上所有的村民的地址,搬到了所有地址的中心地带,并且在海报上更新了自己的布道点的位置。 牧师每一次移动不可能离所有人都更近,有的人发现A牧师移动以后自己还不如去B牧师处听课更近,于是每个村民又去了离自己最近的布道点…… 就这样,牧师每个礼拜更新自己的位置,村民根据自己的情况选择布道点,最终稳定了下来。

我们可以看到该牧师的目的是为了让每个村民到其最近中心点的距离和最小。

截屏2022-02-23 下午3.52.39
  1. 选择初始化的 k 个样本作为初始聚类中心;(图1中红蓝绿三个点)
  2. 针对数据集中每个样本点计算它到 k 个聚类中心的距离并将其分到距离最小的聚类中心所对应的类中;(图2中红蓝绿三片区域中的点就是划分后的情况)
  3. 针对每个类别 ,重新计算它的聚类中心 (即属于该类的所有样本的质心);(图三中箭头移动的位置)
  4. 重复上面 2 3 两步操作,直到达到某个中止条件(迭代次数、最小误差变化等)。

在此,我们对1.2中四个步骤进行逐一分解进行说明:

1.选择初始化的 k 个样本作为初始聚类中心;(下图中红蓝绿三个点)

截屏2022-02-24 下午1.15.38

这个问题中设计了几个子问题:

要使用K-Means聚类算法,我们就应该明确样本点应该用什么数据类型来描述,通常对于聚类问题来说,我们使用到的样本点抽象为用笛卡尔坐标系来确定样本点之间的关系,然后用欧式距离(就是我们初中学过的坐标系上如何求两点之间的距离,如下图所示)来确定两点之间的距离。

截屏2022-02-23 下午4.19.01

因为样本之间的关系被抽象为坐标点之间的关系,所以我们使用二维数组来存储各个样本点。

当然,通常来说数据集都应该来自于某个具体的我们要求解的案例,但是对于我们对算法进行分析实现来说,为了快速获取原始数据,我们将直接使用随机函数(产生0到100之间的数)对这个二维数组进行填充。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
//随机生成volume个样本
void sample_generator(int data_set[volume][dimension])
{
    int i;
    int j;
    for (i = 0; i < volume; i++) {
        for (j = 0; j < dimension; j++) {
            data_set[i][j] = rand()%100;
        }
    }
}

1.设置几个中心点

为了方便展示,我们仅设置3个中心点,预计将50个点分为3类

2.中心点如何赋值

因为中心点要来回运动,每次聚类完中心点点位置都不一样,每次分类结果都有自己的中心点,所以我们干脆再数据集data_set[volume] []、data_set[volume+1] []、data_set[volume+2] []中设置旧的中心点(上一次的中心点)

在data_set[volume+3] []、data_set[volume+4] []、data_set[volume+5] []设置新的中新点(本次分类后的中心点)

3.中心点点初值如何设置

旧中心点的初值可以从数据集中随机抽取,新中心点的初值暂时设置为(0,0,0),(0,0,0),(0,0,0)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
/**
 * 初始化中心点位置
 * @param data_set 
 */
void centroid_Init(int data_set[volume+6][dimension])
{
    int i;
    int j;
    for (i = 0; i < 3; i++) {
        int random = rand()%volume;
        for (j = 0; j < 3; j++) {
            data_set[volume+i][j]=data_set[random][j];
            data_set[volume+i+3][j]=0;
        }
    }
}

2.针对数据集中每个样本xi计算它到 k 个聚类中心的距离并将其分到距离最小的聚类中心所对应的类中;(图2中红蓝绿三片区域中的点就是划分后的情况)

截屏2022-02-24 下午1.16.24

在每个样本xi计算它到 k 个聚类中心的距离中,我们要对每个样本点到中心点点距离进行欧式距离计算,所以要用到距离计算函数:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
/**
 * 计算数据集中任意两个点之间点距离
 * @param pointer1 想要计算的一个点 注意从第一个点从1开始
 * @param pointer2 想要计算点另一个点 注意第一个点从1开始
 * @param date_set 计算所需要点数据集
 * @return
 */
double Find_Distance(int pointer1 , int pointer2 ,int date_set[volume+6][dimension])
{
   double rij = sqrt((((date_set[pointer1-1][0])-(date_set[pointer2-1][0]))*((date_set[pointer1-1][0])-(date_set[pointer2-1][0])))+
           ((date_set[pointer1-1][1])-(date_set[pointer2-1][1]))*((date_set[pointer1-1][1])-(date_set[pointer2-1][1]))
           +((date_set[pointer1-1][2])-(date_set[pointer2-1][2]))*((date_set[pointer1-1][2])-(date_set[pointer2-1][2])));
    // 四舍五入,取整
    return rij;
}

归类问题我们可以转化为比较距离的问题,因为一个点到每个中心点的距离是不一样的,我们按照离哪个中心点距离最短的原则,把这个点化为离他最近的中心点那一个类。

我们定义A类——第51个点 B类——第52个点 C类——第53个点

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
/**
 * 归类算法
 * @param distance1 某样本点离A类的距离
 * @param distance2 某样本点离B类的距离
 * @param distance3 某样本点离C类的距离
 * @return
 */
int sort(double distance1,double distance2,double distance3)
{
    //如果离A类最近
    if (distance1<=distance2&&distance1<=distance3)
    {
        return 0;
    }
    //如果离B类最近
    else if (distance2<=distance1&&distance2<=distance3)
    {
        return 1;
    }
    //如果离C类最近
    else
    {
        return 2;
    }
}

我们使用三个数组来存储样本点的标号,每个数组大小设置为volume,将三个数组都初始化为-1。

在对样本点进行读取时,如果读取到的样本点序号为-1,说明读取完毕。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
/**
 * 给三个存储归类后样本点的数组并就行初始化
 * @param A
 * @param B
 * @param C
 */
void Centroid_Store_Init(int A[volume],int B[volume],int C[volume])
{
    int i;
    for (i = 0; i < volume; i++) {
        A[i]=-1;
        B[i]=-1;
        C[i]=-1;
    }
}

3.针对每个类别ai ,重新计算它的聚类中心 (即属于该类的所有样本的质心);(图三中箭头移动的位置)

截屏2022-02-24 下午1.16.45

每一轮所有点都完成分类后,三个类(A,B,C三个数组)中已经存放了分类好的点,那么如何计算新的中心点? 拿A类来说,最直接的办法就是把A类中所有点的x坐标相加求平均值,得到新的x坐标,所有点的y坐标相加求平均值,得到新的y坐标,所有点的z坐标相加求平均值,得到新的z坐标; 同样的方法,也可以求出B和C两个类的新中心点:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/**
 * 求新的中心
 * @param D 说明是哪个数据集
 * @param class A数据集用0 B数据集用1 C数据集用2表示
 * @param data_set 给定一个数据集
 */
void Centroid_recal(int D[volume],int class, double data_set[volume+6][dimension])
{
    int i = 0;
    double x=0;
    double y=0;
    double z=0;
    while (D[i]!=-1)
    {
        x+=data_set[D[i]][0];
        y+=data_set[D[i]][1];
        z+=data_set[D[i]][2];
        i++;
    }

    double average_x = x/i;
    double average_y = y/i;
    double average_z = z/i;

    //将点存储到新中心点中

        data_set[volume+3+class][0] = average_x;
        data_set[volume+3+class][1] = average_y;
        data_set[volume+3+class][2] = average_z;
}

4.重复上面 2 3 两步操作,直到达到某个中止条件(迭代次数、最小误差变化等)。

截屏2022-02-24 下午1.17.38

对于A类,计算A类新中心点与旧中心点的距离(即50与53求距离),若距离小于0.0005,则A类已经聚合; 对于B类,计算B类新中心点与旧中心点的距离(即51与54求距离),若距离小于0.0005,则B类已经聚合; 对于C类,计算C类新中心点与旧中心点的距离(即52与55求距离),若距离小于0.0005,则C类已经聚合; 当上边三个类都聚合后,数据集的分类已经完成,否则就开启下一轮迭代。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
/**
 * 如果返回0 说明聚合还未完成 需要二次聚类
 * @param date_set
 * @return
 */
int judge(double date_set[volume+6][dimension])
{
    double dis1 = Find_Distance(volume,volume+3,date_set);
    double dis2 = Find_Distance(volume+1,volume+1+3,date_set);
    double dis3 = Find_Distance(volume+2,volume+2+3,date_set);

    if (dis1<=0.0005&&dis2<=0.0005&&dis3<=0.0005)
    {
        return 1;
    }
    //聚合未完成 将新中心赋值给旧中心
    int i;
    int j;
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            date_set[volume+i][j] = date_set[volume+3+i][j];
        }
    }

    return 0;
}

计算出了新的中心点后,若还需进行下一轮迭代,则需要把这几个中心点赋值到旧中心点的位置上(一方面是计算需要,另一方面也是53,54,55位置要空出来留给下一次产生中心点的需要),通俗的说,就是50 = 53, 51 = 54, 52 = 55

1
2
3
4
5
6
7
8
    //聚合未完成 将新中心赋值给旧中心
    int i;
    int j;
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            date_set[volume+i][j] = date_set[volume+3+i][j];
        }
    }

到这里所有的重要函数就已经说明结束,下面给出整体的代码实现。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

//用来表述样本点的纬度 dimension==3表示坐标点为3维空间中点点
#define dimension 3

//用来定义样本点的数量 volume==50表示有50个样本点
#define volume 50

/**
 * 随机生成volume个样本
 * @param data_set
 */
void sample_generator(double data_set[volume][dimension])
{
    int i;
    int j;
    for (i = 0; i < volume; i++) {
        for (j = 0; j < dimension; j++) {
            data_set[i][j] = rand()%100;
        }
    }
}
/**
 * 初始化中心点位置
 * @param data_set
 */
void centroid_Init(double data_set[volume+6][dimension])
{
    int i;
    int j;
    int random;
    for (i = 0; i < 3; i++) {
        random = (rand())%volume;
        for (j = 0; j < 3; j++) {
            data_set[volume+i][j] = data_set[random][j];
            data_set[volume+i+3][j]=0;
        }
    }
}

/**
 * 计算数据集中任意两个点之间点距离
 * 但如果dimension变化,这个函数将失效,所以还需要进行修改
 * @param pointer1 想要计算的一个点 注意从第一个点从1开始
 * @param pointer2 想要计算点另一个点 注意第一个点从1开始
 * @param date_set 计算所需要点数据集
 * @return
 */
double Find_Distance(int pointer1 , int pointer2 ,double date_set[volume+6][dimension])
{
   double rij = sqrt((((date_set[pointer1][0])-(date_set[pointer2][0]))*((date_set[pointer1][0])-(date_set[pointer2][0])))+
           ((date_set[pointer1][1])-(date_set[pointer2][1]))*((date_set[pointer1][1])-(date_set[pointer2][1]))
           +((date_set[pointer1][2])-(date_set[pointer2][2]))*((date_set[pointer1][2])-(date_set[pointer2][2])));
    // 四舍五入,取整
    return rij;
}

/**
 * 归类算法
 * @param distance1 某样本点离A类的距离
 * @param distance2 某样本点离B类的距离
 * @param distance3 某样本点离C类的距离
 * @return
 */
int sort(double distance1,double distance2,double distance3)
{
    //如果离A类最近
    if (distance1<=distance2&&distance1<=distance3)
    {
        return 0;
    }
    //如果离B类最近
    else if (distance2<=distance1&&distance2<=distance3)
    {
        return 1;
    }
    //如果离C类最近
    else
    {
        return 2;
    }
}

/**
 * 给三个存储归类后样本点的数组并就行初始化
 * @param A
 * @param B
 * @param C
 */
void Centroid_Store_Init(int A[volume],int B[volume],int C[volume])
{
    int i;
    for (i = 0; i < volume; i++) {
        A[i]=-1;
        B[i]=-1;
        C[i]=-1;
    }
}

/**
 * 求新的中心
 * @param D 说明是哪个数据集
 * @param class A数据集用0 B数据集用1 C数据集用2表示
 * @param data_set 给定一个数据集
 */
void Centroid_recal(int D[volume],int class, double data_set[volume+6][dimension])
{
    int i = 0;
    double x=0;
    double y=0;
    double z=0;
    while (D[i]!=-1)
    {
        x+=data_set[D[i]][0];
        y+=data_set[D[i]][1];
        z+=data_set[D[i]][2];
        i++;
    }

    double average_x = x/i;
    double average_y = y/i;
    double average_z = z/i;

    //将点存储到新中心点中

        data_set[volume+3+class][0] = average_x;
        data_set[volume+3+class][1] = average_y;
        data_set[volume+3+class][2] = average_z;
}

/**
 * 打印数据集
 * @param data_set
 */
void DataPrint(double data_set[volume][dimension])
{
    int i;
    int j;
    printf("序 X     Y    Z\n");
    for (i = 0; i < volume+6; i++) {
            printf("%d ",i+1);
        for (j = 0; j < dimension; j++) {
            printf("%.1f ",data_set[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

void categorize(int A[volume],int B[volume],int C[volume],double data_set[volume+6][dimension]) {
    //初始化样本点存储数组
    Centroid_Store_Init(A, B, C);
//    int i;
//    for (i = 0; i < volume; i++) {
//        printf("%d",C[i]);
//    }
    //分别计算所有样本点到三个中心点的距离
    int i;
    int j = 0;
    int k = 0;
    int l =0;
    for (i = 0; i < volume; i++) {
        //计算所有样本点i到A类中心点的距离
        double a = Find_Distance(i, volume, data_set);
        //计算样本点i到B类中心点的距离
        double b = Find_Distance(i, volume + 1, data_set);
        //计算样本点i到C类中心点的距离
        double c = Find_Distance(i, volume + 2, data_set);

        //判断该样本点应该属于哪一类
        int class = sort(a, b, c);
        //该样本点属于A类
        if (class == 0) {
            A[j] = i;
            j++;
        }
            //样本点属于B类
        else if (class == 1) {
            B[k] = i;
            k++;
        }
            //样本点属于C类
        else {
            C[l] = i;
            l++;
        }
    }
}


/**
 * 打印每一类归类后的数据集元素
 * @param Class
 */
void Class_Print(int Class[volume],int class,double data_set[volume+6][dimension])
{

    //该样本点属于A类
    if (class == 0) {
        printf("A类样本点为\n");
    }
        //样本点属于B类
    else if (class == 1) {
        printf("B类样本点为\n");
    }
        //样本点属于C类
    else {
        printf("C类样本点为\n");
    }
    int i = 0;
    while (Class[i]!=-1)
    {
        printf("数据集中的第%d个样本点:",Class[i]+1);
        int j;
        for (j = 0; j < 3; j++) {
            printf("%.2f ",data_set[Class[i]][j]);
        }
        printf("\n");
        i++;
    }
    printf("共计%d个样本点\n",i);


}


/**
 * 如果返回0 说明聚合还未完成 需要二次聚类
 * @param date_set
 * @return
 */
int judge(double date_set[volume+6][dimension])
{
    double dis1 = Find_Distance(volume,volume+3,date_set);
    double dis2 = Find_Distance(volume+1,volume+1+3,date_set);
    double dis3 = Find_Distance(volume+2,volume+2+3,date_set);

    if (dis1<=0.0005&&dis2<=0.0005&&dis3<=0.0005)
    {
        return 1;
    }
    //聚合未完成 将新中心赋值给旧中心
    int i;
    int j;
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            date_set[volume+i][j] = date_set[volume+3+i][j];
        }
    }
    return 0;
}

/**
 * 用于测试聚类效果
 */
void test()
{

    //播种
    srand(time( NULL));
    //用来存放3个分类结果
    int A[volume];
    int B[volume];
    int C[volume];
    //声明数据集
    double data_set[volume+6][dimension];
    //生成一个数据集
    sample_generator(data_set);
    //初始化数据集的中心位置
    centroid_Init(data_set);

    //对结果进行分类
    categorize(A,B,C,data_set);
    //计算A类新的中心点 保存到数据集的 53 54 55 三个位置
    Centroid_recal(A,0,data_set);

    //计算B类新点中心点
    Centroid_recal(B,1,data_set);
    //DataPrint(data_set);
    //计算C类新的中心点
    Centroid_recal(C,2,data_set);
    DataPrint(data_set);

    printf("第1次分类后的结果为\n");
    //打印分类后各类的样本点
    Class_Print(A,0,data_set);
    Class_Print(B,1,data_set);
    Class_Print(C,2,data_set);
    printf("\n");

    int k = 2;
    //如果聚合不满足要求 进行再次聚合
    while (!judge(data_set))
    {
        //对结果进行分类
        categorize(A,B,C,data_set);
        // DataPrint(data_set);
        //计算A类新的中心点
        Centroid_recal(A,0,data_set);
        //DataPrint(data_set);
        //计算B类新点中心点
        Centroid_recal(B,1,data_set);
        //DataPrint(data_set);
        //计算C类新的中心点
        Centroid_recal(C,2,data_set);
        //DataPrint(data_set);

        printf("第%d次分类后的结果为\n",k);
        //打印分类后各类的样本点
        Class_Print(A,0,data_set);
        Class_Print(B,1,data_set);
        Class_Print(C,2,data_set);
        printf("\n");
        k++;
    }
}

int main() {
    test();
    return 0;
}

https://en.wikipedia.org/wiki/K-means_clustering

https://blog.csdn.net/Lanyan9/article/details/121993064