模糊聚类算法的C语言实现六.docx

模糊聚类算法的C语言实现六.docx

《模糊聚类算法的C语言实现六.docx》由会员分享，可在线阅读，更多相关《模糊聚类算法的C语言实现六.docx（19页珍藏版）》请在冰豆网上搜索。

模糊聚类算法的C语言实现六

模糊聚类算法的C语言实现六

摘要

聚类分析在生产生活中有着非常广泛的应用比如学生成绩的划分粉煤灰分类等。

而模糊数学的理论也给我们提了更多的聚类方法如基于模糊等价关系的聚类分析方法和基于模糊划分的聚类方法等。

综合考虑二者的优劣势我们给出了一种较合理的算法即首先通过传递闭包法建立模糊等价矩阵针对每个阈值给出相应的分类；再应用F-统计量最大确定最佳阈值；将该阈值对应的分类作为初始值进行ISODATA修正。

根据此算法写出完整的C语言程序方便大量数据的处理。

并给出几个应用实例说明该算法在实际应用中有其合理性且利用计算机可以便捷的得到分类结果。

关键词

聚类分析模糊相似矩阵模糊等价矩阵F-统计量ISODATA方法阈值

模糊聚类的算法及程序实现

在模糊数学中聚类是一个很现实的应用。

常用的聚类方法有基于模糊等价关系的聚类分析方法和基于模糊划分的模糊聚类方法等两种算法各有优势和弱点。

传递闭包法可以根据阈值的取定给出不同的分类而阈值步长的任意性又对分类的合理性提出了挑战虽然我们可以用F-同计量方法确定最佳阈值但毕竟在决定阈值范围时就有很大的主观因素；ISODATA方法能够在退化情形得到最理想的c-模糊划分但方法的收敛快慢和最终分类的好坏过于依赖初始值的选取并且类数的选择对结果也很有影响。

结合两种算法我们给出了修正后较合理的方法即先由模糊等价关系聚类图中选择一最佳值以对应该阈值的分类作为初始值用ISODATA方法加以修正。

C语言程序：

#include<stdio.h>

#include<math.h>

#defineN100

#defineM100

#defineL110

tyedefstruct//定义结构体记录grou：

{

doubleline,F;//Line：

取定阈值F：

相应该种分类的F-统计量（吴震）

intnum,contain[N],mark[N],classify[N][N];//Num：

分类总数（组号为123……num）

}DATA;//Contain[i]：

第i类包含的素数（i=1,2,…,num）

//Mark[j]：

第j个样本相应的组号（j=1,2,…,n）

//Classify[i][j]：

第i类包含的第j素的下标

intn,m,ossible;

doublex[N][M],cArr[N][N],zero[M];//x:

每个对象由一组特征数据表征,cArr:

存储模糊相似（等价）矩阵

doublecoycArr[N][N];//zero:

辅助量每个指标均为零coycArr:

传递闭包法中的中转矩阵

intcBoolArr[N][N];

//取阈值limitLine的截关系存储至cBoolArr

doublecenter[N][M],overallcenter[M];//center:

聚类中心overallcenter:

总体样本的中心向量

doubleA0[N][M],A1[N][M],v[N][M];//A0:

c-模糊划分矩阵A1:

迭代所得c-模糊划分矩阵v:

聚类中心矩阵

DATAgrou[L],otGrou,reviseGrou;//grou[L]:

对于每一个阈值limitLine有一个分类otGrou:

最优分类

//reviseGrou:

ISODATA修正后得到的分类

voidinut（）//输入样本数据

{

inti,j;

scanf（"%d%d",&am;n,&am;m）;//n:

样本个数m:

每个样本的特征

for（i=1;i<=n;i++）

for（j=1;j<=m;j++）

scanf（"%lf",&am;x[i][j]）;//具体数据输入

return;

}

voidinitial（）//设置辅助量zero方便直接利用距离公式求得cArr初始值

{

intk;

for（k=1;k<=m;k++）

zero[k]=0;

return;

}

doubledistance2（doublea[],doubleb[]）//求向量a,b之间的距离

