食品质量安全抽检数据分析论文.doc

数学建模校内竞赛论文 论文题目：食品质量安全抽检数据分析组号：284#成员：王林炜 陈小红 葛敏雪 选题：A题 姓名 学院 年级 专业 学号 联系电话 数学分析 高等代数 高等数学 线性代数 概率统计 数学实验 数学模型 CET4 CET6 王林炜 自动化 2011 物流工程 20115172 15922768472 / / 86 86 / 85 / 567 / 陈小红 自动化 2011 物流工程 20115223 13637985570 / / 95 86 / 92 / / / 葛敏雪 自动化 2011 物流工程 20115174 13512354364 / / 91 90 82 / / 519 / 2013-05-21 食品质量安全抽检数据分析摘要“民以食为天”，食品安全关系千家万户，本文根据深圳市给出2010年、2011年、2012年的食品质量安全抽检数据，针对深圳市的食品质量安全做出了分析。通过对数据的研究，我们对数据进行了处理，将主要食品分为六大类，分别是：蔬果类、肉蛋类、水产类（、淀粉类、调味类、豆乳类。在建立数学模型时，应用了数据筛选和均值统计法、层次分析法、回归分析法、数据拟合法、概率统计分析。针对问题一，我们将影响食品安全的因素分为微生物、重金属、添加剂、其他这四大类。先用数据筛选和均值法对数据进行处理，从而得到2010年、2011年和2012年每个季度食品因微生物、重金属、添加剂、其他而导致的不合格率，并用图表表示出每一类三年的变化趋势。然后用层次分析法，定量地求出三年食品安全系数。三年内食品安全质量有波动，总体上说，安全情况有下降趋势，食品安全有待提高。年份201020112012食品安全系数0.19320.56760.2238等级针对问题二，由于微生物、重金属、添加剂、其他这四方面的不合格导致了食品质量的不过关，但是这四类直接因素又由其它因素影响，通过对数据的分析，我们将影响食品安全的这些因素归结为检测地点、食品产地、抽查季节、生产环节、流通环节、餐饮环节六个要素。通过建立回归模型，运用MATLAB软件，求出六大要素与食品质量间的线性关系，定量求出各因素在影响食品安全中所占影响度大小和影响规律。最终得到，在食品质量影响因素中食品产地影响最大，食品加工环节影响次之，季节影响最小，抽查地点几乎无影响。针对问题三，根据第一问和第二问，我们可以得到不同因素的影响程度是不同的。为了科学有效地反映食品安全情况，合理简化抽检方法，我们减少稳定数据的测量量，增加不稳定数据的测量量，使测量环节有效减少，重要环节的监测力度加大。通过使用SPSS软件，求解出新选出的各项因素与食品安全的线性关系，并通过回归曲线进行检验，论证了新的检验方法求得出的数据与原检验方法下的数据误差趋于0，新方法具有科学性。最终得到，新的抽检方法为减少调料类、豆乳类、蔬果类、淀粉类食品抽检次数，增加水产类、肉蛋类抽检次数。关键词：统计 层次分析 回归分析 拟合 残差分析 蒙特卡洛法 一、 问题重述 “民以食为天”，食品安全关系到千家万户的生活与健康。随着人们对生活质量的追求和安全意思的提高，食品安全已成为社会关注的热点，也是政府民生工程的一个主题。城市食品的来源越来越广泛，人们消费加工好的食品的比例也越来越高，因此除食材的生产收获外，食品的运输、加工、包装、贮存、销售以及餐饮等每一个环节皆可能影响食品的质量与安全。另一方面，食品质量与安全又是一个专业性很强的问题，其标准的制定和抽样检测及评价都需要科学有效的方法。深圳是食品抽检、监督最统一、最规范、最公开的城市之一。请下载2010年、2011年和2012年深圳市的食品抽检数据（注意蔬菜、鱼类、鸡鸭等抽检数据的获取），并根据这些资料来讨论：1. 如何评价深圳市这三年各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势；2. 从这些数据中能否找出某些规律性的东西：如食品产地与食品质量的关系；食品销售地点（即抽检地点）与食品质量的关系；季节因素等等；3. 能否改进食品抽检的办法，使之更科学更有效地反映食品质量状况且不过分增加监管成本（食品抽检是需要费用的），例如对于抽检结果稳定且抽检频次过高的食品领域该作怎样的调整？二、问题分析对于问题一，由于调查中食物种类复杂，为了方便问题的讨论，我们将食品分为六大类，即淀粉类，水产类，蔬果类，肉蛋类豆乳类和调味品，对数据进行筛选。在考虑食品安全情况的影响因素时，根据数据分析，我们主要分析四大主要因素，分别为微生物、重金属、添加剂和其他。从附录中得到深圳市2010，,2011,2012三年抽检样本数据，并用EXCEL制表进行统计和整理，求出每个因素导致的每个季度的食品不合格率。根据相应的数据制作出每个因素安全情况三年的变化趋势。再用层次分析法，从定量的角度求出三年的总体变化趋势。对于问题二，第一问中涉及的影响因素是食品是否安全的直接因素，但这些因素又是由其它相关因素决定。分析数据，我们将影响食品安全的这些因素归结为检测地点、食品产地、抽查季节、生产环节、流通环节、餐饮环节六个要素。这一问要求的是这些要素与食品质量的规律性关系，即线性关系。通过建立回归模型，来定量求出各因素在影响食品安全中所占影响度大小和影响规律。对于各因素对食品安全的影响程度，利用偏回归平方和衡量每个变量在回归中所起的作用大小（即影响程度），对每个变量的偏回归平方和进行了计算，最后把影响程度由深到浅的各因素进行了排序。最后进行回归方程的显著性检验，将总的平方和分解为回归平方和和剩余平方和，运用F检验法来判定线性回归方程的显著性。对于问题三，要改进原有抽检方法，找到一个新方法，使其在成本相同的情况下，测得数据结果更准确，在测得数据结果和老方法接近下，成本更低。根据第一问和第二问得到的数据，得到不同因素的影响程度是不同的。减少稳定数据的测量量，加大不稳定数据的测量量，减少测量环节，加重必要环节。由于我们获取的数据均是在原抽检方法下得到的数据，而没有真正标准真实的数据，所以无论我们用什么样的新方法，只能在原抽检下的数据（或有关预测、推测数据）进行验证。故新方法在测量不稳定数据的结果与原方法不具可比性。但稳定数据的结果两种方法下的测量数据基本一致。新方法得到后，检验新方法的科学性、合理性。拟合曲线确定其未来检测值，采用蒙特卡洛法抽检，带入回归方程进行验证新方法是否合理。三、问题假设1、 假设模型求解过程中所用的数据都是合理的2、假设不考虑抽检时抽检人员的人为影响3、假设每次抽检时遵循公平性原则4、假设所有食品取主要食品即六大类，其余食品对质量安全分析无影响5、假设抽检时所有环节、所有因素的成本和工时相同四、 参数及符号说明符号表示的意义备注Ai 目标层指标i=1Bi 准则层指标i=1,2,3,4Ci方案层指标i=1,2,3 第 i层权重i=1,2,3标准化后的第 i层权重i=1,2,3入一致性指标过度数C.I一致性指标R.I.平均随机一次性指标C.R.一致性比例指标bij 每个因子的评分i=1,2,3j=1,2,3,4Ci综合评价指数i=1,2,3PI食品安全质量综合评价指数影响食品质量安全的因素。i=1,2,.6各食品安全系数a=1,2.6相互独立且服从同一正态分布的随机变量变量的偏回归平方和i=1,2,.6回归平方和剩余平方和待估计系数参数的回归值参数的最小二乘估计五、 模型的分析、建立与求解（1）问题一的分析与模型的建立1、统计、计算、绘制图表来表示出六大类食品，在2010、2011、2012三年中各季度下抽检样品中不合格食品里，微生物、重金属、添加剂的变化趋势；通过数据筛选和统计得出附表一，根据附表一，计算求解出各个季度下 微生物、重金属、添加剂、其它等因素不合格的概率年份201020112012季度二三四一二三四一二三四微生物0.0330.0510.1410.1050.020.1870.3420.0670.0590.0360.074重金属0.0390.0300.3890.1490.1780.16200.1220.0770.134添加剂0.0380.0300.1170.0380.1780.1190.0980.0860.0370.092其它0.050.0480.1410.0410.01500.0270.1020.0670.0230.028后根据表格，绘制图表来表示三年中各季度下抽检样品中不合格食品里，微生物、重金属、添加剂的变化趋势：各因素不合格的概率年平均201020112012微生物0.0750.1640.059重金属0.0230.2200.083添加剂0.0230.1130.078其它0.0800.0210.055从图表得出结论：三年来，因微生物 、重金属、添加剂导致的不合格率变化趋势有波动，食品安全整体有下降，从总体上说，食品安全情况没有改善。2、 根据三年中六类食品的食品安全情况，采用层次分析法，定量结合定性的评价出由于微生物、重金属、添加剂超标而导致食品不安全的年度食品安全情况。对深圳市这三年各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势做出定量的综合评价，首先对3年的数据处理按季节分为12个子样本点的抽检值进行分析，得到各子样本食品质量安全情况。运用层次分析法决定食品安全单位指标在综合评价中权重。 建立层次结构模型（1）构造各季节食品安全质量综合评价的阶梯层次结构根据层次分析法的基本步骤,建立食品质量安全评价递阶层次模型,如下图所示。食品安全综合评价A1 AA 微生物B1 重金属B2 添加剂B3 其它B42010年食品安全情况C12012年食品安全情况C32011年食品安全情况C2目标层A准则层B方案层C 图：食品综合评价递阶层次结构（2）建立两两比较的A-B层判断矩阵并进行一致性检验。建立A-B层判断矩阵的依据是：考虑B层的四个因素对食品质量安全综合评价的重要性时，根据在查找数据知年均指标微生物最重要，且B1比重金属、添加剂、其它含量等值重要。而在重金属、添加剂、其它含量等值中，而各类食品都用到添加剂，所以，添加剂比重金属、其它重要。而其它比重金属重要。利用19标度法(见下表)进行成对比较,同时参考专家意见,确定各因素之间的相对重要性并赋以相应的分值,构造出各层次中的所有判断矩阵,并计算权向量和一致性检验。表2 19 标度的意义标度aij 意义1Ci 与Cj 的影响相同3 Ci 比Cj 的影响稍强5Ci 比Cj 的影响强7Ci 比Cj 的影响明显地强9Ci 比Cj 的影响绝对地强2 ,4 ,6 ,8为上述两判断级的中间值1 ,1/ 2 , ,1/ 9 Ci 较Cj 的影响之比与上述说明相反 1、 根据判断矩阵标度内容及其含义，可得A-B层判断矩阵B： B1 B2 B3 B4 1 4 2 3 B1 1/4 1 1/2 3/4 B2 B= 1/2 2 1 3/2 B3 1/3 4/3 2/3 1 B4 其中B矩阵中的每一个元素Bij根据调查相关数据定性判断出：2、 用方根法计算出判断矩阵B的权重值并进行归一化： 方根法（几何平均法）： Wi= i=1,2,3,.,n归一化： =计算步骤：1、B中的元素按行相乘得一新向量；2、将新向量的每一个分向量开n次方；3、将所有向量归一化即为权重向量。3、计算判断矩阵B的最大特征值并进行一致性检验. 计算一致性指标C.I 其中，n为判断矩阵的阶数计算一致性比例其中为平均一致性指标，当时判断矩阵具有好的一致性。找查相应的平均随机一致性指标R.I。平均随机一致性指标n12345678910R.I000.520.891.121.261.361.411.461.49 经计算得图表：AB1B2B3B4=4C.I.=0R.I.=0.89C.R.<0.1B114232.210.484.00B21/411/23/40.550.124.02B31/2213/21.110.243.99B41/34/32/310.740.163.99(3)建立两两比较的B-C层判断矩阵并进行一致性检验。根据途径一的表格中数据 （表1-1）可完成表格：方法和计算B判断矩阵一致。B1C1C2C3=3C.I.=0R.I.=0.52C.R.<0.1 C110.4571.2710.8340.253.00C22.18712.7801.8250.553.00C30.7870.36010.6570.203.00B2C1C2C3=3C.I.=0R.I.=0.52C.R.<0.1C110.1050.2770.3080.073.00C29.56512.6512.9380.673.00C33.6090.37711.1080.253.00B3C1C2C3=3C.I.=0R.I.=0.52C.R.<0.1C110.2040.2950.3920.113.00C24.91311.4491.9240.533.00C33.3910.69011.3280.363.00 B4C1C2C3=3C.I.=0R.I.=0.52C.R.<0.1 C110.3981.4550.8340.243.00C22.51313.6552.0940.603.00C30.6880.27410.5730.163.00(4)三年中食品安全评价结果表 1、列表计算 Bi biB1B2B3B40.480.120.240.16C10.250.070.110.240.1932C20.550.670.530.600.5676C30.200.250.360.160.22382、 2、深圳市食品质量系数（PI）分级 （表2-2）(5）市场食品质量安全的综合评价指数的计算 PI=根据深圳市食品安检局标准处推荐的评价标准（GB38382002食品安检质量标准以及食品质量分级）为参考依据，见下表1.3表1.3深圳市食品质量系数（PI）分级质量安全值级别食品质量评价<0.1良好0.10.25较好0.250.4一般0.40.55轻劣食品0.550.80重劣食品>0.80劣严重劣食品这样就得到了2013年第一至四期食品安全抽样检验产品抽查样本的食品安全的优劣状况，并做出定量的综合评价。其它年份的食品安全情况同理可求。求解程序见附录1。表1.4为求得的各个抽查样本的食品质量安全状况和定量综合评价。其中C.R.为一致性指标(当C.R.<0.1时有好的一致性)详细的计算结果见附录2。经计算最终得到了三年的食品安全评价等级表：年份201020112012食品安全系数0.19320.56760.2238等级（2）问题二的分析与模型的建立1、 处理数据 根据题设问题，我们将数据表一进行分析、筛选、整理，得题设二所需数据。我们将表中的数据按年度和季分为12子样本，对每个子样本根据回归模型需要，计算出的样本均值，以此类推，最终得到88个子样本.（回归模型需要：因为各因素表现基本成0-1类型，除季节（1-4）与抽检地(1-8)外，所以我们可以先按季节每年分为4个子样本空间，在每个子样本空间中在按抽检地分类，得到次子样本，后将每个次子样本进行均值处理，即可得回归模型所需数据。但计算之前需要利用残差向量进行数据筛选）。我们先通过SPSS软件（程序见附录1）对原始数据进行检验，对残差向量进行分析，得到了残差向量分析图，剔除其中的异常点。假设已得到筛选后的88组数据。2、设随机变量假如变量与另外6个变量，, 的内在联系是线性的，它的第次试验数据是 =1,2，,6 （1）那么这一组数据可以假设有如下的结构式： （2）其中，是7个待估计参数，, 是8个可以精确测量的一般变量，是88个相互独立且服从同一正态分布的随机变量，这就是多元线性回归的数学模型。 令， ， ， 那么多元线性回归的数学模型（2）可以写成矩阵形式 （3）其中是72维随机向量，它的分量是相互独立的。3、参数的最小二乘估计为了估计参数，我们采用最小二乘估计法。设 分别是参数，的最小二乘估计，则回归方程为 （4）由最小二乘法知道， 应使得全部观察值与回归值的偏差平方和达到最小，即使 (5)所以是 的非负二次式，最小值一定存在。根据微积分学中的极值原理。 应是下列正规方程组的解： (6)显然，正规方程组的系数矩阵是对称矩阵，用来表示，则，且其右端常数项矩阵亦可采用矩阵和来表示：。所以可以得到回归方程的回归系数： (7)4、由于利用偏回归平方和可以衡量每个变量在回归中所起的作用大小（即影响程度），设是p个变量所引起的回归平方和，是p-1个变量所引起的回归平方和（即除去），则偏回归平方和为：=-=-= （8）就是去掉变量后，回归平方和所减少的量。1、数据筛选通过MATLAB（程序见附录1）作图如下： 经过残差分析后我们从原先96组数据中筛选出可靠88组数据，数据表格见附录1-1 2、回归方程的求解用筛选后的数据，根据回归模型编程求解（程序如附表） 解得：回归系数为：= -1,.73 = 0.016 = O.001 =0.973 =1.003 =1.004 =0.044回归方程为：3、偏回归平方和的比较运用SPSS进行编程（程序见附录2），得到各因素的偏回方和：0.7670.9730.000.000.000.13结论：在食品质量影响因素中食品产地影响最大，食品加工次之，季节影响最小，抽查地点几乎无影响。（3）问题三的分析与模型的建立1、用数据拟合曲线来验证其 抽减次数减少后的结果，利用残差和来看。2、假如变量与另外5个变量，, 的内在联系是线性的，它的第次试验数据是 =1,2，,5 （1）那么这一组数据可以假设有如下的结构式： （2）其中，是7个待估计参数，, 是8个可以精确测量的一般变量，是96个相互独立且服从同一正态分布的随机变量，这就是多元线性回归的数学模型。 令， ， ， 那么多元线性回归的数学模型（2）可以写成矩阵形式 （3）其中是72维随机向量，它的分量是相互独立的。3、参数的最小二乘估计为了估计参数，我们采用最小二乘估计法。设 分别是参数，的最小二乘估计，则回归方程为 （4）由最小二乘法知道， 应使得全部观察值与回归值的偏差平方和达到最小，即使 (5)所以是 的非负二次式，最小值一定存在。根据微积分学中的极值原理。 应是下列正规方程组的解： (6)显然，正规方程组的系数矩阵是对称矩阵，用来表示，则，且其右端常数项矩阵亦可采用矩阵和来表示：。所以可以得到回归方程的回归系数： (7)4、由于利用偏回归平方和可以衡量每个变量在回归中所起的作用大小（即影响程度），设是p个变量所引起的回归平方和，是p-1个变量所引起的回归平方和（即除去），则偏回归平方和为：=-=-= （8）就是去掉变量后，回归平方和所减少的量。2、由于表格数据限制，新方法分析看见附表中的数据分析。六、对于模型的评价6.1模型的优点 本文对于各种因素对于食品质量安全的影响建立了多元线性回归模型，全面综合考虑了各个方面的因素，避免了单一因素分析的不准确性，得出了合理的数学模型。并且通过各因素的显著性分析，找到了影响食品质量安全的主要因素，较符合实际情况，模型可靠，并且模型相对简单，利于操作；该方法不仅适用于本题，也适用于其他方面的数据预测，有实际背景，可运用于实践，具有广泛适用性。6.2模型的缺点本文忽略了除了所给因素之外的因素对食品安全的影响，与实际问题存在偏差。同时是在假设各因素相互独立的情况下对食品质量的影响进行分析，可能会导致误差七、参考文献1 马新民，王逸迅. 概率与数理统计M. 北京：机械工业出版社，2010.2 刘卫国. MATLAB程序设计与应用M. 北京：高等教育出版社，2009.3 茆诗松. 回归分析及其试验设计M. 上海：华东师范大学出版社，1986.4 姜启源，数学模型M, 北京，高等教育出纳版社，2006附录一、经过残差分析，我们把96个组数据筛选为88组可靠数据数据，数据如表1-1编程如下X1X2X3X4X5X6Y2102022222313623711692406016256116726013022700000280000031016773216451233271116213406910153518912203608131422370203338000004100000420000043800884410011453003346400444700000480000011041301712183964813134601297140426330151225811811607312381715203261801001210000022002122300424240031125003132600313270022228000003102224321512733011321434021133511227153604115370100138050354106218420116417430201323344014311745016184344601031134705217480311411021131203114130100110140211315140151601001170000018010012102103220686142305137182405741225012108222600969270210328032253100111320021233832113341223535741092136132263700101381010241002224201093194302414133844051064502617174346098817470424648043271、筛选数据程序%data(14,:)=;%data(28,:)=;%data(8,:)=;n=35;m=8;alpha=0.05;y=data(:,9);x1=data(:,1);x2=data(:,2);x3=data(:,3);x4=data(:,4);x5=data(:,5);x6=data(:,6);X=ones(n,1),x1,x2,x3,x4,x5,x6;b,bint,r,rint,s=regress(y,X,alpha);% b 回归系数% bint 回归系数的区间估计% r 残差% rint 残差置信区间% stats 用于检验回归模型的统计量，有三个数值：相关系数R2、F值、与F对应的概率p,相关系数R2越接近1，说明回归方程越显著；%F > F1-（k，n-k-1）时拒绝H0，F越大，说明回归方程越显著；与F对应的概率p 时拒绝H0，回归模型成立。% Y为n*1的矩阵；% X为（ones(n,1),x1,xm）的矩阵；% alpha显著性水平s2=sum(r.2)/(n-m-1);b,bint,s,s2rcoplot(r,rint); %用这个图来来做参差及其置信区间的图，如果数据的置信区间不包含零点，则可认为这个数据是异常的，应把它剔除利用第一个多元回归方程求出来的理论值1234567895612.35177.24392.85625.65133.74517.04552.8NAN4092.11011121314151617184521.45590.14991.25385.5NAN4273.75060.75389.95187.01920212223242526274665.65516.64964.05177.85883.24665.05083.34559.05763.5282930313233343536NAN5693.94801.05304.85392.14680.95378.56118.14876.937385245.05315.5利用最终的回归方程求出来的理论值1234567895396.85211.34469.55686.95266.84518.04558.8NAN3984.41011121314151617184530.95542.45072.95320.4NAN4187.65053.05270.05106.61920212223242526274939.95417.54890.65329.46012.64679.35136.44752.25618.2282930313233343536NAN5734.44688.45228.35412.64666.75430.95990.74914.137385260.75308.521