数学建模实战快速入门

标签：实战 编码 入门 数据 建模 item 附件 单品

数学建模实战快速入门

---

前言：

数学建模竞赛是一个对参赛者的发散思维、数理逻辑能力以及资料搜集和处理问题的能力有着较高要求的综和性竞赛，希望获得好的奖项的参赛者需要对次提前进行充足的准备，包括提前找到最合适自己的小伙伴，提前了解和准备好各类需要的软件工具，并提前了解常见的算法及其适用范围，提前评估自己组建的队伍更擅长处理那一类的赛题，并提前确定好参赛期间的时间分配和组内分工。作为一个过去有过多次参赛经验的过来人，希望可以为参赛新人弥补信息差，并快速达到可以应对一定难度的建模竞赛赛题的程度。

关于几个常见的问题：

1.三个人组成的队伍需要明确的进行职责的划分吗？比如一个专门负责建模，一个专门负责编程求解，一个专门写论文？

答：一个参赛队伍存在职责的划分是有必要的，理想情况下，一个建模队伍里三个人都应该具备一定的建模竞赛所需的数学素养，并且三个人所擅长的部分有所区别（这个也是可以通过一定时间的训练来使自己可以胜任自己在队伍中的职责）并足够熟练，可以对应建模、变成、论文撰写三个分工。但往往很多队伍并不能从一开始就具备如此合适的人选，如果时间充足其实也可以考虑提前商量好分工，并且分别去为自己所负责的领域进行提前准备，不能错过的一个环节还是要一个队伍尝试去寻找一些题目（不需要太困难）来提前演练一下分工，不要在小问题上争吵太多，学会沟通甚至必要时候为了整体的进度学会妥协。

2.数学建模的负责不同职责的选手应当分别注意些什么？

熟悉数学建模个赛算法并且懂得去讲具体问题转化为数学模型的人适合负责队伍的建模工作，负责建模工作的人往往是最需要拥有灵活的逻辑思维且有强大的数学基础的人来负责，负责建模的人应当平时在思维上训练如何将显示问题转化为数学模型并且考虑模型是否具有可求解性，巧妙地模型远胜于复杂的模型。

负责编程求解的人应当提前花费时间和精力去了解常用的数学建模软件，比如最常用的MATLAB、Python；当然R语言和lingo也是可以使用的工具，如果是设计到统计分析建模的赛题，学习如何使用Eviews、spss等工具也是有必要的。有很多队伍往往是卡在缺少一个有能力实现和求解队伍所建立的模型的编程手因而极大限制了队伍的选题范围以及模型上限。

负责论文撰写选手的往往需要提前准备好一个排版严格符合要求并且在审美上比较整洁和美观的模版，这样会提前降低不必要的工作量，避免在比赛时手慌脚乱，写论文无处下手，并且负责撰写论文的选手也是也是需要把握整体的进度（相当于在学生时代体验一下当项目经理的感受），因为建模比赛最后的成果都需要论文撰写者来负责展示。当队伍确定了基本的选题、方向以及初步建立起了模型或者得出了一些初步的结论，论文撰写手就可以提前在论文中展示你们的成果，不要等到建模和求解工作完成的差不多的时候才开始动笔。（除非你对你们的模型和求解相当自信，认为可以提起前得到理想的求解结果，并且你习惯一气呵成的完成论文，并且对自己的写作速度以及排版熟练度相当自信）一般是问题分析、模型建立部分可以提前写；模型求解和结论部分等队友进度，论文摘要需要留下充足的时间在最后认真反复雕琢。

3.第一天比赛要不要趁热打铁赶进度？

一般是第一天的下午进行放题，在得到题目后第一天晚上的工作主要是确定好选题并且有一个初步的想法就大致完成任务了，第一天晚上我是建议所有的选手应当休息的，第一天不好好休息，后面的整体效率会不可避免的下降，这就导致后面也需要继续熬夜，得不偿失。

下面我们就进入正题了，先帮各位简单入门一下建模竞赛所需要的工具吧，论文和绘图的工具就不多做介绍了，在这里我会先以实际的赛题为案例简单介绍一下各种数学建模代码工具的实用用法。

