数学建模——食品质量安全抽检数据分析

1、承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的如果引用别人的成果或其他公开的资料（包括网上查到的资料） ，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B 中选择一项填写）：*队 A 题我们的参赛队员为：1、姓名：*学院化学化工学号2012*2、姓名：*学院化学化

2、工学号 2012*3、姓名：日期：2013年5月1日评阅编号（由组委会评阅前进行编号）：评阅记录（可供评阅时使用）：评阅人评 分食品质量安全抽检数据分析摘要食品质量与安全是一个专业性很强的问题，其标准的制定和抽样检测及评价都需要科学有效的方法。 食品的质量和卫生问题是关系到民生的大问题，所以本文结合当今社会实际以及文章提出的不同问题，在阅读大量相关文献后，应用AHP 方法、方差分析方法、回归分析原理、分层抽样和线性目标规化等方法，建立相应的数学模型。利用TableCurve、MATLAB 软件工具进行求解， 对文章提出的相关问题有一个比较精确合理的解答。具体解答如下 :对于问题一： 首先将蔬菜

3、、 鱼类、鸡鸭等主要食品的安全影响因素按主要成分分析法划分为三大类：生物性污染（微生物等） 、化学性污染（食品添加剂等） 、物理性污染（重金属等）。对深圳市这 3 年 3 大安全影响因素的安全情况的变化趋势做出定量的综合评价可以转化为利用 AHP 法分析决定食品安全单位指标在综合指标中的权重的变化情况，可以得出食品安全状况明显改善。如图一：食品安全系数的等级在稳步提升。3.5食品安全系数3等级2.521.510.50201020112012图一对于问题二：对食品质量影响的因素如食品产地，食品销售地点和季节性等的关联度做一个比较， 以找到其中的规律。 为此我们经过商量决定采用数学上通用的方差以及

4、回归统计分析方法， 建立多元线性回归的数学模型，以找到一个随机变量与一组变量的相关性。同时运用有关灰色系统模型的相关理论。利用递推补集的思想，用总的m个变量的回归平方和减去（m-1）个变量的回归平方和来求出变量xi 的偏回归平方和1Qi，并且利用 TableCurve 软件进行多项式拟合，找到精确的回归方程。得到xi 的偏回归平方和如图二所示。图二所得结果为：与食品质量关联度最大的是食品产地， 其次是食品加工，季节性影响最小，食品抽查地点几乎无影响。对于问题三： 食品质量的安全和卫生在现实社会中起着非常重要的作用 , 食品进行抽检也需要一定人力、物力和财力 ( 即成本费用 ), 抽检的越多检测

5、效果就越好，但需要的时间就越长 , 其成本费用也就越高。 通过建立完全随机抽样方案 , 对某一代表食品的每一品牌的各批次各检测项目进行抽样检验。 以期望和方差等相应比较来检测模型可靠性。借助统计学方法 , 用样本均值与配方标准值的差距来刻画样本的代表性。用样本方差来衡量抽检方案的随机性 , 两者合为检测可靠性分析。 对模型进行修正 , 既考虑各个因素之间的相互影响 , 也要准确地体现出反馈效应对模型系统的影响。同时以蒙特卡罗法1n 一样对抽检的全过程进行模拟，我们能够期望观测的平均值趋于零的变化速率与快，以至于使 n时. Y1 Y2YnEY 。 还可以进行灵敏度测验。n关键词： AHP分析法；

6、多元函数回归模型；方差分析法；分层划分抽样；随机抽样；蒙特卡罗法。2一 问题重述“民以食为天”，食品安全关系到千家万户的生活与健康。随着人们对生活质量的追求和安全意思的提高， 食品安全已成为社会关注的热点， 也是政府民生工程的一个主题。城市食品的来源越来越广泛， 人们消费加工好的食品的比例也越来越高， 因此除食材的生产收获外，食品的运输、加工、包装、贮存、销售以及餐饮等每一个环节皆可能影响食品的质量与安全。 另一方面，食品质量与安全又是一个专业性很强的问题， 其标准的制定和抽样检测及评价都需要科学有效的方法。深圳是食品抽检、监督最统一、最规范、最公开的城市之一。请下载 2010 年、 2011

