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

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

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

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

#include <stdio.h>#include <stdlib.h>#include <math.h>#include <iostream>using namespace std;#define NA60/* 数据维数 */#define K6/* 聚类数 */#define Psize36/* 种群大小 */#define T30/* 最大迭代数 */#define ED0.0000001/* 结束条件 */#define Min 1000000 /*最小值*/#define MinCmp(a,b) (a<b?a:b)#define INF300000000//记录每个点的坐标已经到K个中心点的距离typedef struct {double p[NA];double distance[K];}Point;//记录整个种群聚类的相关信息typedef struct {Point clu_cent[K];/* 即cluster_center 簇类中心 */int cluster[K][Psize];/* 簇类数组 */int cluster_num[K];/* 簇类中一组数据的编号 */double fitness;/* 样本适应度值,用于判断结束条件 */double old_fitness;/* 前一次迭代的适应度值 */double Je;/* 所有样本的平方误差和 */}Pop;/* 声明函数 */int Is_equal(int a[], int n, int b);double Euclid1(double x, double y);double dtw(int x, int y);void input_data();void Init_center();void calculate_distance();void Make_new_cluster();void Make_new_center();void output_info(int flag);Point all_data[Psize];/* 数据大小 */Pop pop;/************************************************ * 从外部文件导入数据，对于没有数据文件将报错， * * 数据文件的格式根据 NA 决定，例如NA = 4时，测 * * 试数据为四维，则test.data 为：* * 1234* * 1.01.21.31.4* *......* * ......* ***********************************************/double Dtwdistance(Point x, Point y){double distance[NA+1][NA+1];double output[NA+1][NA+1];int i,j;memset(distance,0,sizeof(distance));memset(output,0,sizeof(output));for (i=1;i<=NA;i++)for (j=1;j<=NA;j++)distance[i][j]=Euclid1(x.p[i-1],y.p[j-1]);output[1][1]=distance[1][1];for (i=0;i<=NA;i++){output[i][0]=INF;output[0][i]=INF;}for (i=2;i<=NA;i++)output[i][1]=output[i-1][1]+distance[i][1];for (i=2;i<=NA;i++)output[1][i]=output[1][i-1]+distance[1][i];for(i=2;i<=NA;i++)for(j=2;j<=NA;j++)output[i][j]=Mintt(Mintt(output[i-1][j-1],output[i][j-1]),output[i-1][j])+distance[i][j];//DP过程，计算DTW距离return output[NA][NA];}double dtw(int a, int b){if (a==b)return 0;return Dtwdistance(all_data[a],all_data[b]);}void input_data(){int i,j,temp; freopen("data.txt","r",stdin);for(i = 0; i < Psize; i++){cin>>temp;for (j=0;j<NA;j++){cin>>all_data[i].p[j];}}}/*************************************************** * 随机初始化聚类质心，当随机数有相同时跳过继续执行* **************************************************/void Init_center(){int i, j;int num = 0;int rand_num;int rand_num_tmp[K];/* 随机产生三个0～Psize的数 */while(num < K){rand_num = rand() % Psize;if(!Is_equal(rand_num_tmp, num, rand_num))rand_num_tmp[num++] = rand_num;}freopen("output.txt","w",stdout);cout<<"初始化重心为:"<<endl;for(i = 0; i < K; i++){cout<<rand_num_tmp[i]<<endl;for(j = 0; j < NA; j++){cout<<all_data[rand_num_tmp[i]].p[j]<<'/t';pop.clu_cent[i].p[j] = all_data[rand_num_tmp[i]].p[j];}cout<<endl;}}/********************************** * 检查数据是否有相等,相等返回1 * *********************************/int Is_equal(int a[], int n, int b){int i;for(i = 0; i < n; i++)if(a[i] == b) return 1;return 0;}/******************************************** * 计算Psize组数据到K个质心的欧几里德距离* *******************************************/void calculate_distance(){int i, j;for(i = 0; i < Psize; i++)for(j = 0; j < K; j++){all_data[i].distance[j] = dtw(i, j);}}/************************************************ * 此函数为欧几里德距离公式函数,此处用于计算* * 一组数据到对应簇中心的欧几里德距离。* ***********************************************/double Euclid1(double x, double y){return (y-x)*(y-x);}/************************ * 将数据进行簇集归类 * ***********************/void Make_new_cluster(){int i, j;double min;for(i = 0; i < K; i++)/* 初始化编号 */pop.cluster_num[i] = 0;for(i = 0; i < Psize; i++){int index = 0;min = all_data[i].distance[0];for(j = 1; j < K; j++){/* 筛选到簇心欧几里德最小的 */if(all_data[i].distance[j] < min){min = all_data[i].distance[j];index = j;}}/* 划分簇集 */pop.cluster[index][pop.cluster_num[index]++] = i;}/* 计算所有样本的平方误差和 */pop.Je = 0;for(i = 0; i < K; i++)for(j = 0; j < pop.cluster_num[i]; j++){/* 样本到簇心的值即为其欧几里德距离 */pop.Je +=pow(all_data[pop.cluster[i][j]].distance[i],2);}pop.old_fitness = pop.fitness;/* 前一次迭代适应度值 *///printf("old_fitness = %lf/n", pop.old_fitness);pop.fitness = pop.Je;/* 所有样本平方误差和即为适应度值 */}/************************************************* * 更新簇心，即求其群类的平均距离为新的簇心 * ************************************************/void Make_new_center(){int i, j, t,n;double tmp_sum;int index;double min=Min;for(i = 0; i < K; i++){tmp_sum = 0;for(j= 0;j< pop.cluster_num[i];j++)for (t=0;t<pop.cluster_num[i];t++){tmp_sum =dtw(pop.cluster[i][j],pop.cluster[i][t]); if(tmp_sum<min)min=tmp_sum;index=j;}for (n=0;n<NA;n++)pop.clu_cent[i].p[n]=all_data[pop.cluster[i][index]].p[n];}}/******************************** *输出信息函数* * 显示格式为：* * 质心K:* * NA维的质心数据* * 簇类K：* * NA维属于簇类K的数据* *......* *......* *******************************/void output_info(int flag){int i, j, n;for(i = 0; i < K; i++){if(flag == 0){cout<<"初始化质心:"<<i+1<<'/t'<<pop.cluster_num[i]<<endl;for(n = 0; n < NA; n++)cout<<pop.clu_cent[i].p[n]<<'/t';cout<<endl;} else if(flag == 1){ cout<<"最终质心:"<<i+1<<endl;for(n = 0; n < NA; n++)cout<<pop.clu_cent[i].p[n]<<'/t';cout<<endl;}cout<<"簇类:" <<i + 1<<endl;for(j = 0; j < pop.cluster_num[i]; j++){cout<<pop.cluster[i][j]<<'/t';for(n = 0; n < NA; n++)cout<<all_data[pop.cluster[i][j]].p[n]<<'/t';cout<<endl;}}}/******************************** *主函数* *******************************/int main(){int i,j;double differ = 1;/* 适应度差值 */int flag = 0;/* 用来标示是显示初始数据还是聚类后的数据 */freopen("data.txt","r",stdin);freopen("output.txt","w",stdout);input_data();/* 导入数据 */Init_center();/* 初始化质心 */for(i = 0; (i < T) && (differ > ED); i++){calculate_distance();/* 计算欧几里德距离 */Make_new_cluster();/* 产生新的聚类 */if(flag == 0){output_info(flag);/* 输出第一次聚类后的结果 */flag = 1;}Make_new_center();/* 对新的聚类产生新的质心 */differ = pop.old_fitness - pop.fitness;/* 判断条件 */differ = fabs(differ);//printf("differ = %lf/n", differ);//printf("old_fitness = %lf/n", pop.old_fitness);//printf("fitness = %lf/n", pop.fitness);}printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++/n");output_info(flag);/* 聚类后显示结果 */return 0;}