---

具体参见下图为2023年高教社杯全国大学生数学建模竞赛的C题：

赛题分析：

---

鉴于近几年人工智能的兴起，学会善用先进的生产工具会给各项工作的效率带来一定的提升；当然人工智能并不是万能的，它实际上不能帮助你建立一个靠谱的模型；但可以帮助你免去大量搜集资料和整理建模思路的时间。现在我将赛题的全文输入到Kimi中，并要求它帮助我们分析一下赛题，Kimi给出的回答如下：

这个问题是一个典型的运筹学和统计学问题，涉及到需求预测、库存管理和定价策略。下面我将分别对每个问题进行分析，并给出可能的解决方案：

问题1: 分析蔬菜品类及单品销售量的分布规律及相互关系

• 数据分析：首先，需要对附件2中的历史销售数据进行分析，包括销售量、销售时间等。可以使用时间序列分析、相关性分析等方法来识别销售数据的趋势和周期性。
• 统计方法：应用统计模型，如ARIMA模型，来预测未来销售量。同时，使用聚类分析或主成分分析来识别不同品类和单品之间的关联性。

问题2: 分析销售总量与成本加成定价的关系，并给出补货和定价策略

• 成本加成定价模型：根据附件3中的批发价格，计算每个品类的成本加成定价。这通常涉及确定一个加成比例，以覆盖成本并实现利润。
• 补货策略：结合问题1中的销售预测，使用库存优化模型（如EOQ模型）来确定最优的补货量，以最大化利润并最小化库存成本。

问题3: 制定单品的补货计划

• 单品补货量：根据附件4中的损耗率和问题2中的销售预测，确定每个单品的补货量。需要满足最小陈列量的要求。
• 定价策略：对于每个单品，根据其销售趋势和损耗率，制定灵活的定价策略，以吸引消费者并减少库存积压。

问题4: 需要采集的相关数据及帮助

• 市场数据：收集消费者偏好、季节性因素、竞争对手定价等市场数据，以更好地理解市场需求。
• 供应链数据：了解供应商的供货稳定性和成本变化，以优化采购决策。
• 损耗数据：持续监控损耗率，以便及时调整补货和定价策略。

总结

这个问题需要综合运用统计分析、运筹学模型和市场研究来解决。实际的解决方案将依赖于具体的数据和商业环境。如果需要进一步的数学建模或编程实现，可以提供更详细的数据和要求。

我们继续追问要求给出详细的建模方案：

问题1: 销售量的分布规律及相互关系分析
数据预处理：清洗数据，处理缺失值和异常值。
描述性统计分析：了解各品类和单品的销售量、价格、时间等基本统计特性。
时间序列分析：对销售数据进行时间序列分析，识别趋势、季节性和周期性。
相关性分析：使用皮尔逊相关系数等方法分析不同品类和单品之间的相关性。
多元回归分析：建立多元回归模型，分析不同因素对销售量的影响。
问题2: 销售总量与成本加成定价的关系分析及补货策略
成本加成定价模型：设定定价公式，如
设定定价公式 , 如 P = C + m × C , 其中 P 是售价 , C 是成本 , m 是加成比例。 \text{设定定价公式},\text{如}P=C+m\times C, \\ \text{其中}P\text{是售价},C\text{是成本},m\text{是加成比例。} 设定定价公式,如P=C+m×C,其中P是售价,C是成本,m是加成比例。
需求预测：使用时间序列分析或机器学习方法预测未来一周的需求量。
库存优化：使用经济订货量(EOQ)模型或更复杂的库存模型来确定最优补货量。
利润最大化：建立目标函数，考虑成本、收入和库存持有成本，使用线性规划或非线性规划方法求解。
问题3: 单品补货计划
单品销售预测：对每个单品使用问题1中的方法进行销售预测。
陈列量约束：确保每个单品的补货量满足最小陈列量的要求。
补货量优化：结合预测销售量和损耗率，使用整数规划或混合整数规划来确定每个单品的补货量。
定价策略：根据单品的供需情况和竞争环境，制定动态定价策略。
问题4: 数据采集建议
消费者行为数据：收集消费者的购买历史、偏好和反馈，用于市场细分和个性化推荐。
市场趋势数据：监控宏观经济指标、行业趋势和季节性因素。
供应链数据：收集供应商的供货能力、价格波动和交货时间数据。
损耗和退货数据：持续跟踪损耗率和退货率，分析原因并优化管理流程。
实施步骤
数据收集与整合：收集所有必要的数据，并进行整合和清洗。
模型建立与校准：根据收集的数据建立数学模型，并通过历史数据进行校准。
模型验证与测试：使用交叉验证等方法测试模型的准确性和鲁棒性。
策略制定与实施：根据模型结果制定补货和定价策略，并在实际运营中实施。
监控与调整：持续监控实施效果，并根据反馈进行策略调整。
请注意，这只是一个高层次的建模方案，具体的模型和算法选择将取决于实际数据的特性和业务需求。此外，实施这些策略可能需要专业的统计软件和优化工具。

