引言

本篇文章将具体讲解笔者是如何复现2023年数学建模国赛c题第一问的，这里以论文C50为例。

C050.pdf

正文

首先我们打开数学建模官方网站，下载题目明确需求，下面是第一问的问题：

论文中首先进行了描述性统计分析，这要求我们首先需要对数据进行预处理，首先按照论文中所说，我先将每日销售量小于0的数据全部剔除，然后我使用excel对附件1和附件2做了合并处理，主要使用的方法是“筛选”“数据透视表”“排序”，最后经过处理之后，我得到了六大类蔬菜的每日销售量和每个单品的每日销售量，但是这时候我发现表中缺少了一些时间，发现的步骤也很简单，我首先计算了2020年7月1日到2022年6月30日的相隔天数，对比发现表中行数较少，得出结论缺少了10天，于是我使用spss对数据进行了插值填充，（这里只对六大类进行了填充，单品并没有进行填充，原因相信大家也能理解。）最后得出的数据格式如下所示：

同样的，我把相关数据放在下面：

单品每日销售量.xlsx六大类蔬菜每日销售量.xlsx

然后我通过matlab代码使用这些数据做了描述性统计，代码见下：

Test=xlsread('单品蔬菜日销售量.xlsx','sheet1','B2:EM1096')
%% 统计描述
MIN = min(Test); % 每一列最小值
MAX = max(Test); % 每一列最大值
MEAN = mean(Test); % 每一列均值
MEDIAN = median(Test); %每一列中位数
SKEWNESS = skewness(Test); %每一列偏度
KURTOSIS = kurtosis(Test); %每一列峰度
STD = std(Test); % 每一列标准差
RESULT = [MIN;MAX;MEAN;MEDIAN;SKEWNESS;KURTOSIS;STD] %矩阵表示统计量
sum=sum(Test);
Test=xlsread('各类每日销售量.xlsx','sheet2','B2:G1096')
%% 统计描述
MIN = min(Test); % 每一列最小值
MAX = max(Test); % 每一列最大值
MEAN = mean(Test); % 每一列均值
MEDIAN = median(Test); %每一列中位数
SKEWNESS = skewness(Test); %每一列偏度
KURTOSIS = kurtosis(Test); %每一列峰度
STD = std(Test); % 每一列标准差
RESULT = [MIN;MAX;MEAN;MEDIAN;SKEWNESS;KURTOSIS;STD] %矩阵表示统计量

运行之后就可以得出论文中的表格的数据：

然后根据处理之后的excel可以画出折线图，人为的说一下各大品类和单品的分布规律即可。

接下来我们开始计算各大品类的相关性。

这里的思路是首先看一下数据是不是符合Pearson的要求：即数据是不是线性的以及数据是不是符合正态分布的，首先我使用spss画出了数据的矩阵散点图,如下所示：

可以看到数据勉强还是有线性关系的，接下来我使用JB检验验证数据是否符合正态分布，这里使用matlab实现，代码运行之后的结果如下：

可以看到数据不符合正态分布，因此不能选用Pearson，故使用Spearman,依旧通过matlab运行代码实现，代码如下所示：

RX = xlsread('各类每日销售量.xlsx', 'sheet2', 'B2:G1096');
[R, P] = fun_spearman(RX, 1);
function [p] = calculate_p(r, m, kind)
    z = abs(r) * sqrt(m - 1); % 计算检验值
    p = (1 - normcdf(z)) * kind; % 计算 p 值
end
function [r] = calculate_r(X, Y)
    RX = rank_data(X); % 计算 X 的等级
    RY = rank_data(Y); % 计算 Y 的等级
    d = RX - RY; % 计算 X 和 Y 等级差
    n = size(X, 1); % 计算样本个数 n
    r = 1 - (6 * sum(d .* d)) / (n * (n^2 - 1)); % 计算斯皮尔曼相关系数
end
function [R, P] = fun_spearman(X, kind)
    if nargin == 1 % 判断用户输入的参数
        kind = 2;
    end
    [m, n] = size(X); % 计算样本个数和指标个数
    if m < 30 % 判断是否样本数量
        disp('样本个数少于 30，请直接查临界值表进行假设检验');
    elseif n < 2 % 判断是否指标数太少
        disp('指标个数太少，无法计算');
    elseif kind ~= 1 && kind ~= 2 % 判断 kind 是否为 1 或者 2
        disp('kind 只能取 1 或者 2');
    else
        R = ones(n); % 初始化 R 矩阵
        P = ones(n); % 初始化 P 矩阵
        for i = 1:n
            for j = (i + 1):n
                r = calculate_r(X(:, i), X(:, j)); % 计算 i 和 j 两列的相关系数 r
                p = calculate_p(r, m, kind); % 计算 p 值
                R(i, j) = r; R(j, i) = r; % 求得相关系数
                P(i, j) = p; P(j, i) = p; % 求得检验 p 值
            end
        end
    end
end

代码运行之后我们可以得到两个表格，一个是p表格，一个是r（相关系数），为了美化表格，我在SPSSPRO上面花了相关系数的热力图，如下：

根据这两个表格我们就分析完成了各大种类的相互关系，接下来我们继续分析单品之间的相互关系。

这里首先需要使用excel处理一下得出单品对应的总销售量、每日最大销售量、平均每日销售量，我把数据放在下面的附件中：

问题1待聚类单品数据.xlsx

然后拿着这些数据我首先使用SPSS进行了系统聚类，这么做的目的是进行肘部法则从而有个大致判断，可以知道分成几类最好，下面是肘部法则的折线图：

对图片进行缩放拉伸，可以看到：

分成4类效果就很好，所以接下来使用SPSS做K-means聚类，可以得到聚类中心：

当然也可以得到每个单品的类别：

最后给这四个类别取个适合的类别名称即可。

结语

至此，论文C50的第一问的求解就复现完成了，希望这篇文章可以帮到大家，如有疑惑，欢迎留言！