灰色关联度分析，灰色预测——基于matlab

1.引言

灰色系统理论建模的主要任务是根据具体灰色系统的行为特征数据，利用不多的数据中的显信息和隐信息，寻找因素间或因素本身的数学关系

2. 灰色关联度分析

背景：

通过对某健将级女子铅球运动员的跟踪调查，获得其 1982 年至 1986 年每年最好成绩及 16 项专项素质和身体素质的时间序列资料，见下表 ，试对此铅球运动员的专项成绩进行因素分析。

数据：

文本形式数据：

13.6 14.01 14.54 15.64 15.69
11.50 13.00 15.15 15.30 15.02
13.76 16.36 16.90 16.56 17.30
12.41 12.70 13.96 14.04 13.46
2.48 2.49 2.56 2.64 2.59
85 85 90 100 105
55 65 75 80 80
65 70 75 85 90
12.80 15.30 16.24 16.40 17.05
15.30 18.40 18.75 17.95 19.30
12.71 14.50 14.66 15.88 15.70
14.78 15.54 16.03 16.87 17.82
7.64 7.56 7.76 7.54 7.70
120 125 130 140 140
80 85 90 90 95
4.2 4.25 4.1 4.06 3.99
13.1 13.42 12.85 12.72 12.56

代码：

clc,clear 
load shuju.txt %把原始数据存放在纯文本文件 x.txt 中，其中把数据的"替换替换成. 
for i=1:15 
 shuju(i,:)=shuju(i,:)/shuju(i,1); %标准化数据
end 
for i=16:17 
 shuju(i,:)=shuju(i,1)./shuju(i,:); %标准化数据
end 
data=shuju; 
n=size(data,2); %求矩阵的列数，即观测时刻的个数
ck=data(1,:); %提出参考数列
bj=data(2:end,:); %提出比较数列
m2=size(bj,1); %求比较数列的个数
for j=1:m2 
 t(j,:)=bj(j,:)-ck; 
end 
mn=min(min(abs(t'))); %求最小差
mx=max(max(abs(t'))); %求最大差
rho=0.5; %分辨系数设置
ksi=(mn+rho*mx)./(abs(t)+rho*mx); %求关联系数
r=sum(ksi')/n %求关联度
[rs,rind]=sort(r,'descend') %对关联度进行排序

得到的关联度结果如下：

从上表可看出，全蹲（X13）与铅球专项（X0）关联度最大，为0.933。应该指出的是，公式中的| x0 (k) − xi(k) |不能区别因素关联是正关联还是负关联，这是灰色关联度分析的局限性所在。

3. 灰色预测

介绍

灰色预测是指利用 GM 模型对系统行为特征的发展变化规律进行估计预测，同时也可以对行为特征的异常情况发生的时刻进行估计计算，以及对在特定时区内发生事件的未来时间分布情况做出研究等等。
这些工作实质上是将“随机过程”当作“灰色过程”，“随机变量”当作“灰变量”，并主要以灰色系统理论中的 GM(1,1)模型来进行处理。
灰色预测在工业、农业、商业等经济领域，以及环境、社会和军事等领域中都有广泛的应用。特别是依据目前已有的数据对未来的发展趋势做出预测分析。

理论

代码

x0=[6093
6128
6144
6163
6228
6120
6110
6118
6135
6131
5957
5968
5988
6030
6083
6144
6196
6255
6324
6366
];
n=length(x0); 
lamda=x0(1:n-1)./x0(2:n); %计算级比
range=minmax(lamda');%计算级比的范围
%计算可容覆盖(a,b)
%c=exp(-2/(n+1));
%d=exp(2/(n+2));
%建立模型
x1=cumsum(x0) %累加运算
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)]; 
Y=x0(2:n); 
u=(inv(B'*B))*B'*Y
%n越大，外推越多
for t=0:n
    x(t+1)=(x0(1)-u(2)/u(1))*exp(-u(1)*t)+u(2)/u(1)%微分方程求解公式
    
end

for i=1:n
    z(1)=x(1)
    z(i+1)=x(i+1)-x(i)%微分方程还原，得到预测值
    
    %残差检验
    error(i)=abs(z(i)-x0(i))/x0(i)
    
end
plot([0:19],error(i),'*')%误差散点图
%y=vpa(z,6) %其中的 6 表示显示 6 位数字
%yuce=[x0(1),diff(yuce1)] %差分运算，还原数据
%误差