你当然也可以继续追问具体的代码案例，但问题是人工智能给出的方案太过宽泛，并且给出的回答有时候也不够准确，AI并不能给出一个具体的用于解决问题的方案，因此数学建模的各种技巧和建模本身的难度并没有因为人工智能的出现而大幅降低，反而因为人工智能使得信息扁平化，想要取得好的成绩，对模型建立的完整度和水平反而提出了更高的要求。

下面我要结合自己对于赛题的理解，以及人工智能给出的建议对赛题进行进一步的分析，并且形成属于自己的建模方案。

我们首先重新阅读和梳理题目，通过对关键词的提取和对逻辑的梳理得到下面的流程图关系图：

而题目对附件内容也进行了说明：

附件1给出了某商超经销的6个蔬菜品类的商品信息；

附件2和附件3分别给出了该 商超2020年7月1日至2023年6月30日各商品的销售流水明细与批发价格的相关数据；

附件 4 给出了各商品近期的损耗率数据。

现在我们开始去查看和分析具体的数据，并试图从中获得更多的信息以及切入点。

当然我们也可以直接以问题一为切入点去分析数据。

下图是附件2中的数据，在研究销量的分布情况，可以从销售日期这个时间维度来进行分析，具体可以分析各个单品或者品类在每日、每周、每月甚至以年为单位的销售分布。

使用Excel对数据进行一些初步的探索也是很方便的选择：

下面的部分，我会展示使用多种数学建模代码工具对C题的数据进行初步的一些数据分析和处理；在下一篇文章中，我会将C题的全部思路和代码公开，并且将主要使用更为主流和方便使用的Python语言。

Python代码

---

# 数学建模国赛的C题偏向数据分析和优化求解。这类赛题的原始数据一般数据量比较大，并且当你尝试使用MATLAB读取题目给出的excel文件时，你会发现字符串型的数据读入MATLAB时是无法使用MATLAB最为强大的矩阵去存储的，这时候使用python进行数据的处理和各种分析求解工作或许是更为方便的选项。


#2023年数学世界杯全国大学生数学建模竞赛C题
#加载数据分析所需的类库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']#设置中文字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号

from pandas import DataFrame,Series 

#读取数据
vegetables = pd.read_excel('附件1.xlsx')
details  = pd.read_excel('附件2.xlsx')
price = pd.read_excel('附件3.xlsx')
attrition_rate = pd.read_excel('附件4.xlsx')

# 数据预处理
print(np.any(vegetables.isnull()))  # 只要有一个空值便会返回True，否则返回False
print(np.any(details.isnull()))  # 只要有一个空值便会返回True，否则返回False  

#统计各个单品的销量总和
single_item=details['单品编码'].unique().tolist() #记录不重复的单品编码
df_singlesales=pd.DataFrame(columns=['单品编码','销量'])#创建新的空dataframe
for item in single_item:
    sales=details[details['单品编码']==item]['销量(千克)'].sum()
    df_singlesales.loc[len(df_singlesales)] = [item, sales]
print(df_singlesales)   

---

