一、灰度关联分析

• 在数学建模中，灰色关联分析模型（Grey Relational Analysis, GRA） 是一种用于分析 “部分信息已知、部分信息未知” 的 “灰色系统” 中各因素间关联程度的量化方法。它通过比较参考序列与比较序列之间的 “几何相似性”，来判断因素对系统的影响强度，适用于数据量少、信息不完整或非线性关系的场景。对样本量的多少和样本有无规律都同样适用。
• 若两个序列在变化趋势上越接近（即曲线形状越相似），则它们的关联度越高，说明对应的因素对系统的影响越显著。

二、关联分析步骤

1. 母序列 (又称参考序列、母指标)
   能反映系统行为特征的数据序列，类似于因变量(Y)，记为：
   $$Y=[y_1,y_2,\cdots ,y_n{]}^T$$

2. 子序列（又称比较序列、子指标）
   影响系统行为的因素组成的数据序列，类似于自变量(X)，记为：
   $$X_{nm}=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1m}\\ x_{21}& x_{22}& \cdots & x_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nm}\end{array}\right]$$

3. 数据预处理
   由于不同要素具有不同量纲和数据范围，因此需要对其进行预处理以去除量纲，统一到近似范围内
   处理方法为：先求出每个指标的均值，再用指标中的元素除以其均值
   $$\widetilde{y_k}=\frac{y_k}{\overline{y_i}},\overline{y_i}=\frac{1}{n}\sum _{k=1}^n{y}_k$$
   $$\widetilde{x_{ki}}=\frac{x_{ki}}{\overline{x_i}},\overline{x_i}=\frac{1}{n}\sum _{k=1}^n{x}_{ki}(i=1,2,\cdots ,m)$$

4. 计算灰色关联系数
   计算子序列中各个指标与母序列的关联系数
   记为： $a={\min }_i{\min }_k\left|x_o(k)-x_i(k)\right|,b={\max }_i{\max }_k\left|x_o(k)-x_i(k)\right|$
   为两极最小差和最大差
   构造： ${\xi }_i(k)=y(x_o(k),x_i(k))=\frac{a+\rho b}{|x_o(k)-x_i(k)|+\rho b}$
   其中$\rho$为分辨系数，一般取(0.5)

5. 计算关联度
   $$r_i=\frac{1}{n}\sum _{k=1}^n{\xi }_i(k)=\frac{1}{n}\sum _{k=1}^{n}y(x_0(k),x_i(k))$$

6. 比较
   根据关联度 $r_i$ 的大小排序，关联度越高的因素，对母序列的影响越显著。

三、MATLAB实现

% 灰色关联分析用于系统分析例题
clear;clc
% A=[55 24 10;65 38 22;75 40 18;100 50 20]	%计算示例
A=input('请输入初始矩阵=') % 输入初始矩阵
Mean = mean(A); % 求出每一列的均值以提供后续的数据预处理
A = A ./ repmat(Mean,size(A,1),1); % size(A,1)=4, repmat(Mean,4,1)可以将矩阵按维度复制
disp('预处理后的矩阵为：'); disp(A)
Y = A(:,1); % 母序列 ，取第一行
X = A(:,2:end); % 子序列 ，取第二行至最后一行
absX0_Xi = abs(X - repmat(Y,1,size(X,2))) % 计算|X0 - Xi|矩阵(在这里我们把X0定义为了Y)
a = min(min(absX0_Xi)) % 计算两级最小差a 
b = max(max(absX0_Xi)) % 计算两级最大差b
rho = 0.5; % 分辨系数取0.5
gamma = (a+rho*b) ./ (absX0_Xi + rho*b) % 计算子序列中各个指标与母序列的关联系数
disp('子序列中各个指标的灰色关联度分别为：')
disp(mean(gamma))

四、Python实现

