食品质量安全抽检数据分析的模型探究.docx

.数学建模校内竞赛论文 论文题目：食品质量安全抽检数据分析的模型探究 组号：68#成员： 选题：A题 2013-06-05摘要 对于问题1，针对该题我们大量的收集和整理了2010年到2012年深圳市的食品安全情况抽检数据；为了能综合考虑生产领域、流通领域、餐饮领域以及在微生物、添加剂、重金属方面对3年来食品安全的变化趋势，我们采用了层次分析法：将微生物、添加剂、重金属作为第一准则层，将生产领域、流通领域、餐饮领域作为第二准则层，将2010年、2011年、2012年作为方案层。最终得到了3年中，2011年的食品安全问题最为严重，2010年和2012年与2011年相比较为缓和，2012相对于2010年、2011年的食品安全问题有所改善。对于问题2中的食品质量与食品产地的关系，我们选取了2011年数据表中深圳、福建、北京等产地，统计出其不合格率，建立了产地距离与产品不合格率之间的线性回归模型，通过Matlab软件求出产地距离与产品不合格率之间二次回归方程，且F概率p满足小于置信水平的成立条件，由此得出食品产地与质量产地成二次的波峰关系。对于食品质量与抽检地点的关系，我们选取了了深圳市内的8个抽检地点，并分别统计了2011年、2012年8个抽检地点的不合格率，运用单因素方差分析法，利用Matlab软件，以Anoval函数求解。求出p-value，显著性水平取0.05作为标准来判断是否有显著性，最终求得，得出抽检地点与食品质量无显著性关系的结论。对于问题3，从以上两个问题可以看出，,食品的抽检涉及到很多方面，如果在抽检过程中方向性目的性不强，会导致在投入大量资金、人力、物力后往往得不到准确的抽检结果，所以在模型建立求解中，我们选择一个具体食品领域（餐饮），通过对各个因素进行量化分析，采用方差分析及matlab逐步线性回归求解的方法分析，最后得出实际在2011年到2013年间深圳市餐饮食品20多次大大小小的抽检中，其抽检样本量从几十到上千，但其对最后的食品不合格率变化的影响较小。在逐步线性回归分析中又得到抽样样本量、食品添加剂、重金属、菌落四种因素和抽检不合格率的回归方程，通过对回归方程系数的分析可得，食品添加剂对食品不合格率变化的影响近乎是重金属和菌落影响之和，由此结论可对抽检方案进行合理调整。关键字：层次分析法 线性回归 方差分析，逐步线性回归分析二、模型假设(1)假设食品产地、抽检地点对食品质量存在着某种关联；(2)假设抽样过程是随机的，数据是准确的，地理或人为因素对抽样食品不合格性的影响忽略不计。（3）假设在数据采取的几年间，深圳市的食品抽检合格标准未变化。三、符号说明2010年、2011年、2012年因微生物食品不合格的不合格率的平均值2010年、2011年、2012年因添加剂食品不合格的不合格率的平均值2010年、2011年、2012年因重金属食品不合格的不合格率的平均值2010年、2011年、2012年的微生物在生产领域不合格率的平均值2010年、2011年、2012年的微生物在流通领域不合格率的平均值2010年、2011年、2012年的微生物在餐饮领域不合格率的平均值2010年、2011年、2012年的添加剂在生产领域不合格率的平均值2010年、2011年、2012年的添加剂在流通领域不合格率的平均值2010年、2011年、2012年的添加剂在餐饮领域不合格率的平均值2010年、2011年、2012年的重金属在生产领域不合格率的平均值2010年、2011年、2012年的重金属在流通领域不合格率的平均值2010年、2011年、2012年的重金属在餐饮领域不合格率的平均值生产领域中年不合格率占3年不合格率和的比重流通领域中年不合格率占3年不合格率和的比重餐饮领域中年不合格率占3年不合格率和的比重准则层1对目标层的成对比较矩阵准则层2对准则层1的成对比较矩阵方案层2对准则层2的成对比较矩阵 成对比较矩阵最大特征值不一致比率第层对第层的权向量第层对第层的权向量构成的矩阵相关系数置信概率置信区间FF统计量值偏回归系数剩余方差因素各类数据源的平方和各类数据相应的自由度各类的均方值大于F的概率各组均值对总方差的偏差平方和各组数据对均值偏差平方和的总和四、问题分析 对于问题1，为了针对深圳市近3年来食品安全变化情况的的整体趋势，同时考虑微生物、添加剂、重金属等在生产领域、流通领域、餐饮领域的造成的影响。我们根据统计好数据，将各个领域2010年、2011年、2012年合格抽检数据，以及微生物、添加剂、重金属等不合格抽检数据进行整理和分析，考虑微生物、添加剂、重金属等对各个领域的影响，以及各个领域对这3年的影响，我们小组利用层次分析法，利用相互关联的不合格率之间的比重进行重要性的强弱对比，构造成对比较矩阵，然后利用Matlab编程求得各个成对比较矩阵的特征值，进行一致性检验，符合一致性检验后，利用各个特征向量求得方案层对目标层的组合权向量，然后作出综合评价。对于问题2中食品产地与食品质量的关系，我们统计了2011年各主要产地的不合格率，并对其进行线性回归进行拟合来判断二者是否有关系，而对于抽检地点已经不能用距离与产地不合格率来进行线性拟合，所以我们统计8个抽检地区的2011、2012的数据量，利用单因素分差分析方法来判断抽检地点与食品质量是否具有显著性关系。对于问题3，要找到一种合适的抽检方法或是对现有的抽检方法进行合理的调整，我们需要对一个抽检较长时间、抽检量较大，方向针对性较强的数据体进行分析。最后得到较为普遍的结论，在对抽检的方法就行调整；餐饮是食品安全的重要环节，对餐饮环节的抽检至关重要，所以我们选取2011年到2013年深圳市的餐饮抽检结果数据进行分析，首先我们针对抽检量对于抽检不合格率的影响进行方差分析，得到其影响程度；其次，综合所有因素进行逐步线性回归分析，筛选出对食品不合格率影响相对较大的几个因素建立线性回归模型；五、模型建立与求解5.1 食品安全变化趋势层次分析模型的建立与求解由于我国居民消费的食物种类非常复杂，为便于数据整理和分析，我们将3年的食品抽查数据按照生产领域、流通领域、餐饮领域进行整理得到如下表格：抽检合格抽检不合格微生物添加剂重金属生产领域2010年264869392832011年3349104402832012年1283211441流通领域2010年296679294142011年5143436154269282012年51422454313411餐饮领域2010年28341035732122011年59534451811271352012年8937152646227着重考虑3大领域中微生物、添加剂、重金属对食品安全问题的影响，为得到对3年的食品安全问题的综合评价，我们综合考虑微生物、添加剂、重金属分别在生产领域、流通领域、餐饮领域的影响，然后通过生产领域、流通领域、餐饮领域对2010年、2011年、2012年的影响，建立相关的层次分析模型。1）食品质量安全综合评价(选择安全问题最严重的年份)的阶梯层次结构A(综合评价)方案层目标层B1（微生物）B2（添加剂）B3（重金属）准则层1C3（餐饮领域）C2（流通领域）C1（生产领域）准则层2D1（2010年）D3（2012年）D2（2011年）2）建立B-A层对比较矩阵建立B-A成对比较矩阵的依据是：要获得B层对A层的重要性的强弱，我们通过2010年、2011年、2012年的因微生物、添加剂、重金属造成食品不合格的不合格率的平均值间的比重来构建成对比较矩阵：0.0154470.0172550.005138得到成对比较矩阵：3）建立C-B成对比较矩阵建立C-B成对比较矩阵的依据是：要获得C层对B层的重要性的强弱，我们通过2010年、2011年、2012年的微生物、添加剂、重金属分别在生产领域、流通领域、餐饮领域的不合格率的平均值间的比重构建成对比较矩阵：0.0122250.0071610.027398得到成对比较矩阵：0.0150370.0288520.002801得到成对比较矩阵：0.0182460.0125220.009386得到成对比较矩阵：4）建立D-C成对比较矩阵建立D-C成对比较矩阵的依据是：要获得D层对C层的重要性的强弱，我们通过2010年、2011年、2012年抽检的食品分别在生产领域、流通领域、餐饮领域的不合格率相互比较构建成对比较矩阵：0.3545950.4205430.224648得到成对比较矩阵：0.1780730.5224840.304062得到成对比较矩阵：0.2860020.5672190.136384得到成对比较矩阵： 5）一致性检验 成对比较矩阵通常是不一致阵，但是为了能用它的对应特征根的特征向量作为被比较因素的权向量，其一致程度应在一定允许范围内，因此可做一致性检验来进行讨论。一致性指标：一致性指标： 随机一致性指标： 随机一致性指标的数值1234567891011000.541.321.411.451.491.51 对于的成对比较矩阵A，将它的一致性指标与同阶（指n相同）的随机一致性指标之比称为一致性比率，当时认为A的不一致程度在允许范围内，可用其特征向量作为权向量。6)利用Matlab编程求解得到成对比较矩阵的特征向量，即权向量 第1层只有一个因素，第2,3层分别有个因素，第2,3,层对第1,2层的权向量分别为以为列向量构成矩阵则第3层对第1层的组合权向量为 一般地，若共有 层，则第层对第1层（设只有一个因素）的组合权向量满足：其中是以第层对第层的权向量为列向量组成的矩阵。于是最下层（第层）对最上层的组合权向量为：二、模型求解利用Matlab编程求解成对比较矩阵的特征向量B的权向量其最大特征值，不一致性比率：，因此可用其特征向量作为权向量。综合评价问题第3层的计算结果1230.27100.31870.44340.15300.61530.38740.57600.06600.16923.00373.00123.07350.00320.00100.0633综合评价问题第4层的计算结果1230.35900.16560.28570.41380.51460.57140.22710.31980.14293.00013.001230.000090.00100则方案层对目标层的组合权向量为：即：5.2.1食品产地与食品质量的关系模型建立和求解 我们统计了深圳、东莞、中山、广州、珠海、佛山、汕头、广西、湖南、福建、上海、重庆、四川、北京14个主要产地2011年产品的合格数和不合格数，并上网统计各产地与深圳市的距离，得到各产地2011年产品的不合格率，如表5.1表5.1产地合格不合格不合格率距离深圳41061710.0399810东莞615280.04354657中山150110.068323121广州398240.056872147珠海5950.078125161佛山163160.089385169汕头174180.09375494广西127150.105634736湖南202250.110132790福建347390.101036846上海205250.1086961680重庆101100.090091844四川154170.0994152048北京71100.1234572372对于二者是否存在一定的相关性，我们采取了线性回归模型，由于距离与不合格率的量级相差较大，所以将距离统一缩小1000倍，通过Matlab软件得到产地与食品不合格率的散点图：图5.1 散点图由散点图可以看出随着距离增大，食品不合格率增大，到了一定幅度有下降，故该曲线用的是二次模型即，其中为偏回归系数，通过Matlab程序中regress函数拟合线性回归函数，求得 置信区间为 置信区间为 置信区间为相关系数 即拟合的二次回归方程为，取置信水平，由于，其中相关系数指因变量y（不合格率的）可由模型确定， 查表得F检验临界值，F值超过F检验临界值，故模型可用。运用函数rcoplot进行残差分析：图5.2 残差分析图从残差图可以看出除第14个数据外其余数据的残差离零点均较近，且残差的置信区间均包含零点，这说明回归模型能较好的符合原始数据而第14个数据可视为异常点(而剔除)。对得到的回归模型进行预测及作图，得到预测比较图如图所示：图5.3 预测比较图 结果分析：有线性拟合的结果可以看出，食品产地与食品质量存在先上升后下降的关系，对拟合曲线进行用diff函数求导，并用solve函数求解该二次函数的极值点，得该拟合函数的极值点对应的横坐标，即距离为，有拟合结果可以看出，当距离小于1600米时食品质量的不合格率随着产地距离的增大而增大，但超过1600以后，食品质量的不合格率随着产地距离的增大而减小，由于该模型只是一个粗略的模型，距离的精确度不一定精确，但是食品产地与食品质量存在先上升后下降的关系成立5.2.2食品抽检地点与食品质量的关系模型建立及其求解 我们选取了深圳市内抽检的8个地点，罗湖区、福田区、南山区、盐田区、宝安区、龙岗区、光明新区、坪山新区，要判断食品抽检地点与食品质量是否具有显著性关系，我们采用单因素方差分析法去解决。单因素方差分析法: 只考虑一个因素A 对所关心的指标的影响，A 取几个水平，在每个水平上作若干个试验，试验过程中除A 外其它影响指标的因素都保持不变（只有随机因素存在),我们的任务是从试验结果推断，因素A 对指标有无显著影响，即当A 取不同水平时指标有无显著差别。A 取某个水平下的指标视为随机变量,判断A 取不同水平时指标有无显著差别，相当于检验若干总体的均值是否相等。设 A取n 个水平,在水平下总体 服从正态分步N(,),i=1,.,n,这里,未知，可以互不相同，但假定有相同的方差，又设在每个水平下作了次独立试验，即从中抽取容量为的样本，记作服从N(,)，i=1,n,j=1, 且且相互独立。将这些数据列成表1表5.2 单因素试验数据表因素水平12jn.根据上述理论，首先我们对数据进行处理,我们统计了2011、2012年分别在8个抽检地区的合格与不合格数，采用单因素方差分析法，我们将食品的不合格率唯一考虑的因素A，运用matlab软件编程每个抽检不合格率的平均分，见下表：表5.4 2011、2012年个抽检地区不合格率表 时间地区20112012不合格率不合格率罗湖区0.061056110.028846福田区0.075174830.023881南山区0.070996980.016308盐田区0.060606060.018939宝安区0.081065090.029535龙岗区0.108877720.027037光明新区0.084398980.029891坪山新区0.066666670.038741运用matlab软件对数据处理编程得出以下结果，标准ANOVA表分析见下表：图5.4 抽检地盒型（box）图图5.3 抽检地ANOVA表可构造方差表来完成计算： 表5.3 方差表方差来源平方和自由度均方比值显著性因素A的影响s-1随机因素的影响n-s总 和n-1总平均：；总离差平方和：；组内平方和（误差平方和）：，随机因素的影响；组间平方和（因素平方和）：，水平差异的影响；H0的拒绝域为：检验结果：高度显著：；显著：；有一定影响：；无显著影响：。通常情况下，实验结果p达到0.05水平或0.01水平，才可以说数据之间具备了差异显著或是极显著。在作结论时，应确实描述方向性（例如显著大于或显著小于）。sig值通常用 P0.05 表示差异性不显著。在此我们以0.05作为显著性水平标准，当,时不拒绝，抽检地区中的ANOVA表中ProbF栏p值为0.9944 0.05,且F值很小，故接受，故无有显著影响。故抽检地点与食品质量并不显著性关系。5.3.1问题三模型的建立 一、 方差分析（讨论不同的抽样样本量对于抽样结果即食品不合格率的影响，这里采用方差分析法。）（1）总平均不合格率X=x1+x2+x3+x4+x55;（2）总变差平方和S总=ST=i=1Pj=1r(Xij-X);（3）组间变差平方和SA=i=1Pj=1r(Xi-X)2;（4）组内变差平方和：S=i=1Pj=1rXij-X2=SA1+SA2+SA3+SA4+SA5（5）自由度：fA=P-1，f= P.r-1，fT=fA+f（6）平均平方（均方）S=SP（r-1）；SA=SAP-1；（7）显著性检验（F检验）：F=SAS；二、逐步线性规划考虑所有变量对于食品抽检不合格率的影响进行逐步线性回归分析，同样抽样近年来的餐饮环节的抽检数据作为分析对象，考虑包括抽样样本(x1)、添加剂(x2)、重金属(x3)、菌落(x4)四种因素对食品不合格率(Y)的影响，进行逐步线性回归分析影响；建立线性模型：Y=a1.x1+a2.x2+a3.x3+a4.x4+a5;(a1a5是方程系数)；5.3.2问题三模型的求解和结果分析一、方程分析求解； 20112013年餐饮抽检样本量与食品不合格率关系表A抽样样本大小A1=250250 =A2 =400400=A3 =500600=A4 =9001000=A5 =2000不合格率0.0681820.0090634440.125786160.0173661360.04430929600.0029411760.125786160.0127551020.0292397660.0683760.0607734810.280172410.0209059230.0391547550.1322030.0208877280.005154640.0115667720.043577982平均不合格率0.067190.0234164570.134224850.0156484830.03907045由表可得：（1）总平均不合格率：X=0.06719+0.02341+0.13422+0.015648+0.039075=0.055910111；（2）总变差平方和：S总=ST=i=1Pj=1r(Xij-X)=(0.06818-0.05591)2+(0-0.05591)2+(0.06837-0.05591)2(0.04357-0.05591)2=0.08595;（3）组间变差平方和：SA=i=1Pj=1r(Xi-X)2=ri=1P(Xi-X)2=4(0.06719-0.05591)2+(0.02341-0.05591)2+(0.13422-0.05591)2+(0.01564-0.05591)2+(0.03907-0.05591)2=0.03688；（4）组内变差平方和：S=i=1Pj=1rXij-X2=SA1+SA2+SA3+SA4+SA5；SA1=0.068182-0.067192+0-0.067192+0.068376-0.067192+0.132203-0.067192=0.00874SA2=0.00906-0.023412+0.00294-0.023412+0.06077-0.023412+0.02088-0.023412=0.00202；SA3=0.03810; SA4=0.00005; SA5=0.00014;所以：S=SA1+SA2+SA3+SA4+SA5=0.04907； （5）自由度：SA的自由度是fA=P-1=4；S的自由度是f= P.r-1=15；ST自由度是fT=fA+f=P.r-1=19；（6）平均平方(均方)SA=SAP-1=0.036884=0.00922；S=SP（r-1）=0.0490715=0.00327；（7）显著性检验（F检验）F=SAS=0.009220.00327=2.81851；通过查表得到：=0.05时Fa=3.29，=0.01时Fa=5.42；因为： F=2.8185513.29;所以，可以得出结论：抽样样本大小在2000以下变化时对食品的不合格率有一定影响，但不是很显著，故在许多抽样合格率较为稳定的食品环节适当减少抽检的数量对抽检结果影响不大，同时也可以减少抽检费用。二、逐步线性规划的求解过程： 2011年至2013年餐饮食品抽检表序号抽样样本量添加剂重金属菌落不合格率12085384570.04357822184930.06837636796050.017366434014230.06077351546108440.0391556256231180.13220373282000.00906383758000.020888913808060.00933210823030.0681821184361200.020906125792010.0051551494010100.011567157745050.012755162179181590.01891171328100300.02924183391000.002941191100133280.04430920417158390.12578621417158390.125786逐步线性回归编程所得结果：由上表可知：x1、x2、x3、x4即抽样样本（x1）、添加剂（x2）、重金属(x3)和菌落（x4）四种因素对食品不合格率（Y）变化影响显著，根据表中的数据得到其线性模型为：且a1=-6.0232210-5；a2=0.00217956；a3=0.00143211；a4=-0.00122579；a5=0.0458531；故建立线性方程为：Y=-6.0232210-5.x1+0.00217956.x2+0.00143211.x3-0.00122579.x4+0.0458531由此方程可以看出添加剂、重金属、菌落这三个影响因素中，菌落和重金属对食品不合格率变化的影响相对较低，而添加剂的影响则较为显著，其相关系数近似是前面两者系数的和，所以，在抽检时可以将一半的抽检费用花费在添加剂的检验上，重金属和菌落的抽检可各占14左右。六、模型的评价6.1 模型的优点（1） 模型建立的原始数据均来自深圳市市场监督管理局网站，准确率高，有很好的权威性，对各个领域进行宏观数据统计，科学合理，有说服力。（2）模型1的建立，使得能够综合考虑各个因素对食品安全评价的影响。通过对数据的科学处理和分析，将对综合评价的贡献的重要性进行合理的量化；此次模型成对比较矩阵的获得，是通过数据的比值得到，最终得到2011-2013年食品不合格严重程度的比重，得到食品安全变化的整体趋势，考虑问题全面和合理。6.1 模型的缺点（1） 模型在数据统计时未考虑除添加剂、微生物、重金属以外物质对食品安全问题的影响，对模型有一定的影响，但是不明显。七、参考文献1姜启源，谢金鑫，叶俊.数学建模，北京：高等教育出版社，20112杨虎，刘琼荪，钟波，概率论与数理统计，北京：重庆大学出版社，2011.53matlab多元与非线性回归即拟合问题regress、nlinfit /view/0fb58b18fc4ffe473368abe0.html 2013.5.234 matlab方差分析教案/view/c1c630ea81c758f5f61f6716.html 2013.5.245 数据下载网站：（深圳市市场监督管理局网站）附录1.食品产地与食品质量的关系模型的matlab程序：x=0 57 121 147 161 169 494 736 790 846 1680 1844 2048 2372/1000;y=0.039981295 0.043545879 0.068322981 0.056872038 0.078125 0.089385475 0.09375 0.105633803 0.110132159 0.101036269 0.108695652 0.09009009 0.099415205 0.12345679;plot(x,y, k+) %做散点图X=ones(length(y),1),x,(x.2); b,bint,r,rint,stats=regress(y,X) rcoplot(r,rint) %进行残差分析z=b(1)+b(2)*x+b(3)*(x.2 );plot(x,y, k+,x,z, r) %比较作图2.求极值点的matlab程序：syms x;y=0.0566+0.0722.*x-0.0225*(x.2);Y1=diff(y) %求一阶导Xmin=solve(Y1) % 求极值点 A=0.061056106 0.028846154;0.075174825 0.023880597 ;0.070996979 0.016307513;0