以下面这道题为例：

把题目中的数据转换为文本文件Wm.txt（代码中会涉及到，也可以转化为其他格式，但我习惯用文本，比较简单方便）：

function [idxbest, Cbest,totsumDbest,Zbest] = myKmeans(k,start)
%idxbest:最优解的索引
%Cbest:最优解的聚类中心
%totsumDbest:最优解的所有类内距离和的总和，即准则函数的值
%X为样本参数矩阵
%k为聚类个数
%start人为选择的初始聚类中心的样本点的编号向量
if (k<1||k~=size(start,2)||k>30)
    errordlg('输入的参数不合法，请重新输入！','错误提示');
end
load Wm.txt;
X=Wm(:,2:3);
maxit=100;%所允许的最大迭代次数
[n, p] = size(X);%矩阵X的行列数
totsumD = 0;%最优解的所有类内距离和的总和，即准则函数的值
iter = 0;%初始迭代次数
C = X(start,:);%作为初始聚类中心的样本点的参数信息矩阵
ClusterBest= cell(3,1);%定义元胞数组，存放最后聚类的结果
ClusterBest=loopBody();%获取最后结果（调用了循环体函数）

%所输出的结果对应以上所定义的元胞数组的内容
idxbest = ClusterBest{3};
Cbest = ClusterBest{2};
totsumDbest = ClusterBest{1};
Zbest= Iterresult(idxbest,iter,Cbest,totsumDbest);


function cellout = loopBody()%循环体函数（类似于主函数）

cellout = cell(3,1);         
                              % cellout{1} = 总距离(准则函数的值）
                              % cellout{2} = 聚类中心
                              % cellout{3} = 每个样本点归属到某一个聚类中心的索引向量

                %计算从每个点到每个聚类中心的距离以及点到聚类的初始化
        D = distfun(X, C);%样本点距离每个聚类中心的距离矩阵
        [d, idx] = min(D, [], 2);%d记录D每行的最小值，idx记录最小值的列号
        m = accumarray(idx,1,[k,1]);%m为一个k行1列的矩阵，向量idx记录的位置的数值加1，即k个类其中每个所含样本点的个数
                totsumD = sum(D((idx-1)*n + (1:n)'));
            converged = batchUpdate();%调用更新类的函数（即主要迭代过程）


            if ~converged%类更新失败
                  errordlg('迭代过程出现异常！','错误提示');
            end

            %迭代结束后计算最后相关的参数
            nonempties = find(m>0);%找到不空的类的位置信息
            D(:,nonempties) = distfun(X, C(nonempties,:));%迭代完成后更新距离矩阵
            d = D((idx-1)*n + (1:n)');%每个样本点离其最近的聚类中心的距离向量
            sumD = accumarray(idx,d,[k,1]);%每个类的类内距离和
            totsumD = sum(sumD(nonempties));%准则函数的值                                                                                                                        

             cellout = {totsumD,C,idx}'; % 保存到目前为止的最佳解决方案
               
               
        function converged = batchUpdate()%更新类函数 完成迭代

            moved = 1:n;%先假设所有样本点都会发生移动
            changed = 1:k;%假设所有聚类中心都会发生改变
                        %idx为记录每个样本点离目前所有聚类中心距离最近的那个中心的编号（主函数里出现过，那时的聚类中心为人为输入的初始聚类中心，还没有开始迭代）
                        %idx会随着迭代发生变化，previdx则记录上一次的数据
                        previdx = zeros(n,1);%先用零矩阵对其进行初始化
            prevtotsumD = Inf;%Inf在matlab里表示无穷大
            iter = 0;%初始迭代次数为0

            Z=Iterresult(idx,iter,C,totsumD);

            converged = false;%初始返回值，迭代成功即变为true
            while true
                iter = iter + 1;

                %计算新的聚类中心以及其包含的样本点数量 并且更新每一个样本点到这些中心点的距离
                                %X为样本点的参数信息，idx为记录距离矩阵D每行最小值的列号。m表示每个类所含样本点的个数的变化情况，C为聚类中心的变化情况
                [C(changed,:), m(changed)] = gcentroids(X,idx,changed);
                D(:,changed) = distfun(X, C(changed,:));%重新计算每个样本点和新的聚类中心点之间的距离
                totsumD = sum(D((idx-1)*n + (1:n)'));%准则函数的值

                if prevtotsumD <= totsumD%判断准则函数是否收敛
                    idx = previdx;
                    [C(changed,:), m(changed)] = gcentroids(X, idx, changed);
                    iter = iter - 1;
                    break;
                end

                if iter >= maxit%迭代次数超过预定义次数，结束迭代
                    break;
                end

                % 确定离每个样本点最近的聚类中心，并将样本点重新分配给聚类中心
                previdx = idx;
                prevtotsumD = totsumD;
                [d, nidx] = min(D, [], 2);

                % 确定移动的样本点
                moved = find(nidx ~= previdx);
                if isempty(moved)
                    converged = true;%表示聚类不再发生变化，迭代成功，返回true，并结束迭代，回到循环体函数
                    break;
                end
                idx(moved) = nidx(moved);

                Z=Iterresult(idx,iter,C,totsumD);

                changed = unique([idx(moved); previdx(moved)])'; %查找获得或失去的样本点

            end

        end
end

function D = distfun(X, C)%计算距离矩阵（每个样本点到目前的每个聚类中心的距离）
[n,p] = size(X);
D = zeros(n,size(C,1));%D为n行k列的距离矩阵
nclusts = size(C,1);%值为k
%计算每个样本点和所传入的每个聚类中心的距离
        for i = 1:nclusts
            D(:,i) = (X(:,1) - C(i,1)).^2;
            for j = 2:p
                D(:,i) = D(:,i) + (X(:,j) - C(i,j)).^2;
            end
        end
end
function [centroids, counts] = gcentroids(X, index, clusts)%重新计算新的聚类的中心
p = size(X,2);%返回矩阵X的列数（本题应为2）
num = length(clusts);%返回聚类的个数（k）
centroids = NaN(num,p);%即k行2列的非数矩阵
counts = zeros(num,1);%k行1列的零矩阵
for i = 1:num
    members = (index == clusts(i));%回一个和index行列相同的矩阵，双等号成立则相应的位置为1，否则为零
    if any(members)%如果有非零的数存在于menmbers里
        counts(i) = sum(members);%计算到第i个聚类的距离最短的样本点的个数
        centroids(i,:) = sum(X(members,:),1) / counts(i);%重新计算新的聚类的中心
    end
end
end
function Z= Iterresult(idx,iter,C,totsumD)%输出每一次迭代结果，Z存放每次迭代的每个聚类中心所包含的样本点编号
   fprintf('第%d次迭代结果:\n',iter);
        Z=cell(k,1);%定义元胞数组，存放最后所得的每个聚类中心所包含的样本点的编号
for j=1:k
    q=1;
    for i=1:n
        if(idx(i)==j)
            Z{j}(q)=i;
            q=q+1;
        end
    end
end
for j=1:k
    fprintf('Z{%d}:',j);
    disp(Z{j});
end
   disp('聚类中心：');
   disp(C);
   disp('准则函数E的值：');
   disp(totsumD);
end
end