经典算法之k-means聚类

原理

给定一系列样本点，想要分成k类

（1）首先从中随机选取k个点作为质心

（2）分别计算每个点到这k个质心的欧式距离（为例），到哪个质心距离最小，这个点就属于这个质心，这样把所有点分成k类

（3）再分别计算k类中每一类新的质心，根据公式2

（4）迭代（2）（3），直到达到最优条件，条件如最小距离平方和公式3；

公式1 

公式2

公式3

close all;clc;%第一类数据mu1=[0 0 0];  %均值S1=[0.3 0 0;0 0.35 0;0 0 0.3];  %协方差data1=mvnrnd(mu1,S1,100);   %产生高斯分布数据%%第二类数据mu2=[1.25 1.25 1.25];S2=[0.3 0 0;0 0.35 0;0 0 0.3];data2=mvnrnd(mu2,S2,100);%第三个类数据mu3=[-1.25 1.25 -1.25];S3=[0.3 0 0;0 0.35 0;0 0 0.3];data3=mvnrnd(mu3,S3,100);%显示数据plot3(data1(:,1),data1(:,2),data1(:,3),'+');hold on;plot3(data2(:,1),data2(:,2),data2(:,3),'r+');plot3(data3(:,1),data3(:,2),data3(:,3),'g+');grid on;%三类数据合成一个不带标号的数据类data=[data1;data2;data3];   %这里的data是不带标号的%k-means聚类[u re]=KMeans(data,3);  %最后产生带标号的数据，标号在所有数据的最后，意思就是数据再加一维度[m n]=size(re);%最后显示聚类后的数据figure;hold on;for i=1:m     if re(i,4)==1            plot3(re(i,1),re(i,2),re(i,3),'ro');     elseif re(i,4)==2         plot3(re(i,1),re(i,2),re(i,3),'go');     else          plot3(re(i,1),re(i,2),re(i,3),'bo');     endendgrid on;

KMeans.m

%N是数据一共分多少类%data是输入的不带分类标号的数据%u是每一类的中心%re是返回的带分类标号的数据function [u re]=KMeans(data,N)       [m n]=size(data);   %m是数据个数，n是数据维数    ma=zeros(n);        %每一维最大的数    mi=zeros(n);        %每一维最小的数    u=zeros(N,n);       %随机初始化，最终迭代到每一类的中心位置    for i=1:n       ma(i)=max(data(:,i));    %每一维最大的数       mi(i)=min(data(:,i));    %每一维最小的数       for j=1:N            u(j,i)=ma(i)+(mi(i)-ma(i))*rand();  %随机初始化，不过还是在每一维[min max]中初始化好些       end          end       while 1        pre_u=u;            %上一次求得的中心位置        for i=1:N            tmp{i}=[];      % 公式一中的x(i)-uj,为公式一实现做准备            for j=1:m                tmp{i}=[tmp{i};data(j,:)-u(i,:)];            end        end                quan=zeros(m,N);        for i=1:m        %公式一的实现            c=[];            for j=1:N                c=[c norm(tmp{j}(i,:))];            end            [junk index]=min(c);            quan(i,index)=norm(tmp{index}(i,:));                   end                for i=1:N            %公式二的实现           for j=1:n                u(i,j)=sum(quan(:,i).*data(:,j))/sum(quan(:,i));           end                   end                if norm(pre_u-u)<0.1  %不断迭代直到位置不再变化            break;        end    end        re=[];    for i=1:m        tmp=[];        for j=1:N            tmp=[tmp norm(data(i,:)-u(j,:))];        end        [junk index]=min(tmp);        re=[re;data(i,:) index];    end    end

C++实现