7、 年和 2012 年深圳市的食品抽检数据（注意蔬菜、鱼类、鸡鸭等抽检数据的获取） ，并根据这些资料来讨论：1. 如何评价深圳市这三年各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势；2. 从这些数据中能否找出某些规律性的东西：如食品产地与食品质量的关系；食品销售地点（即抽检地点）与食品质量的关系；季节因素等等；3. 能否改进食品抽检的办法， 使之更科学更有效地反映食品质量状况且不过分增加监管成本（食品抽检是需要费用的） ，例如对于抽检结果稳定且抽检频次过高的食品领域该作怎样的调整？3二、问题分析在当今时代， 食品安全情况越来越让人堪忧，频发的食品安全事件让社会出现了不和谐的音符。为

8、此进行相关的研究十分有意义的。1 、食品安全生产涉及许多过程，我们关心的是食品安全的影响因素，为此我们运用统计学分析中的方差确定出哪些因素是主要的， 哪些是次要的。 利用层次分析法分析不同因素在对食品安全影响不同的权重变化得知其内在的规律。 同时也可以利用方差分析原理到找对食品安全系数有显著影响的因素。 在对 2010 年、2011 年、2012 年、2013年的抽样数据的分析基础上， 考虑到影响食品安全系数的环节和抽样方法带来的影响因子，运用层次分析法和方差分析法判断出评价深圳市这三年各主要食品领域微生物、 重金属、添加剂含量等安全情况的变化趋势。2、在对原始数据进行处理， 剔除其中的异常点

9、。 建立多元函数的线性回归方程，利用最小二乘法估计参数。 把模型写成矩阵的形式， 化简整理得其正规方程组参数。 通过对正规方程组的求解， 最后得到回归方程。 遵照建立回归方程的程序进行回归方程显著性检验，问题转化为建立的待检假设：H0 :0 ，若能通过检验拒12n绝 H 0 ,则 Y 与 m 个变量 xi (i 1,2,3m) 之间存在线性相关关系。3、为了建立合理的抽检模型，我们在对现有的统计信息进行定量分析，确定主要的食品和抽检项目，利用分层抽样模型确定这些主要食品中每种品牌所要抽检批次的最优值，根据实际情况计算出每种品牌中要抽的批次数以及每个批次抽检的项目数，对其优劣进行综合评价。三、问

10、题假设1、假设影响食品安全性因素仅分为三大类，其它没有被分类的因素对食品安全性所造成的影响忽略不计。2、假设调查样本在一定的范围内是均一的。3、假设所有食品生产厂商的信誉度均相同。4、环境中自然因素土壤对食品取材状况基本相同，其它因素对食品安全系数的生成没有影响。45、除食材的生产收获外，食品的运输、加工、包装、贮存、销售以及餐饮等每一个环节对食品的安全影响程度相同。6、食品的抽查过程对食品的安全系数影响程度相同。7、假设抽检不受国家相关政策的影响。8、每次抽检时都公平对待每一食品品牌，每一批次，且由计算机随机选出，不存在人为的干扰。9、假设模型求解过程中所用的数据都是合理的，且每次检测得结果

11、真实可靠。四、参数及符号说明O：目标层；C：准则层；P：方案层；aij ：准则层两个因素Ci 和 Cj 对目标层的影响程度之比；C.I.:判断矩阵一致性指标 ;bij:方案层两个因素Pi 和 Pj 对准则层的影响程度之比;W，W1，W2Wn :权重向量 ;max:最大特征值；R.I.:平均随机一致性指标 ;PI; 食品安全质量综合评价指数；C.R.: 一致性比例指标；Hi ：第 i 年食品不合格数；Ri: 第 i 年需要处理的不合格数；Xi ：影响食品安全的因素（i=1,2.8）;Y : 各食品安全系数（ 1,2,.,38） ;S : 离差平方总和；TS : 组内平方和；ES : 组间平方和；

12、A5：待估计系数参数（ 1,2,.,8）；iy?： Y 的回归函数值；i ：参数的最小二乘估计（ 1,2,.,8）；bi五、模型的分析、建立与求解将因素 A 分成 r 个水平，因素 B分成 s 个水平，对 A,B 的每一个水平的一对组合Ai ,B j 只进行一次试验，列出实验结果记录在下表中。因素 AB1B2Bs因素 BA1X11X 12X1 sAX 21X 22X 2S2ArX r 1X r 2X r s列 平 均值X ?1X ?2X?sX ? j行 平 均值X i ?X 1?X 2?X r ?X其中，Xij 表示因素 A 第 i 个视频与因素 B 第 j 个水平的一对组合Ai，Bj 进行一