{

doubledis;

intk;

dis=0;

for（k=1;k<=m;k++）

dis+=（a[k]-b[k]）*（a[k]-b[k]）;

returndis;

}

voidcreat_conformArr（）//建立模糊相似矩阵

{

inti,j,k;

intnegative=0;//指示所得矩阵是否含有负值

doubledis_to_origin[N],tem;

for（i=1;i<=n;i++）

dis_to_origin[i]=sqrt（distance2（x[i],zero））;

for（i=1;i<=n;i++）//夹角余弦法确定x[i]与x[j]的相似程度cArr[i][j]

{

for（j=1;j<=n;j++）

{

tem=0;

for（k=1;k<=m;k++）

tem+=x[i][k]*x[j][k];

cArr[i][j]=tem/（dis_to_origin[i]*dis_to_origin[j]）;

}

}

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

if（cArr[i][j]<0）

{

negative=1;

break;

}

if（negative）//若cArr中出现负值则将矩阵全体重新调整为非负

{

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

{

cArr[i][j]++;

cArr[i][j]/=2;

}

}

return;

}

voidenclosure_algorithm（）//传递闭包法求模糊等价矩阵

{

inti,j,k;

intunchanged=0;//判断是否求得传递闭包

doublemergeVal;

while（unchanged==0）

{

unchanged=1;

for（i=1;i<=n;i++）//矩阵之间的自乘算

法

for（j=1;j<=n;j++）

{

mergeVal=0;

for（k=1;k<=n;k++）

if（cArr[i][k]<=cArr[k][j]&am;&am;cArr[i][k]>mergeVal）

mergeVal=cArr[i][k];

elseif（cArr[i][k]>cArr[k][j]&am;&am;cArr[k][j]>mergeVal）

mergeVal=cArr[k][j];

if（mergeVal>cArr[i][j]）

{

coycArr[i][j]=mergeVal;

unchanged=0;

}

else

coycArr[i][j]=cArr[i][j];

}

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

cArr[i][j]=coycArr[i][j];

}

return;

}

voidaroximate_closure（）//简化模糊等价矩阵

{

inti,j;

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

cArr[i][j]=floor（cArr[i][j]*10+0.5）/10;//每个四舍五入到小数点第一位

return;

}

inttwo_ower（int）//求2的整数次幂

{

inti,add;

add=1;

for（i=1;i<=;i++）

add*=2;

returnadd;

}

voidcluster（）//基于传递闭包的聚类过程

{

inti,j,k;

doublelimitLine,limitSte,minLimit;

intgrouNum;

intcmArr[N],grouIndentify[N],grouDelegate[N];

intflag,count;

inttem;

grouNum=1;

limitLine=1;//聚类阈值起始值

limitSte=0.1;//聚类阈值步长

minLimit=0.1;//聚类阈值下限

k=0;

while（limitLine>=minLimit）

{

k++;

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

if（cArr[i][j]>limitLine||fabs（cArr[i][j]-limitLine）<1e-5）

cBoolArr[i][j]=1;//凡是cArr[i][j]>=limitLine的位置在0-1矩阵中均置1

else

cBoolArr[i][j]=0;//而cArr[i][j]<limitLine的位置置为0

for（j=1;j<=n;j++）

{

tem=0;

for（i=1;i<=n;i++）

tem+=cBoolArr[i][j]*two_ower（i-1）;

cmArr[j]=tem;//cmArr:

将各列的值用唯一的二进制数表示方便比较

}

grouIndentify[1]=1;//当前该轮聚类分析中各素所在的组号

grouDelegate[1]=1;//当前该轮聚类分析中各组的代表素

grouNum=1;//当前该轮聚类分析中的分类数即最大组号

for（i=2;i<=n;i++）//算法的主要部分

{

flag=0;//当前素是否属于已分出的组

for（count=1;count<=grouNum;count++）//如果当前的素属于已分出的组则并出其所在的组

if（cmArr[grouDelegate[count]]==cmArr[i]）

{

grouIndentify[i]=count;

flag=1;

break;

}

if（flag==0）//如果当前的素不属于任何一份出的组则将其置入新组并选该作为新组的代表

{

grouNum++;

grouDelegate[grouNum]=i;

grouIndentify[i]=grouNum;

}

}

grou[k].line=limitLine;//对于每个聚类阈值确定其分组情况并存入相应grou[k]

grou[k].num=grouNum;

for（j=1;j<=n;j++）

grou[k].mark[j]=grouIndentify[j];

for（i=1;i<=grouNum;i++）

{

grou[k].contain[i]=0;

for（j=1;j<=n;j++）

if（grouIndentify[j]==i）

{

grou[k].classify[i][++grou[k].contain[i]]=j;

}

}

if（grouNum==1）

break;

limitLine-=limitSte;

}

ossible=k;

return;

}

