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

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

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

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

[cpp] view plaincopy

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <iostream>

using namespace std;

#define NA 60 /* 数据维数 */

#define K 6 /* 聚类数 */

#define Psize 36 /* 种群大小 */

#define T 30 /* 最大迭代数 */

#define ED 0.0000001 /* 结束条件 */

#define Min 1000000 /*最小值*/

#define MinCmp(a,b) (a<b?a:b)

#define INF 300000000

//记录每个点的坐标已经到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 为： *

* 1 2 3 4 *

* 1.0 1.2 1.3 1.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
=all_data[pop.cluster[i][index]].p
;

}

}

/********************************

* 输出信息函数 *

* 显示格式为： *

* 质心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
<<'/t';

cout<<endl;

}

else if(flag == 1){

cout<<"最终质心:"<<i+1<<endl;

for(n = 0; n < NA; n++)

cout<<pop.clu_cent[i].p
<<'/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
<<'/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;

}