首页 > 代码库 > 通过K-MEDOIDS算法对时间序列进行聚类的实现

通过K-MEDOIDS算法对时间序列进行聚类的实现

2024-11-02 06:43:02   208人阅读

最近做数据挖掘相关的工作，题目是时间序列聚类研究，目前对于这方面的研究都还只是在起步阶段，被广泛使用的还是基于K-MEDOIDS的聚类，放弃K-MEANS的主要原因还是时间序列之间序列的计算难度，对于这方面我们也已经有了一定的进展，不过也还是有很多的问题。

把基于DTW与K-MEDOIDS的时间序列聚类的算法贴出来，希望对大家有些帮助吧。

这份代码是我在以前的代码的基础上直接改的，所以C和C++有些混用。

 

[cpp] view plaincopy
 
  1. #include <stdio.h>  
  2. #include <stdlib.h>  
  3. #include <math.h>  
  4. #include <iostream>  
  5. using namespace std;  
  6. #define NA  60      /* 数据维数 */  
  7. #define K   6       /* 聚类数 */  
  8. #define Psize   36      /* 种群大小 */  
  9. #define T   30      /* 最大迭代数 */  
  10. #define ED  0.0000001   /* 结束条件 */  
  11. #define Min 1000000   /*最小值*/  
  12. #define MinCmp(a,b) (a<b?a:b)  
  13. #define INF 300000000         
  14. //记录每个点的坐标已经到K个中心点的距离  
  15. typedef struct {  
  16.     double p[NA];  
  17.     double distance[K];  
  18. }Point;  
  19. //记录整个种群聚类的相关信息  
  20. typedef struct {  
  21.     Point clu_cent[K];  /* 即cluster_center 簇类中心 */  
  22.     int cluster[K][Psize];  /* 簇类数组 */  
  23.     int cluster_num[K]; /* 簇类中一组数据的编号 */  
  24.     double fitness;     /* 样本适应度值,用于判断结束条件 */  
  25.     double old_fitness; /* 前一次迭代的适应度值 */  
  26.     double Je;      /* 所有样本的平方误差和 */  
  27. }Pop;  
  28.   
  29. /* 声明函数 */  
  30. int Is_equal(int a[], int n, int b);  
  31. double Euclid1(double x, double y);  
  32. double dtw(int x, int y);  
  33. void input_data();  
  34. void Init_center();  
  35. void calculate_distance();  
  36. void Make_new_cluster();  
  37. void Make_new_center();  
  38. void output_info(int flag);  
  39. Point all_data[Psize];      /* 数据大小 */  
  40. Pop pop;  
  41. /************************************************ 
  42.  * 从外部文件导入数据，对于没有数据文件将报错， *  
  43.  * 数据文件的格式根据 NA 决定，例如NA = 4时，测 * 
  44.  * 试数据为四维，则test.data 为：     * 
  45.  *  1   2   3   4       * 
  46.  *  1.0 1.2 1.3 1.4     * 
  47.  *      ......              * 
  48.  *      ......              * 
  49.  ***********************************************/  
  50. double Dtwdistance(Point x, Point y)  
  51.     {  
  52.     double distance[NA+1][NA+1];  
  53.     double output[NA+1][NA+1];  
  54.     int i,j;  
  55.     memset(distance,0,sizeof(distance));  
  56.     memset(output,0,sizeof(output));  
  57.     for (i=1;i<=NA;i++)  
  58.         for (j=1;j<=NA;j++)  
  59.             distance[i][j]=Euclid1(x.p[i-1],y.p[j-1]);  
  60.     output[1][1]=distance[1][1];  
  61.     for (i=0;i<=NA;i++)  
  62.         {  
  63.         output[i][0]=INF;  
  64.         output[0][i]=INF;  
  65.         }  
  66.     for (i=2;i<=NA;i++)  
  67.         output[i][1]=output[i-1][1]+distance[i][1];  
  68.     for (i=2;i<=NA;i++)  
  69.         output[1][i]=output[1][i-1]+distance[1][i];  
  70.     for(i=2;i<=NA;i++)  
  71.         for(j=2;j<=NA;j++)  
  72.             output[i][j]=Mintt(Mintt(output[i-1][j-1],output[i][j-1]),output[i-1][j])+distance[i][j];  
  73.     //DP过程，计算DTW距离  
  74.     return output[NA][NA];  
  75.     }  
  76.   
  77. double dtw(int a, int b)  
  78.     {  
  79.     if (a==b)  
  80.     return 0;  
  81.     return Dtwdistance(all_data[a],all_data[b]);  
  82.     }  
  83.   
  84. void input_data()  
  85. {  
  86.         int i,j,temp;  
  87.         freopen("data.txt","r",stdin);  
  88.         for(i = 0; i < Psize; i++)  
  89.         {  
  90.             cin>>temp;  
  91.             for (j=0;j<NA;j++)  
  92.             {  
  93.             cin>>all_data[i].p[j];  
  94.             }             
  95.         }  
  96. }  
  97.   
  98. /*************************************************** 
  99.  * 随机初始化聚类质心，当随机数有相同时跳过继续执行* 
  100.  **************************************************/  
  101. void Init_center()  
  102. {  
  103.     int i, j;  
  104.     int num = 0;  
  105.     int rand_num;  
  106.     int rand_num_tmp[K];  
  107.     /* 随机产生三个0～Psize的数 */  
  108.     while(num < K){  
  109.         rand_num = rand() % Psize;  
  110.         if(!Is_equal(rand_num_tmp, num, rand_num))  
  111.             rand_num_tmp[num++] = rand_num;  
  112.     }  
  113.     freopen("output.txt","w",stdout);  
  114.     cout<<"初始化重心为:"<<endl;  
  115.     for(i = 0; i < K; i++){  
  116.     cout<<rand_num_tmp[i]<<endl;  
  117.         for(j = 0; j < NA; j++)  
  118.         {  
  119.                 cout<<all_data[rand_num_tmp[i]].p[j]<<‘/t‘;  
  120.                 pop.clu_cent[i].p[j] = all_data[rand_num_tmp[i]].p[j];  
  121.         }  
  122.         cout<<endl;  
  123.     }  
  124. }  
  125. /**********************************  
  126.  * 检查数据是否有相等,相等返回1 * 
  127.  *********************************/  
  128. int Is_equal(int a[], int n, int b)  
  129. {  
  130.     int i;  
  131.     for(i = 0; i < n; i++)  
  132.         if(a[i] == b) return 1;  
  133.     return 0;  
  134. }  
  135.   
  136. /******************************************** 
  137.  * 计算Psize组数据到K个质心的欧几里德距离* 
  138.  *******************************************/  
  139. void calculate_distance()  
  140. {  
  141.     int i, j;  
  142.     for(i = 0; i < Psize; i++)  
  143.         for(j = 0; j < K; j++){  
  144.             all_data[i].distance[j] = dtw(i, j);  
  145.         }  
  146. }  
  147.   
  148. /************************************************ 
  149.  * 此函数为欧几里德距离公式函数,此处用于计算* 
  150.  * 一组数据到对应簇中心的欧几里德距离。   * 
  151.  ***********************************************/  
  152. double Euclid1(double x, double y)  
  153.     {  
  154.     return (y-x)*(y-x);  
  155.     }  
  156. /************************ 
  157.  * 将数据进行簇集归类 * 
  158.  ***********************/  
  159. void Make_new_cluster()  
  160. {  
  161.     int i, j;  
  162.     double min;  
  163.       
  164.     for(i = 0; i < K; i++)       /* 初始化编号 */  
  165.         pop.cluster_num[i] = 0;  
  166.     for(i = 0; i < Psize; i++){    
  167.         int index = 0;  
  168.         min = all_data[i].distance[0];  
  169.         for(j = 1; j < K; j++){      /* 筛选到簇心欧几里德最小的 */  
  170.             if(all_data[i].distance[j] < min){  
  171.                 min = all_data[i].distance[j];  
  172.                 index = j;  
  173.             }  
  174.         }  
  175.         /* 划分簇集 */  
  176.         pop.cluster[index][pop.cluster_num[index]++] = i;     
  177.     }  
  178.     /* 计算所有样本的平方误差和 */  
  179.     pop.Je = 0;  
  180.     for(i = 0; i < K; i++)  
  181.         for(j = 0; j < pop.cluster_num[i]; j++){  
  182.         /* 样本到簇心的值即为其欧几里德距离 */  
  183.             pop.Je +=pow(all_data[pop.cluster[i][j]].distance[i],2);  
  184.         }  
  185.     pop.old_fitness = pop.fitness;  /* 前一次迭代适应度值 */  
  186. //  printf("old_fitness = %lf/n", pop.old_fitness);  
  187.     pop.fitness = pop.Je;   /* 所有样本平方误差和即为适应度值 */  
  188. }  
  189. /************************************************* 
  190.  * 更新簇心，即求其群类的平均距离为新的簇心  * 
  191.  ************************************************/  
  192. void Make_new_center()  
  193. {  
  194.     int i, j, t,n;  
  195.     double tmp_sum;  
  196.     int index;  
  197.     double min=Min;  
  198.     for(i = 0; i < K; i++)  
  199.     {  
  200.         tmp_sum = 0;  
  201.         for(j= 0;j< pop.cluster_num[i];j++)  
  202.         for (t=0;t<pop.cluster_num[i];t++)  
  203.         {  
  204.             tmp_sum =dtw(pop.cluster[i][j],pop.cluster[i][t]);   
  205.             if(tmp_sum<min)  
  206.                 min=tmp_sum;  
  207.             index=j;  
  208.         }  
  209.         for (n=0;n<NA;n++)  
  210.             pop.clu_cent[i].p[n]=all_data[pop.cluster[i][index]].p[n];  
  211.           
  212.     }  
  213. }  
  214.   
  215. /******************************** 
  216.  *  输出信息函数      *                
  217.  * 显示格式为：       * 
  218.  * 质心K:         * 
  219.  * NA维的质心数据     * 
  220.  * 簇类K：         * 
  221.  * NA维属于簇类K的数据      * 
  222.  *  ......          * 
  223.  *  ......          * 
  224.  *******************************/  
  225. void output_info(int flag)  
  226. {  
  227.     int i, j, n;  
  228.     for(i = 0; i < K; i++){    
  229.         if(flag == 0){  
  230.             cout<<"初始化质心:"<<i+1<<‘/t‘<<pop.cluster_num[i]<<endl;  
  231.             for(n = 0; n < NA; n++)  
  232.                 cout<<pop.clu_cent[i].p[n]<<‘/t‘;  
  233.             cout<<endl;  
  234.         }   
  235.         else if(flag == 1){   
  236.             cout<<"最终质心:"<<i+1<<endl;  
  237.             for(n = 0; n < NA; n++)  
  238.                 cout<<pop.clu_cent[i].p[n]<<‘/t‘;  
  239.             cout<<endl;  
  240.         }  
  241.         cout<<"簇类:" <<i + 1<<endl;  
  242.         for(j = 0; j < pop.cluster_num[i]; j++){       
  243.             cout<<pop.cluster[i][j]<<‘/t‘;  
  244.             for(n = 0; n < NA; n++)  
  245.                 cout<<all_data[pop.cluster[i][j]].p[n]<<‘/t‘;  
  246.                 cout<<endl;  
  247.             }  
  248.     }  
  249. }  
  250. /******************************** 
  251.  *      主函数     * 
  252.  *******************************/  
  253. int main()  
  254. {  
  255.     int i,j;  
  256.     double differ = 1;  /* 适应度差值 */  
  257.     int flag = 0;       /* 用来标示是显示初始数据还是聚类后的数据 */  
  258.     freopen("data.txt","r",stdin);  
  259.     freopen("output.txt","w",stdout);  
  260.     input_data();       /* 导入数据 */  
  261.     Init_center();      /* 初始化质心 */  
  262.     for(i = 0; (i < T) && (differ > ED); i++){  
  263.         calculate_distance();       /* 计算欧几里德距离 */  
  264.         Make_new_cluster();     /* 产生新的聚类 */  
  265.   
  266.         if(flag == 0){  
  267.             output_info(flag);  /* 输出第一次聚类后的结果 */  
  268.             flag = 1;     
  269.         }  
  270.         Make_new_center();      /* 对新的聚类产生新的质心 */  
  271.         differ = pop.old_fitness - pop.fitness; /* 判断条件 */  
  272.         differ = fabs(differ);  
  273.     //  printf("differ = %lf/n", differ);  
  274.     //  printf("old_fitness = %lf/n", pop.old_fitness);  
  275.     //  printf("fitness = %lf/n", pop.fitness);  
  276.     }  
  277.     printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++/n");  
  278.     output_info(flag);  /* 聚类后显示结果 */  
  279.     return 0;  
  280. }

 

通过K-MEDOIDS算法对时间序列进行聚类的实现


联系
我们