word linkage 选择合适的聚类个数matlab code

clear
load fisheriris
X = meas;
m  = size(X,2);
% load machine
%  load census
% % X = meas;
% X=X(1:2000,:);
d = pdist(X,'euclidean');

Z = linkage(d,'ward');%Create a hierarchical binary cluster tree using linkage
% [Hr,Tr]=dendrogram(Z, 0);  %generates a dendrogram plot of the hierarchical binary cluster tree
% h_gca = gca;
% h_gca.TickDir = 'out';
% h_gca.TickLength = [.002 0];
% h_gca.XTickLabel = [];
% c = cluster

val0=eval(vpa(Z(:,3),5));

t=1; n(1)=0; val=val0(1); f(1)=val;
for i=1:length(val0)
if val0(i)==val
n(t)=n(t)+1;
continue
else
t=t+1;
val=val0(i);
n(t)=1;
f(t)=val;
continue
end
end
figure(2)
fsum=cumsum(n);
x= length(val0)-fsum;
plot(x, f, 'o-','LineWidth', 2)
xlabel('Number of Clusters')
ylabel('Merge Distance')
title('A sample evalution graph')

% for i = 1: length(x)
%     y= cluster(Z, x(i));
%     B = Amean(X, y)
% end

x =x(end:-1:1);
f=f(end:-1:1);

i=0;
flag = 1;  b =  length(x); minerror_c_old = b-1;
while flag ~=0 & b>=20
RMSE = [];
for c = 2: (b-2)
x1 = x(1:c); y1 = f(1 : c);
x2 = x(c:b); y2 = f(c : b);
[~, s1] = polyfit(x1, y1,1);
[~, s2] = polyfit(x2, y2,1);
RMSE(c-1) = ((c-1)/(b-1))*s1.normr + ( (b-c)/(b-1))*s2.normr;
end
[~, minerror_c ]= min(RMSE);
minerror_c= minerror_c+1;
if  minerror_c >= minerror_c_old
break
else
minerror_c_old = minerror_c;
%           b = b -1;
b = 2*minerror_c ;
end
i=i+1
end

hold on
plot(x( minerror_c), f( minerror_c),'ro')
c =x( minerror_c)

labels = cluster(Z, c);

labs = unique(labels);
CORR = zeros();
for i = 1:length(labs)
corri = cov(X(labels == i, :));
CORR = CORR +  corri;
end