13、次实验的结果。6记 n=rs , 则1均值 : Xn平方和：rs1X ij , X i ?i 1j 1ss1X ij , Xj 1? jrrX iji1rsSTX ijXi1j 1rsSEX ijXi1j 1r22i ?X ? jX2s2SA S X i ? X ,SB rX ? j Xi 1j 1可以得到平方和关系 :ST SE SA SB .(1)判断因素 A 的影响是否显著，就是要检验假设H 0 A： 1 j2 jrj? j , j1,2, s选取统计量 :AAE, 得到拒绝域Fs 1 S / S F r 1, (r 1)( s 1)W(s1)SAF r1(r1)( s1)FASE.(2

14、)若 FAF ,拒绝 H0A;否则，接收 H0A判断因素 B 的影响是否显著，就是要检验假设H 0 B： i 1i 2iri ? ,i1,2, r选取统计量 :FB( s 1)SB / SE F s 1, (r 1)( s 1) ,得到拒绝域7s1 SBF (r 1,(r 1)( s 1)(3)WFBSE若FBF ,拒绝 H0B;否则，接收 H0B.对影响食品质量的三大因素两两进行比较共六次后，根据其影响因素程度的变化得出食品的质量综合质量评价。有时如果要考察两个因素A,B 之间是否存在交互作用的影响，需要对两个因素各种水平的组合进行重复试验，以得到比较精确的解答。问题一的模型1、将食品安全综

15、合评价分解为3 个层次，最上层为目标层O，中间层为因素层 C，最下层为方案层 P.O:食品安全C:微生物重 金食品 添环 境其它影污染属加剂污染响P:2012 年2010 年2011 年2、针对问题一，我们建立如下的成对比较矩阵。假设要比较某一层 n 个因素 C1，C2， ,Cn对上一层一个因素 O的影响，每次取两个因素 Ci,Cj, 用 aij表示 Ci 和 Cj 对 O的影响之比，全部结果可以用成对比较a11a12a1n矩阵a21a22a2nA=an1an2ann8aij>0,aji=1/aij.表示，其中 aij=1,我们称此矩阵为正互反矩阵。依据深圳市 2010-2012年第一期

16、到第五期统计数据得知1143322175511111A= 4723112113511311353 、一致矩阵检验及权重向量的的确定。易知成对比较矩阵通常不是一致矩阵，应先转化为一致矩阵，为了用A 的最大特征根（记作max）的特征向量（归一化后）作为权向量W,即全向量满足AWmax W .(4)对于给出的 A 可以计算出max =5.073.同时其不一致性程度衡量需要用到：一致性指标 C.I.=maxn （）n=5.073-5/1.12=0.018;1C.I .一致性比率 C.R.= R.I . =0.018/1.12=0.016其中 R.I.为随机一致性指标。其数值如下表：n123456789

17、1011R.I.000.580.91.121.241.321.411.451.491.51其一致性比率 C.R.<0.1 时，认为 A 的不一致程度在容许的范围之内， 可用其归一化的特征向量TW= 0.263,0.475,0.055，0.099,0.110作为权向量。4、组合权向量对于 3 个层次的决策问题，记第2 层对第 1 层的权向量为W 2（1 2),（22), , n( 2 )T第 3 层对第 2 层的权向量为Twk3wk31,wk32, ,wkm3 , k 1,2, , n.( 3 )以k为列向量构成矩阵9W（ 3 )w1( 3 ) , w2( 3 ) , , wn( 3 )

18、,则第 3 层对第 1 层的组合权向量为w( 3 )W ( 3 ) w( 3 ) , .(5)把第 2 层对第 1 层的权向量记作w( 2 )w1( 2 ) , w2( 2 ) , , wn( 2 )T(0.263,0.475,0.055,0.099,0.110) T用同样的方法构造第三层对第二层的每一个因素的成对比阵，不妨设它们为125113B1112，1132B3111111523311113438B2311B41118333111141114B51114441,5 中的元素 bi jk其中矩阵 Bkk 1,2,是方案（食品安全） Pi 与 Pj 对于因素 Ck（微生物污染等）的优越性的比

19、较尺度。1,2, ,5 计算出权向量 wk3由第三层的成对比较矩阵 Bk k，最大特征根k 和一致性检验 CRk 列表如下：k123450.5950.0820.4290.6330.16630.2770.2360.4290.1930.166wk0.,1290.6820.1420.1750.688k3.0053.00233.0093CRk0.0050.00100.0080有列表可以看出成对比较矩阵B k 1,2, ,5 均通过一致性检验。k可计算组合权向量10w( 3 )W ( 3 ) w( 3 ) , =0.2630.5950.0820.4290.6330.1660.475T0.2770.236

