%% I. 清除环境变量
clear all
clc     
%% III. 导入数据  
%% 读入数据
file='表1.xlsx';          
%引号内为文件名，根据需要更改

%下列为读取文件中不同的数据，根据文件中的实际数据种类自行吗，命名。
%主要更改‘=’之前和‘()’内 B2:B3001 部分内容
T=xlsread(file,1,'C2:C1380');              %温度
V=xlsread(file,1,'P2:P1380');              %变形速率
R=xlsread(file,1,'Q2:Q1380');              %变形程度
P=xlsread(file,1,'O2:O1380');              
%将上述列矩阵合并为二维矩阵
input = [T,V,R];                %把以上输入参数放在同一个矩阵
output = [P];                %把以上输入参数放在同一个矩阵
[IP,TS]=mapminmax(input',0,1);
input=IP';   
x=input;                                    %读入数据
[n,p]=size(input);                          %确定输入参数的行数和列数，对应工作区的n值和p值
kmax=n;                                   %设定循环结束的次数，计算每行的伪F值和R2值，分成最大的类
format long                                 %设置显示字符串格式
pm=zeros(kmax,1);                           %给矩阵赋值为0
pm(1)=1;                                    %给定初始值，把第一行初始值赋值为1，否则分母为零，报错
for k=2:kmax                                %循环，分类最少是2类，
    d=pdist(x);                             %计算欧氏距离
    z1=linkage(d,'average');                %设置为类平均距离
    julei=cluster(z1,k);                    %将数据分割成簇进行聚类
    for t=1:k                               %计算k个类的类间偏差平方和的总合
        index_t=find(julei==t);
        size_t=length(index_t);
        a=x(index_t,:);
        pm(k)=sum((size_t-1)*var(a))+pm(k);
    end
end
Tm=sum(kmax*var(x));                       
bm=Tm-pm;
F=zeros(kmax,1);
for kk=2:kmax
    F(kk)=bm(kk)/pm(kk)*(n-kk)/(kk-1); %伪F值的计算，参考P209
end
R2=1-pm./Tm;%R2值的计算  参考P208
%%画图
figure()
plot((2:40),R2(2:40),'*');xlabel('分类数');ylabel('R2值');           %R2图

figure()
plot((2:40),F(2:40),'*-');xlabel('分类数');ylabel('F值');           %伪F值图
FIG=[(2:40)', R2(2:40), F(2:40)];