# 蔬菜品类与蔬菜单品之间的关联统计--统计各个品类的销售总量区别
vegetables = pd.read_excel('附件1.xlsx')
single_item=details['单品编码'].unique().tolist() #记录不重复的单品编码
type_item=vegetables['分类编码'].unique().tolist() #记录不重复的分类编码
#type_names=vegetable['分类名称'].unique().tolist()
typeitem_name = {}
typeitem_single = {}
for item in type_item:
    print(item)
    typeitem_name[item] = vegetables[vegetables['分类编码']==item] ['分类名称'].unique()    #字典可以不用一对一，也可以一个关键字对应一个列表，使用unique方法去除重复项
    print("---")                                     
    typeitem_single[item] = vegetables[vegetables['分类编码']==item]['单品编码']

---

type_sales=pd.DataFrame(columns=['分类名称','总销量','单品数量','平均销量'])#创建新的空dataframe

for item in type_item:
    
    sales = df_singlesales[df_singlesales['单品编码'].isin(typeitem_single[item])]['销量'].sum()

    count = len(typeitem_single[item])

    avg_sales = sales/len(typeitem_single[item])

    name = typeitem_name[item]
        
    type_sales.loc[len(type_sales)] = [name,sales,count,avg_sales]
 
print(type_sales)

#我们通过对每月的平均销量进行统计，进而分析不同品类蔬菜的销量与季节是否存在关联

df = details #附件2的数据
df['销售日期']=pd.to_datetime(df['销售日期'])
df['月份'] = df['销售日期'].dt.month 

month_sales=pd.DataFrame(columns=['月份','总销量'])#创建新的空dataframe
#item = 1011010101
item = 1011010201
#item = 1011010402
#item = 1011010501
#item = 1011010504
#item = 1011010801

items =  df[df['单品编码'].isin(typeitem_single[item])][['月份','销量(千克)']]
for i in range(1,13):
    #sales = df[df['月份']==i ]['销量(千克)'].sum()
    sales = items[items['月份']==i ]['销量(千克)'].sum()
    month_sales.loc[len(month_sales)] = [i,sales]

print(month_sales)

Excel数据处理 & MATLAB代码

---

使用MATLAB直接读取Excel中的数据，你也许会发现数据没办法读入，或者难以被存储到矩阵变量当中；实际上MATLAB作为数学求解软件，更偏向于读入数值型的数据，因此为了减少不必要的麻烦，我建议首先使用python或者Excel对数据进行初步的整理再读入MATLAB当中进行分析和计算。

在上述给出的python代码中你也许会觉得仅仅是计算各个品类的月销量就需要经过这么复杂的操作，如果对数据表本身进行修改，比如附件2中增加一列用以表示每个单品编码对应的品类编码或者说是品类名称，对于后续的分析与计算也许会方便非常多。那如何实现这个想法呢？

如果你有过学习oracle或者mysql的基础，你也许能想到在sql查询语句中或许可以通过表连接的方式将附件1和附件2中的内容拼接到一起。如果你想保留附件2的全部信息，那就可以以附件2为主表（左表），采用左外连接的方式实现我们的想法。

假设我们将附件1和附件2的数据导入到了oracle当中，将附件1的表命名为vegetable；将附件2的表命名为detail

左外连接的查询语句如下：

select v.*，d.单品名称，d.分类编码，d.分类名称 from vegetable v left join detail  d on v.分类编码 = d.分类编码；

如果你没有学过sql，python中也有类似的方法：merge，具体的参数意义或者其他使用方法可自行查阅。

import pandas as pd #使用merge方式实现数据拼接

df1 = pd.read_csv('file1.csv')
df2 = pd.read_csv('file2.csv')
result = pd.merge(df1, df2, on='key_column', how='left')#key_column是两个表都有的那一项，用作建立连接关系
result.to_csv('result.csv', index=False)

当然直接使用Excel进行这项操作也是可以的，这也是我要着重介绍的，Excel的功能实际上十分强大，善用Excel可以方便快捷的处理数据，而在数学建模竞赛中，设计数据分析和统计建模的任务，Excel也可以承担相当多的任务。