20、0.4290.1930.166? 0.0550.300，0.246，0.456 .0.1290.6820.1420.1750.6680.0990.1105、组合一致性检验即逐层进行。 定义最下层（第 S 层）对第一层的组利用还原思想进行组合一致性检验，S*P合一致性比率为 CRP2CR ，.(6)其中第 P 层的组合一致性比率为CRP CIRIPP,P3，4， ,S.3332可以算出 CI0.00176, RI0.58, CR0.003,再加上已经有了的 CR 0.016.SP*所以CR P2CR =0.019 0.1，组合一致性通过，前面得到的组合权向量w( 3 )W ( 3 ) w( 3

21、) , 可以作为最终决策的依据。同理可得其它年份的组合权向量，根据影响食品安全系数的因素的权重的变化得出对深圳市这三年来的食品安全综合评价 （此时近似的以组合权向量代替） 。最后定量的表示出来深圳市这三年食品安全状况趋于良好。（二）问题二：“多元线性回归模型”1、食品质量这一随机变量与食品产地，抽检地点，季节等变量之间的相关性问题很自然地想到运用多元函数回归分析的方法来解决。2、设随机变量 Y 与 m个变量 x1, x2 , x3 , , xm 有关系Y01 x12 x2m xm + ·········&

22、#183;·················（7）2其中，为随机项，且i N(0， ) ,y11x1 1x1 m11Yy2, X1x2 1x2 m2,2记，yn1xn1xm nnm11则随机变量 Y 与变量 X 的关系可化为YX3、求回归系数nyi残差平方和 Qi 101 xi 1im xim 2求 Q分别关于， ， ，的一阶偏导数，并令它们等于0，得01m1X T X1T2X Y ···

23、3;········································（8）m利用 MATLAB求得：? =110.4733， ?=0.0826 ， ? =0.0478， ? =0.05

24、28， ? =0.1199， ? =-0.0257 ，012345? =0.1217 ， ? =0.1220， ? =-0.0015.6784、多元线性回归方程的显著性检验y1nnn2yi , syyyiy ,QSyyUi 1i1其中 U称为回归平方和。选取统计量U / mFQ / nm 1 (9)在 H 0 成立的条件下， FF（m-1,n-m-1) 。得到拒绝域WFU / mF (m 1, n m 1)Q / n m 1若 F>F 拒绝 H0即Y与m个变量 X1,X2 ,X m 之间存在线性相关关系； 否则，接受 H0,12即 Y 与 m 个变量 X 1 , X 2 ,Xm 之间不存

25、在线性相关关系。对于食品安全这个模型取显著性水平=0.05，可查表得到Fm1, nm1 =2.4047.由统计数据得到 F（m-1,n-m-1)=5.3966 2.4047 ，所以此模型是显著的。即Y 与 m 个变量 X 1 , X 2 ,X m 之间存在线性相关关系。同样我们也知道在多元线性回归模型中，拒绝假设H 0 ，即回归方程显著。我们还关心 Y 对 X 1 , X 2 ,X m 的回归哪些是主要因素哪些是次要因素， 需要采用偏 F 检验法（求解略）。5、 回归系数的显著性检验：利用一开始的建模部分的方差分析法由（1）式知在模型取显著性水平=0.05 时，可查表得到 Fm 1, n m

26、1 =3.0088.同样由统计数据可以得出F（m-1,n-m-1)=0.2623<3.0088,由此说明模型的拟合是好的，即模型的省略项造成的影响不大。6、还可以进行相关的预测（计算过程略）一个自然的想法是用预测量:y001 x102 x20m xm0来代替，与测量y0 的优劣取决于y0Y0 的大小 。记1 nnxkj xj ,i , j 1,2, , mxjxij ,lijxki xin i1k 1l1 1l1ml1' 1l1' mL, L 1lm1lm mlm' 1lm' md 211nmmQlij' x0 i xi x0 j xj , 21可

27、以证i 1j 1nm13明当 Y0与 Y1 ,Y2 ,Yn 相互独立时，在显著性水平下可得到 Y0 的预测区间 :y0 t0 n m 1d ,y0 t0 n m 1d······················（10）根据统计的相关数据可以定量的计算出在显著性水平=0.05 下即有 95%的可能性得到 Y0 的预测区间为 :（ 24.75,37.00 ）7、基于灰色关联分析的综