#include <iostream>#include <sstream>#include <fstream>#include <vector>#include <math.h>#include <stdlib.h>#define k 3//簇的数目using namespace std;//存放元组的属性信息typedef vector<double> Tuple;//存储每条数据记录int dataNum;//数据集中数据记录数目int dimNum;//每条记录的维数//计算两个元组间的欧几里距离double getDistXY(const Tuple& t1, const Tuple& t2) {double sum = 0;for(int i=1; i<=dimNum; ++i){sum += (t1[i]-t2[i]) * (t1[i]-t2[i]);}return sqrt(sum);}//根据质心，决定当前元组属于哪个簇int clusterOfTuple(Tuple means[],const Tuple& tuple){double dist=getDistXY(means[0],tuple);double tmp;int label=0;//标示属于哪一个簇for(int i=1;i<k;i++){tmp=getDistXY(means[i],tuple);if(tmp<dist) {dist=tmp;label=i;}}return label;}//获得给定簇集的平方误差double getVar(vector<Tuple> clusters[],Tuple means[]){double var = 0;for (int i = 0; i < k; i++){vector<Tuple> t = clusters[i];for (int j = 0; j< t.size(); j++){var += getDistXY(t[j],means[i]);}}//cout<<"sum:"<<sum<<endl;return var;}//获得当前簇的均值（质心）Tuple getMeans(const vector<Tuple>& cluster){int num = cluster.size();Tuple t(dimNum+1, 0);for (int i = 0; i < num; i++){for(int j=1; j<=dimNum; ++j){t[j] += cluster[i][j];}}for(int j=1; j<=dimNum; ++j)t[j] /= num;return t;//cout<<"sum:"<<sum<<endl;}void print(const vector<Tuple> clusters[]){for(int lable=0; lable<k; lable++){cout<<"第"<<lable+1<<"个簇："<<endl;vector<Tuple> t = clusters[lable];for(int i=0; i<t.size(); i++){cout<<i+1<<".(";for(int j=0; j<=dimNum; ++j){cout<<t[i][j]<<", ";}cout<<")\n";}}}void KMeans(vector<Tuple>& tuples){vector<Tuple> clusters[k];//k个簇Tuple means[k];//k个中心点int i=0;//一开始随机选取k条记录的值作为k个簇的质心（均值）srand((unsigned int)time(NULL));for(i=0;i<k;){int iToSelect = rand()%tuples.size();if(means[iToSelect].size() == 0){for(int j=0; j<=dimNum; ++j){means[i].push_back(tuples[iToSelect][j]);}++i;}}int lable=0;//根据默认的质心给簇赋值for(i=0;i!=tuples.size();++i){lable=clusterOfTuple(means,tuples[i]);clusters[lable].push_back(tuples[i]);}double oldVar=-1;double newVar=getVar(clusters,means);cout<<"初始的的整体误差平方和为："<<newVar<<endl; int t = 0;while(abs(newVar - oldVar) >= 1) //当新旧函数值相差不到1即准则函数值不发生明显变化时，算法终止{cout<<"第 "<<++t<<" 次迭代开始："<<endl;for (i = 0; i < k; i++) //更新每个簇的中心点{means[i] = getMeans(clusters[i]);}oldVar = newVar;newVar = getVar(clusters,means); //计算新的准则函数值for (i = 0; i < k; i++) //清空每个簇{clusters[i].clear();}//根据新的质心获得新的簇for(i=0; i!=tuples.size(); ++i){lable=clusterOfTuple(means,tuples[i]);clusters[lable].push_back(tuples[i]);}cout<<"此次迭代之后的整体误差平方和为："<<newVar<<endl; }cout<<"The result is:\n";print(clusters);}int main(){char fname[256];cout<<"请输入存放数据的文件名： ";cin>>fname;cout<<endl<<" 请依次输入: 维数 样本数目"<<endl;cout<<endl<<" 维数dimNum: ";cin>>dimNum;cout<<endl<<" 样本数目dataNum: ";cin>>dataNum;ifstream infile(fname);if(!infile){cout<<"不能打开输入的文件"<<fname<<endl;return 0;}vector<Tuple> tuples;//从文件流中读入数据for(int i=0; i<dataNum && !infile.eof(); ++i){string str;getline(infile, str);istringstream istr(str);Tuple tuple(dimNum+1, 0);//第一个位置存放记录编号，第2到dimNum+1个位置存放实际元素tuple[0] = i+1;for(int j=1; j<=dimNum; ++j){istr>>tuple[j];}tuples.push_back(tuple);}cout<<endl<<"开始聚类"<<endl;KMeans(tuples);return 0;}

引用

http://www.cnblogs.com/tiandsp/archive/2013/04/24/3040883.html

http://blog.csdn.net/qll125596718/article/details/8243404/