import numpy as np
# 请输入初始矩阵，示例为：[[55, 24, 10], [65, 38, 22], [75, 40, 18], [100, 50, 20]]
A = np.array(eval(input("请输入初始矩阵=")))  # 假设输入的形式是合法的Python列表形式
# 求出每一列的均值以供后续的数据预处理
Mean = np.mean(A, axis=0)
# 预处理后的矩阵
A_norm = A / Mean
print('预处理后的矩阵为: ')
print(A_norm)
# 母序列
Y = A_norm[:, 0]
# 子序列
X = A_norm[:, 1:]
# 计算|X0 - Xi|矩阵(在这里我们把X0定义为了Y)
absX0_Xi = np.abs(X - np.tile(Y.reshape(-1, 1), reps=(1, X.shape[1])))
# 计算两级最小差a
a = np.min(absX0_Xi)
# 计算两级最大差b
b = np.max(absX0_Xi)
# 分辨系数取0.5
rho = 0.5
# 计算子序列中各个指标与母序列的关联系数
gamma=(a+rho*b)/(absX0_Xi+rho*b)
print('子序列中各个指标的灰色关联度分别为：')
print(np.mean(gamma,axis=0))

五、练习

clear;clc

% 初始化三维数组
data = zeros(12, 3, 1);

% 读取数据并检查
try
    data(:,:,1) = readmatrix('灰色关联分析.xlsx','Sheet','sheet1','Range','B72:D83');
    disp('数据读取成功！');
    disp(['数据维度: ', num2str(size(data(:,:,1)))]);
catch err
    disp('数据读取失败:');
    disp(err.message);
    return;  % 读取失败则退出
end

% 灰色关联分析
i = 1;  % 简化循环，直接使用i=1

% 计算均值并检查
Mean = mean(data(:,:,i));
disp('每列均值:');
disp(Mean);

% 检查是否有零均值
if any(Mean == 0)
    disp('警告：存在零均值，可能导致归一化错误！');
end

% 数据归一化
data(:,:,i) = data(:,:,i) ./ repmat(Mean, size(data(:,:,i),1), 1);

% 定义母序列和子序列
Y = data(:,1,i);  % 母序列（第一列）

% 正确提取子序列：从第2列到最后一列
num_cols = size(data(:,:,i), 2);
if num_cols >= 2
    X = data(:, 2:num_cols, i);  % 子序列（第2列及以后）
    disp(['子序列维度: ', num2str(size(X))]);
else
    disp('错误：数据列数不足，无法提取子序列！');
    disp(['数据总列数: ', num2str(num_cols)]);
    return;
end

% 计算绝对差值矩阵
absX0_Xi = abs(X - repmat(Y,1,size(X,2)));

% 计算两级最小差和最大差
a = min(min(absX0_Xi));
b = max(max(absX0_Xi));

disp(['两级最小差a: ', num2str(a)]);
disp(['两级最大差b: ', num2str(b)]);

% 计算关联系数
rho = 0.5;  % 分辨系数
gamma = (a + rho*b) ./ (absX0_Xi + rho*b);

% 修复逻辑运算符错误：使用元素级 | 代替 ||
if any(isnan(gamma(:)) | isinf(gamma(:)))
    disp('警告：关联系数中存在NaN或Inf值！');
end

% 输出关联度
disp('子序列中各个指标的灰色关联度分别为：')
disp(mean(gamma))

% 绘图
x = 1:12;
figure(i)
hold on
plot(x, data(:,1,i), 'k', 'LineWidth', 1, 'DisplayName', '水生根茎类季度平均销量')
plot(x, data(:,2,i), 'r', 'LineWidth', 1, 'DisplayName', '季度各蔬菜品类的单品蔬菜销量种数')
plot(x, data(:,3,i), 'b', 'LineWidth', 1, 'DisplayName', '特殊节气期间特征单品蔬菜销量')
set(gca, 'linewidth', 1)
set(gca, 'Box', 'on')
legend()  % 添加图例
title('母序列与子序列对比图')
xlabel('样本点')
ylabel('归一化值')
hold off

运行结果如下图：