现在我们的目标是将附件1和附件2中的数据按照单品编号进行拼接。

我们面临的第一个问题是附件2中的单品编号是文本型数据，我们希望将其转化为数值型。可以在某一列，比如在H列编写公式：=VALUE(C2),回车后再下拉填充就OK了，但问题是附件2有88w多条数据，该如何处理呢？

可以使用ctrl + G 定位范围H2：H878504，再使用快捷键ctrl + D填充就OK了，其他的需要下拉填充的数据都可以这样操作。

excel中的VLOOKUP函数以及index+match操作都可以实现我们的目标；我们这两种方式可以都尝试一下：

添加单品名称列：=VLOOKUP($H2,附件1!A:B,2,)

添加分类编码列：=INDEX(附件1!C:C,MATCH(附件2!H2,附件1!A:A,0))

添加分类名称列：=INDEX(附件1!D:D,MATCH(附件2!H2,附件1!A:A,0))

完成拼接后就可以使用Excel的函数计算我们感兴趣的值，我们计算了不同品类的蔬菜在不同季度的平均销量，绘制折线图可以直观的了解到各个品类的销售量是否存在季节性变化，以及进几年的销量变化趋势。

使用MATLAB一直以来都是数学建模竞赛的主流选择之一，原因就在与MATLAB拥有全面且强大的数学函数库，方便的矩阵操作，大量专业的数学模型求解器以及相对简单易理解的代码风格和友好完善的帮助文档。（早年就流传“MATLAB除了不能生孩子，啥都能干”）

由于C题是属于数据分析类的赛题，使用MATLAB当然也可以实现我们的目标，与此同时这类题目使用python、spss、lingo、R甚至Excel都是不错的选择，因此MATLAB显得没有那么不可或缺。实际上，即使是数据分析类的赛题，如果当你需要建立一个优化模型并需要求解或者使用某些数学建模算法想要去寻找一些现成的代码作参考，你就会发现数学建模其实很难绕开MATLAB的使用。

所以在这里我也想初步的介绍一些有关MATLAB的实战内容。 想要在使用MATLAB中更进一步就需要个人多去尝试使用常见的操作并且在实际解决问题的过程中学会查阅帮助文档等相关学习资料。

我们使用execl的函数功能计算了6个品类在3年内的每日销量数据，并放在了附件2的mat_date工作表内，我们使用MATLAB将该数据读入并绘制折线图。

clear；
clcl；
data=xlsread('附件2.xlsx','mat_date','B2:G1096');
x=[1:1095]; 
x=x'; 
y1=data(:,1); 
y2=data(:,2); 
y3=data(:,3); 
y4=data(:,4); 
y5=data(:,5); 
y6=data(:,6); 
figure(1) 
hold on 
plot(x,y1,'b', 'LineWidth', 1) 
plot(x,y2,'r','LineWidth', 1) 
plot(x,y3,'k','LineWidth', 1) 
set(gca,'linewidth',1) 
set(gca,'Box', 'on') 
figure(2) 
hold on 
plot(x,y4,'g','LineWidth', 1) 
plot(x,y5,'c','LineWidth', 1) 
plot(x,y6,'m','LineWidth', 1) 
set(gca,'linewidth',1) 
set(gca,'Box', 'on') 

同样的，我们也可以对读入的各品类的日销量数据的统计量进行计算，以描述日销量的分布情况。