voidredigest_grou（）//将相同的分组并入一个情况使grou[1]~grou[ossible]为互不相同的分组

{

inti,j,k;

intflag;

for（i=1;;i++）

{

flag=0;

j=i+1;

if（j>ossible）

break;

while（grou[i].num==grou[j].num）

{

j++;

flag=1;

}

if（flag）

{

for（k=j;k<=ossible;k++）

grou[k-j+i+1]=grou[k];

ossible=ossible-j+i+1;

}

if（i==ossible-1）

break;

}

}

doubleF_stat（DATAtemGrou）//计算聚类temGrou的F-统计量方法

{

inti,j,k;

doubleabove,below;

for（i=1;i<=temGrou.num;i++）

for（k=1;k<=m;k++）

{//第i类的聚类中心向量

center[i][k]=0;

for（j=1;j<=temGrou.contain[i];j++）

center[i][k]+=x[temGrou.classify[i][j]][k];

center[i][j]/=temGrou.contain[i];

}

for（k=1;k<=m;k++）//总体样本的中心向量

{

overallcenter[k]=0;

for（j=1;j<=n;j++）

overallcenter[k]+=x[j][k];

overallcenter[k]/=n;

}

if（temGrou.num==1||temGrou.num==n）

return0;

//求F值

above=0;//above:

分子表征类与类之间的距离

for（i=1;i<=temGrou.num;i++）

above+=temGrou.contain[i]*distance2（center[i],overallcenter）;

above/=temGrou.num-1;

below=0;//below:

分母表征类内样本间距离

for（i=1;i<=temGrou.num;i++）

for（j=1;j<=temGrou.contain[i];j++）

below+=distance2（x[temGrou.classify[i][j]],center[i]）;

below/=n-temGrou.num;

returnabove/below;

}

voidotlimitLine_select（）//确定最佳阈值

{

inti,ot;

doublemax;

for（i=1;i<=ossible;i++）

grou[i].F=F_stat（grou[i]）;

max=0;ot=1;

for（i=1;i<=ossible;i++）

if（grou[i].F>max）//F值越大说明分类越合理

{

max=grou[i].F;

ot=i;

}

otGrou=grou[ot];

return;

}

doublematrix_model（doublea[][N],doubleb[][N],ints,intt）//求矩阵范数|a-b|s行t列

{

inti,j;

doubledis;

dis=0;

for（i=1;i<=s;i++）

for（j=1;j<=t;j++）

dis+=（a[i][j]-b[i][j]）*（a[i][j]-b[i][j]）;

returnsqrt（dis）;

}

voidISOdata（DATAtemGrou）//将temGrou作为初始值进行ISODATA迭代