28、合评价， 即将评价问题转化为关联度， 按照一定的准则定义母序列，计算各个子序列与母序列的关联度， 再加以比较得到综合评价结果。经过选定母序列和比较序列， 计算关联度，关联度表示各因素评价值序列和和最理想的评价值序列相关联的程度，其数值越大表示此因素对食品质量的影响越大。通过MATLAB 对综合评价模型的求解 ,可得如下表xix1x2x3x4x5x6x7x80.7920.6800.5730.8500.9250.7900.4590.029综合评价值排名35621478结论：由上表可以得出食品产地影响食品质量最大， 食品加工次之， 再者是季节因素，最后的是抽查地点，影响比较小。8 、由于线性回归常数

29、项的影响存在，在上面处理的过程中是忽略了其影响，为了是数学模型更加精确，可以进行多项式拟合分析，将常数项的影响转化为Xi 变量的影响，这样能更好的拟合实际情况。所以建模如下：nf ( xi )yn2要使 J1,2,mmini 1得值最小，只需利用极值的J0(k0, m) ，得到关于必要条件， ， ， 的线性方程组01mKnmyii 1r 1xik 1k r kxi0nr mximk r kxiyii 1k 10线性最小二乘拟合可由解超定方程TTRRARY来求得。其中14r 1x1r 2 x1r m x1Rr 1x2r 2x2r m x2。r 1xnr 2xnr m xn当 r 1 x ,r 2

30、x , r m x 线性无关时，以下利用 TableCurve进行多项式的拟合。图三TA1, 2,my1,y2,TY, ynTR R 可逆，方程组有唯一解。（三）问题三、模型的模拟与检验1 、要改进食品抽检的办法， 使之更科学更有效地反映食品质量状况且不过分增加监管成本，可以进行分层次划分抽样。假设抽检不合格几率为 (1-P) 。第一种方案是：抽检费用为 C元 / 次。第二种方案是 : 按组抽查 x 次为（ B+x）元, 一发现不合格的就要重新以第一种抽查，所以可以建立如下模型。总费用的期望xEFB(C1)xcxp ······

31、83;·······（ 11）15则每一次抽取的平均检验费用为：BC1-CPXFX.(12)利用 TableCurve 软件数形结合。（1）、若 P=0.997， B=60， C=89;则（ 11) 式的平面图形如下，图四由图可以清楚的看到存在最优解 X，使得抽查过程中的 每一次抽取的平均检验费用 F 最小，最小 X 及最小 F 同样可以通过 TableCurve 软件求得，在图四中可以大致看出所求得最小F=9, 同理由图五可以求得最小X=17，同时也可以进行灵敏度的测试。即sF , q dF ? q dq F=0.16 ，

32、于是 q 的微小变化可能不会导致检验费用大的变化。更一般的稳健性分析要考虑独立性的假设。图五16图六（2）、函数 Z60XX90 89 Y 的三维立体图形如下：这是在抽检的合格几率为变量Y 下所得。在实际问题中我们可以清楚的知道函数 f （x,y ）的最大值或最小值一定在 D 内取得，而函数在 D内只有一个驻点， 那么可以肯定函数在驻点处的函数值就是最大值或最小值，对于多元函数的最值求法可以由下面的定理求得：Z0, Z0YX用 TableCurve 软件做出相应的图形并且求得其解。即Z60xX289 y ln y 0,XZxy89 xy 017图七图八18图九2、最后用蒙特卡洛方法进行模拟与试

33、验，具体步骤如下：初始化 :count=0.步骤一：产生随机数；步骤二：判断随机数是否在指定的区域内，若在，则计数器加一，否则若 count N， count+, 转步骤一，否则转步骤三。步骤三：利用数据数据进行相关统计分析,得出实际问题的解答。四模型的评价4、 1模型的优点1 、在问题一中通过层次分析法进行相关的计算和决策，以比较精确的结果综合评价了深圳市近三年来主要食品领域的变化趋势， 在运用层次分析法中进行了组合权向*SP量一致性检验， CRP 2 CR =0.0190.1使得一致性检验结果通过，得到比较精确的组合权向量，对文章提出的问题有了一个在客观上定量的解答。2 、在问题二中采用多元线性回归的方法考虑影响食品质量的因素，通过对模型进行显著性与拟合性检验可知， 回归方程是显著的， 拟合性也是较好的， 因此该模型比较准确的刻画了食品质量与食品产地， 季节性因素等之间关系。 同时还使用了灰色模型理论的关联度进行食品质量与食品产地， 季节性因素等之间关系的刻画， 得到的结果用来作为食品综合评价的参考物， 是问题具有可比性的定量衡量标准。 同时还对影响食品质量的随机因素转化到对 X 变量