clear;
clc;
Test = xlsread('附件2.xlsx','mat_date','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]  %矩阵表示统计量

第一问中要求我们探讨单品或者品类之间的相关性；我们计算spearman相关系数来实现这一目标：

MATLAB代码如下：

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的等级 
    RX = tiedrank(X);
    RY = tiedrank(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 

clc;
clear;
RX=xlsread('附件2.xlsx','mat_date','B2:G1096');
[ R , P ]= fun_spearman(RX, 1)  

% 绘制热力图
figure; % 创建一个新的图形窗口
imagesc(R); % 使用imagesc绘制相关系数矩阵
colorbar; % 显示颜色条

% 设置颜色映射
colormap('jet'); % 可以使用'hot', 'cool', 'spring', 'summer', 'autumn', 'winter', 'gray', 'bone', 'copper', 'pink'等

% 添加标题和坐标轴标签
title('六大品类的Spearman相关系数');
%xlabel('Variables');
%ylabel('Variables');

% 添加网格线
grid on;

% 调整坐标轴刻度，使其与矩阵的行和列对齐
axis square;
set(gca, 'XTick', 1:size(R, 2), 'YTick', 1:size(R, 1));

% 在每个单元格中显示数值
for i = 1:size(R, 1)
    for j = 1:size(R, 2)
        text(j, i, num2str(R(i, j), '%.3f'), ...
             'HorizontalAlignment', 'center', ...
             'VerticalAlignment', 'middle', ...
             'Color', 'black'); % 显示数值，保留三位小数
    end
end

计算结果如下图：

R代码

---

R语言我就不多赘述了，作为统计与数据分析的开源软件，R语言拥有很多功能强大的包，包括其灵活的语言，美观的作图，能给我们提供很大的帮助，不过想要深入学习R语言也是需要耗费不少精力的。下面就粘贴了一份某优秀论文的R语言第一问代码帮读者感受下R语言的特点

R语言一般推荐使用Rstudio进行文件的编写和管理，在左上部分的窗口编写代码文件后，可以逐行执行，也可以执行整个文件（个人喜欢逐行执行）。

个人的感受是R语言虽然存在大量强大的包，但由于是开源免费的平台，且不像python那样如此受欢迎，因此大量的包是处于无人维护的状态，版本更新后一些包可能就无法使用,而且R语言的语法也与主流的语言在语法上有一些区别，个人建议对R语言不够熟悉的同学在比赛过程中就不要去尝试使用R语言了，以避免出现一些问题，导致代码无法运行下去，当然如果本身熟悉R且经常使用R，R也可以是一款不逊色Python的编程语言。

library("readxl") 
library("forecast") 
a = read_excel(path = "附件1.xlsx", sheet = "Sheet1") 
b = read_excel(path = "附件2.xlsx", sheet = "Sheet1") 
c = read_excel(path = "附件3.xlsx", sheet = "Sheet1") 
d = read_excel(path = "附件4.xlsx", sheet = "Sheet1") 
correspond=a[,c(1,3)] 
library(showtext) 
showtext_auto() 

library(dplyr) 
library("data.table") 
correspond=data.table(correspond) 
setkey(correspond,单品编码) 
bt=data.table(b) 
setkey(bt,单品编码) 
b_1=bt[correspond] 
b_2=aggregate((`销量(千克)`)~(分类编码+销售日期),data=b_1,FUN = sum) 
b_2_11=b_2[b_2$分类编码=="1011010101",] #之后会对比改进前的时间序列分析 

## 对每个单品的异常值处理(数据预处理) 
mylist=unique(b_1$单品编码) #获得所有蔬菜单品的编码 
b_i=b_1[b_1$单品编码==mylist[1],] 
b_i_1=aggregate((`销量(千克)`)~(单品编码+销售日期),data=b_i,FUN = sum) 
#获得第一个蔬菜的每天的售量之和 
remove_outlier=function(x){ 
  names(x)[3]="a" 
  q1=quantile(x$a)[[2]] 
  q2=quantile(x$a)[[3]] 
  q3=quantile(x$a)[[4]] 
  iqr=q3-q1 
  minq=q1-1.5*iqr 
  maxq=q3+1.5*iqr 
  x$a[x$a>maxq|x$a<minq]=NA 
  x$a[is.na(x$a)]=mean(x$a,na.rm=TRUE) 
  names(x)[3]="销量(千克)" 
  return(x) 
} 
b_i_1=remove_outlier(b_i_1) 
bb=b_i_1 
for(ii in 2:251){ 
  b_i=b_1[b_1$单品编码==mylist[ii],] 
  if(nrow(b_i)!=1){ 
    b_i=aggregate((`销量(千克)`)~(单品编码+销售日期),data=b_i,FUN = sum) 
    b_i=remove_outlier(b_i) 
    bb=rbind(bb,b_i) 
  }else{ 
    b_i=b_i[,c(3,1,4)] 
    bb=rbind(bb,b_i) 
  } 
} 
b_i11=aggregate((`销量(千克)`)~(单品编码+销售日期),data=b_1,FUN = sum) 
rm=bb 
rm=na.omit(rm) 
#write.csv(rm,"remove.outlier.csv") 

## 改进后的时间序列分析 
rmt=data.table(rm) 
setkey(rmt,单品编码) 
b_1=rmt[correspond] 
b_2=aggregate((`销量(千克)`)~(分类编码+销售日期),data=b_1,FUN = sum) 

b_2_1=b_2[b_2$分类编码=="1011010101",] 
b_2_2=b_2[b_2$分类编码=="1011010201",] 
b_2_3=b_2[b_2$分类编码=="1011010402",] 
b_2_4=b_2[b_2$分类编码=="1011010501",] 
b_2_5=b_2[b_2$分类编码=="1011010504",] 
b_2_6=b_2[b_2$分类编码=="1011010801",] 

sxfx=function(x){ 
  b_sxfx=x 
  b_sxfx$销售日期1=as.Date(b_sxfx$销售日期) 
  sales <- b_sxfx[, 3] 
  dates <- b_sxfx[, 4] 
  dates <- as.Date(dates) 
  weekly_sales <- data.frame(Week = format(dates, "%Y-%W"), Sales = 
                               sales) %>% 
    group_by(Week) %>% 
    summarise(TotalSales = sum(Sales)) 
  #按周进行加总 再进行ARIMA分析 
  myts <- ts(weekly_sales[, 2], frequency = 52) 
  #plot(myts) 
  #acf(myts) 
  #pacf(myts) 
  fit <- auto.arima(myts,seasonal=TRUE) 
  myts <- ts(weekly_sales[, 2], frequency = 1) 
  plot(myts,xlab="Week",ylab="Total sales per week") 
  summary(fit) 
} 
sxfx(b_2_1) 
sxfx(b_2_11) 
sxfx(b_2_2) 
sxfx(b_2_3) 
sxfx(b_2_4) 
sxfx(b_2_5) 
sxfx(b_2_6) 

## 品类、按日、总销量算corr并画热力图 
b21=b_2_1$'(\'销量(千克)\')'
b22=b_2_2$'(\'销量(千克)\')'
b23=b_2_3$'(\'销量(千克)\')'
b24=b_2_4$'(\'销量(千克)\')'
b26=b_2_6$'(\'销量(千克)\')'

length(b21) 
length(b22) 
b22=c(b22,0) 
length(b23) 
length(b24) 
b24=c(b24,rep(0,35)) 
length(b25) 
length(b26) 

bbb=data.frame(花叶类=b21, 
                  花菜类=b22, 
                  水生根茎类=b23, 
                  茄类=b24, 
                  辣椒类=b25, 
                  食用菌=b26) 
cor(bbb) 
library(corrplot) # 基础可视化-定量变量（相关性热图） 
bbb %>% 
  select_if(is.numeric) %>% 
  cor() %>% 
  corrplot(method = 'color', order = 'AOE', addCoef.col = 'black') 

## 品类、按月、总销量corr并画热力图 
b_2$月份 <- format(b_2$销售日期, "%Y-%m") 
monthly_sales <- aggregate('(\'销量(千克)\')' ~ (分类编码+月份), b_2, sum) 
bb_2=monthly_sales 
bb_2_1=bb_2[bb_2$分类编码=="1011010101",] 
bb_2_2=bb_2[bb_2$分类编码=="1011010201",] 
bb_2_3=bb_2[bb_2$分类编码=="1011010402",] 
bb_2_4=bb_2[bb_2$分类编码=="1011010501",] 
bb_2_5=bb_2[bb_2$分类编码=="1011010504",] 
bb_2_6=bb_2[bb_2$分类编码=="1011010801",] 

bb21=bb_2_1$'(\'销量(千克)\')' 
bb22=bb_2_2$'(\'销量(千克)\')'
bb23=bb_2_3$'(\'销量(千克)\')'
bb24=bb_2_4$'(\'销量(千克)\')'
bb25=bb_2_5$'(\'销量(千克)\')' 
bb26=bb_2_6$'(\'销量(千克)\')'

length(bb21) 
length(bb22) 
#b22=c(bb22,0) 
length(bb23) 
length(bb24) 
#b24=c(bb24,rep(0,35)) 
length(bb25) 
length(bb26) 

bbbb=data.frame(花叶类=bb21, 
                   花菜类=bb22, 
                   水生根茎类=bb23, 
                   茄类=bb24, 
                   辣椒类=bb25, 
                   食用菌=bb26) 
cor(bbbb) 
library(corrplot) # 基础可视化-定量变量（相关性热图） 
bbbb %>% 
  select_if(is.numeric) %>% 
  cor() %>% 
  corrplot(method = 'color', order = 'AOE', addCoef.col = 'black') 

## 不同单品蔬菜的销量计算corr并画热力图 
#单品：将总销量占比低于1%的行去除 
kk=aggregate((`销量(千克)`)~(单品编码),data=b_1,FUN = sum) 
kk$比例=kk$'(\'销量(千克)\')'/sum(kk$'(\'销量(千克)\')') 
b_1_75_list=kk$单品编码[kk$比例<=0.01] 
for(i in 1:length(b_1_75_list)){ 
  if(i==1){ 
    b_1_75=b_1[b_1$单品编码!=b_1_75_list[i],] 
  }else{ 
    b_1_75=b_1_75[b_1_75$单品编码!=b_1_75_list[i],] 
  } 
}   
b_1_75=na.omit(b_1_75) 
#获得不同单品蔬菜的月销量corr 
b_1_75$月份 <- format(b_1_75$销售日期, "%Y-%m") 
b_1_75=data.table(b_1_75) 
setkey(b_1_75,单品编码) 
at=data.table(a) 
setkey(at,单品编码) 
b_1_75=b_1_75[at] 
b_1_75=na.omit(b_1_75) 
monthly_sales <- aggregate(`销量(千克)` ~ (单品名称+月份), b_1_75, sum) 
library(tidyr) 
b_1_75_1 <- pivot_wider(monthly_sales, names_from = 单品名称, values_from = `
                        销量(千克)`) 
for(i in 1:nrow(b_1_75_1)){ 
  b_1_75_1[i,][is.na(b_1_75_1[i,])]=0 
} 
b_1_75_1=b_1_75_1[,-1] 
cor(b_1_75_1) 
# 使用CURE聚类算法对数据进行聚类，生成聚类标签 
library(cluster) 
set.seed(123)  # 设置随机种子以获得可重复的结果 
k <- 3  # 聚类数 
cluster_labels <- pam(cor(b_1_75_1), k) 
# 将聚类标签添加到原始数据中 
data_a <- cbind(cor(b_1_75_1), Cluster = cluster_labels$clustering) 
data_a=data.frame(data_a) 
# 根据聚类标签重新排序数据 
data_sorted <- data_a[order(data_a$Cluster), ] 
# 绘制热力图（使用heatmap函数或其他绘图方法） 
heatmap_data <- data_sorted[, -ncol(data_sorted)]  # 去除聚类标签列 
for(i in nrow(heatmap_data)){ 
  for(j in nrow(heatmap_data)){ 
    heatmap_data[i,j]=as.numeric(heatmap_data[i,j]) 
  } 
} 
heatmap(as.matrix(heatmap_data),  
        #labRow = NA,  # 隐藏行标签 
        #labCol = NA,  # 隐藏列标签 
        col = colorRampPalette(c("blue", "white", "red"))(100),  # 设置颜色范围 
        main = "") 

我会在3~5天内更新2023年高教社杯C题后续问题的完整思路和代码，并且一直到今年9月份的数学建模国赛开始前，我都会持续更新数学建模相关的内容，全部资料全部都会无条件免费分享，希望多多关注收藏！

标签：实战,编码,入门,数据,建模,item,附件,单品
From： https://blog.csdn.net/m0_47498690/article/details/140861451