{

intc;

inti,j,k,rn=n;

intexcet[N],judge[N];

doubletem;

intis_nonemty[N];

intcount;

c=temGrou.num;

for（i=1;i<=n;i++）

excet[i]=0;

j=0;

for（i=1;i<=temGrou.num;i++）

{

if（temGrou.contain[i]==1）//对于仅有一个样本的类聚类前应该排除待聚类后再加上该类

{

excet[temGrou.classify[i][1]]=1;

c--;

rn--;

continue;

}

j++;

for（k=1;k<=n;k++）

if（temGrou.mark[k]==i）

A0[j][k]=1;

else

A0[j][k]=0;

}

for（i=1;i<=c;i++）//根据初始分类取初始c-模糊划分矩阵A0

{

k=0;

for（j=1;j<=n;j++）

if（excet[j]==0）

{

k++;

A0[i][k]=A0[i][j];

if（i==1）

judge[k]=j;

}

}

count=0;

while

（1）

{

count++;

for（i=1;i<=c;i++）

for（j=1;j<=rn;j++）

A1[i][j]=A0[i][j]*A0[i][j];

for（i=1;i<=c;i++）//计算当前迭代步的聚类中心矩阵

{

tem=0;

for（j=1;j<=rn;j++）

tem+=A1[i][j];

for（k=1;k<=m;k++）

{

v[i][k]=0;

for（j=1;j<=rn;j++）

v[i][k]+=A1[i][j]*x[j][k];

v[i][k]/=tem;

}

}

for（i=1;i<=c;i++）//修正A0得到新的A1

for（j=1;j<=rn;j++）

{

tem=0;

for（k=1;k<=c;k++）

tem+=1/distance2（x[j],v[k]）;

tem*=distance2（x[j],v[i]）;

A1[i][j]=1/tem;

}

if（matrix_model（A0,A1,c,rn）<=1e-5）//判断是否达到迭代出口

break;

for（i=1;i<=c;i++）

for（j=1;j<=rn;j++）

A0[i][j]=A1[i][j];

if（count>1000）

{

rintf（"OverFlow\n"）;

break;

}

}

rintf（"ste=%d\n",count）;

for（j=1;j<=rn;j++）//将模糊划分清晰化

{

tem=A0[1][j];k=1;

for（i=2;i<=c;i++）

if（A0[i][j]>tem）

{

tem=A0[i][j];

k=i;

}

reviseGrou.mark[judge[j]]=k;//按最大隶属度原则将x[judge[j]]归入组号为k的类

}

for（i=1;i<=c;i++）

{

is_nonemty[i]=0;//记录已分出类的样本数目

for（j=1;j<=rn;j++）

{

if（reviseGrou.mark[judge[j]]==i）

is_nonemty[i]++;//由于ISODATA给出的最终划分是退化的情形因此空集也可能成为子类

}

}

for（i=1;i<=c;）//重新调整组号使共有c组编号从1至c且每组非空

{

if（is_none

mty[i]==0）

{

for（k=1;k<=rn;k++）

if（reviseGrou.mark[judge[k]]>i）

reviseGrou.mark[judge[k]]--;

for（j=i;j<=c;j++）

is_nonemty[j]=is_nonemty[j+1];

c--;

}

else

i++;

}

k=c;//这时对应每一个组号的类都是非空的

for（i=1;i<=n;i++）

if（excet[i]==1）//将迭代前排除的类加上

{

reviseGrou.mark[i]=++k;

}

reviseGrou.num=k;//修正后得到的最终分类类数

//reviseGrou.mark标记在最终分类中每个样本所在的组

for（i=1;i<=reviseGrou.num;i++）

{

reviseGrou.contain[i]=0;

for（j=1;j<=n;j++）

if（reviseGrou.mark[j]==i）

{

reviseGrou.classify[i][++reviseGrou.contain[i]]=j;

}

}

return;

}

voidoutut（）

{

intj,k;

rintf（"%d\n",reviseGrou.num）;

for（j=1;j<=n;j++）

rintf（"x%d",j）;

rintf（"\n"）;

for（j=1;j<=n;j++）

rintf（"%d",reviseGrou.mark[j]）;

rintf（"\n"）;

for（j=1;j<=reviseGrou.num;j++）

{

rintf（"%d:

%d\n",j,reviseGrou.contain[j]）;

for（k=1;k<=reviseGrou.contain[j];k++）

rintf（"%d",reviseGrou.classify[j][k]）;

rintf（"\n"）;

}

return;

}

intmain（）

{

intk;

inti,j;

freoen（"in.txt","r",stdin）;

freoen（"out.txt","w",stdout）;

inut（）;//读入数据

initial（）;//初始化

creat_conformArr（）;//建立模糊相似矩阵

enclosure_algorithm（）;//求其传递闭包

aroximate_closure（）;//将得到的模糊等价矩阵的素四舍五入到一位小数

cluster（）;//对阈值为00.1~0.91的情形做聚类

redigest_grou（）;//简化分组使其不出现重复的分类方案