模式识别 之 模糊C均值聚类算法

完全照抄公式实现，没什么好讲的，有意的大家可以看： http://wenku.baidu.com/view/ee968c00eff9aef8941e06a2.html

下面是代码:

不知道为什么，在分类数，较多时候，容易出错！  

据论文上讲，算法的性能依赖于初始聚类中心。 所以下面这个算法，最好是，半交互，先指定聚类中心! 

void get_class_mohu_c_junzhi_julei(yangben_t*pYanBeng, int yangbengNum, yangben_t* pCentre,int classNum,double stopThreshold, double m)
{
assert(yangbengNum > classNum);
assert(classNum>=2);
assert(classNum<=20);
assert(stopThreshold>0.00000001);

//特征数量
int teZhengNum = pYanBeng[0].iTeZhengNum;

CString strLog;

//原始数据
// for (int n=0; n<yangbengNum; n++)
// {
// for (int i=0; i<teZhengNum; i++)
// {
// pYanBeng[n].pTeZheng[i];
// CString tmp;
// tmp.Format(" %7.3f ", pYanBeng[n].pTeZheng[i]);
// strLog += tmp;
// }
//
// strLog += "\n\r";
// }

////// 数据预处理 -- 归一化

yangben_t minYanBeng;
yangben_t maxYanBeng;

{
// 1. 取出最大最小值
for (int i=0; i<teZhengNum; i++)
{
double minVal = pYanBeng[0].pTeZheng[i];
double maxVal = pYanBeng[0].pTeZheng[i];

for (int n=1; n<yangbengNum; n++)
{
if (pYanBeng[n].pTeZheng[i] < minVal)
{
minVal = pYanBeng[n].pTeZheng[i];
}
if (pYanBeng[n].pTeZheng[i] > maxVal)
{
maxVal = pYanBeng[n].pTeZheng[i];
}
}
minYanBeng.pTeZheng[i] = minVal;
maxYanBeng.pTeZheng[i] = maxVal;
}

// 2. 归一化
for (int n=0; n<yangbengNum; n++)
{
for (int i=0; i<teZhengNum; i++)
{
pYanBeng[n].pTeZheng[i] = (pYanBeng[n].pTeZheng[i] - minYanBeng.pTeZheng[i])/(maxYanBeng.pTeZheng[i] - minYanBeng.pTeZheng[i]);
}
}

}

//归一化的逆运算
// for (int n=0; n<yangbengNum; n++)
// {
// for (int t=0; t<teZhengNum; t++)
// {
// pYanBeng[n].pTeZheng[t] = pYanBeng[n].pTeZheng[t] * (maxYanBeng.pTeZheng[t] - minYanBeng.pTeZheng[t]) + minYanBeng.pTeZheng[t];
// }
// }

// 归一化后的数据
strLog += "\n\r";
strLog += "归一化数据\n\r";
for (int n=0; n<yangbengNum; n++)
{
for (int i=0; i<teZhengNum; i++)
{
pYanBeng[n].pTeZheng[i];
CString tmp;
tmp.Format(" %7.5f ", pYanBeng[n].pTeZheng[i]);
strLog += tmp;
}

strLog += "\n\r";
}

//聚类中心
yangben_t*pJuLeiCent = pCentre;

//指定任意，聚类初始中心
time_t t;
time(&t);
srand(t);

// 随机聚类中心
for (int i = 0; i<classNum; i++)
{
int index = rand()%yangbengNum;

pJuLeiCent[i].iTeZhengNum = pYanBeng[index].iTeZhengNum;
for (int k=0; k<pYanBeng[index].iTeZhengNum; k++)
{
pJuLeiCent[i].pTeZheng[k] = pYanBeng[index].pTeZheng[k];
}

}

// 申请隶属矩阵
double *pLiShu = new double[yangbengNum * classNum];
if (pLiShu == NULL) return;

// 样本到各聚类中心的距离
double *pJuLi = new double[yangbengNum * classNum];
if (pJuLi == NULL) return;

double newSize = 10.0;
double oldSize = 0.0;

int loop = 0;
// 比较价值 -- 如果价值差异过大，则循环
while( abs(newSize - oldSize) > stopThreshold)
//while( loop < 50)
{
loop++;
//存上次计算价值
oldSize = newSize;

// 更新样本到聚类中心的距离
{
for (int i=0; i<yangbengNum; i++)
{
for (int c = 0; c<classNum; c++)
{
pJuLi[i*classNum + c] = get_distance_ou(&pYanBeng[i], &pJuLeiCent[c]);

if (pJuLi[i*classNum + c]<0.000000001)
{
pJuLi[i*classNum + c] = 0.000000001;
}
}
}

strLog += "\n\r";
strLog += "\n\r";
strLog += "距离列表\n\r";

for (int i=0; i<yangbengNum; i++)
{
for (int c = 0; c<classNum; c++)
{
CString tmp;
tmp.Format(" %7.5f ", pJuLi[i*classNum + c]);
strLog += tmp;
}
strLog += "\n\r";
}

}

{
// 更新隶属矩阵
for (int i=0; i<yangbengNum; i++)
{
for (int c = 0; c<classNum; c++)
{
double lishu = 0.0;

// 并不是完全隶属某一类
for (int k = 0; k<classNum; k++)
{
double tmp = 0.0;

tmp = ( pJuLi[i*classNum + c] / pJuLi[i*classNum + k] );

lishu += pow(tmp,2/(m - 1));
}
pLiShu[i*classNum + c] = 1.0 / lishu;


}
}
// 隶属矩阵不需要做归一化处理

}

// 更新聚类中心
{

for (int c=0; c<classNum; c++)
{
double tmp = 0.0;
double sum0 = 0.0;
yangben_t cenTmp;
for (int i=0; i<yangbengNum; i++)
{
tmp = pLiShu[i*classNum + c];
tmp = pow(tmp,m);
sum0 += tmp;
for (int j=0; j<teZhengNum; j++)
{
cenTmp.pTeZheng[j] += tmp * pYanBeng[i].pTeZheng[j];
}

}

/////////////=-========-------------------
for (int j=0; j<teZhengNum; j++)
{
pJuLeiCent[c].pTeZheng[j] = cenTmp.pTeZheng[j]/sum0;
}

}

strLog = "\n\r";
strLog += "\n\r";
strLog += "隶属矩阵\n\r";

for (int i=0; i<yangbengNum; i++)
{
CString tmp;
double sum = 0.0;
for (int c = 0; c<classNum; c++)
{
tmp.Format(" %14.8f ", pLiShu[i*classNum + c]);
strLog += tmp;

sum += pLiShu[i*classNum + c];
}
tmp.Format(" %14.8f \n\r", sum);
strLog += tmp;

}

strLog += "\n\r";
strLog += "\n\r";
strLog += "聚类中心\n\r";

for (int c=0; c<classNum; c++)
{
for (int i=0; i<teZhengNum; i++)
{
CString tmp;
tmp.Format(" %7.5f ", pJuLeiCent[c].pTeZheng[i]);
strLog += tmp;
}

strLog += "\n\r";
}
//AfxMessageBox(strLog);

}

// 更新价值
{
newSize = 0.0;
for (int c=0; c<classNum; c++)
{
for (int i=0; i<yangbengNum; i++)
{
double tmpLiShu = 0.0;
double tmpJuLi = 0.0;

tmpLiShu = pLiShu[i*classNum + c];
tmpJuLi = pJuLi[i*classNum + c];
newSize += pow(tmpLiShu,m) * tmpJuLi * tmpJuLi;
}
}
}

//break;
}

delete [] pJuLi;
delete [] pLiShu;

// 把聚类中心--做归一化的逆运算
for (int c=0; c<classNum; c++)
{
for (int t = 0; t<teZhengNum; t++)
{
pCentre[c].pTeZheng[t] = pCentre[c].pTeZheng[t] * (maxYanBeng.pTeZheng[t] - minYanBeng.pTeZheng[t]) + minYanBeng.pTeZheng[t];
}
}

strLog = "\n\r";
strLog += "\n\r";
strLog += "聚类中心\n\r";

for (int c=0; c<classNum; c++)
{
for (int i=0; i<teZhengNum; i++)
{
CString tmp;
tmp.Format(" %7.5f ", pJuLeiCent[c].pTeZheng[i]);
strLog += tmp;
}

strLog += "\n\r";
}

//AfxMessageBox(strLog);

}