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

2013高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： B 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 贵州民族大学 参赛队员 (打印并签名) ：1. 沈兴富 2. 苏小钰 3. 冯 彬 指导教师或指导教师组负责人 (打印并签名)： 储昌木 （论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。） 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2013高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：食品质量安全抽检数据分析摘要本文结合深圳市2010年到2012年食品抽检数据, 通过分析各类抽检食品不合格率的变化关系,建立了关于食品安全的综合评价模型. 随后, 对数据进行k-均值聚类, 利用因子分析建立反映食品安全与各因素之间变化规律的模型, 最后针对食品抽检有效性, 提出一套科学有效的抽检方案. 针对问题一, 首先统计食品中含微生物, 重金属, 添加剂和其它有害物质的数量, 计算出不合格率, 再利用时间序列ARIMA模型分别给出微生物, 重金属, 添加剂和其它有害物质的变化趋势, 并建立综合评价矩阵, 给出综合评价指数模型. 最后得出深圳市近三年来食品安全指数的变化呈下降趋势结论, 即深圳市这三年来食品安全质量总体提高了. 针对问题二, 首先对数据进行k-均值聚类, 再利用因子分析的方法来判断季节、生产环节、食品包装、抽检地点、食品种类等因素对食品安全的影响. 结果显示食品种类对食品安全影响最大. 另外, 对这五个因素的食品质量相关性利用卡方统计量进行一致性检验, 得出食品种类的变化是很显著的结论. 此外, 对食品产地与食品质量之间的相关关系利用pearson积距相关系数进行判别, 得出食品运输距离的长短与食品合格率相关性较明显, 运输距离越长越容易出现食品安全问题的结论.针对问题三, 通过对这三年的抽检情况进行分析, 发现每年抽检未能利用上一年数据的有效性, 针对这一情况我们提出, 对于前两年抽检不合格的地方所抽检到不合格的产品未来一年一定进行抽检, 而对于一般情况则采用分层抽样的方法进行抽检. 最终达到更科学、更有效地反映食品质量状况, 减少抽检成本的目的. 关键词 食品安全综合评价 K-均值聚类 因子分析 分层抽样1 问题重述近年来, 我国食品安全领域出现了令人忧虑的问题. 肯德基的“苏丹红”、豆腐中的“吊白块”、水饺中的“霉青菜”更危险的是“三聚氰胺”, 它不仅在牛奶中大量出现, 还出现在鸡蛋中. 这些形形色色的食品安全问题, 就像日前全球爆发的经融危机一样席卷整个大地, 给人们的生命将抗带来了严重影响, 更牵动着大家的心. 城市食品的来源越来越广泛, 人们消费加工好的食品的比例也越来越高, 因此除食材的生产收获外, 食品的加工、包装、贮存、运输、销售等每一个环节皆可能影响食品的质量与安全. 另一方面, 食品质量与安全又是一个专业性很强的问题, 其标准的制定和抽样检测及评价都需要科学有效的方法. 本次研究以2010至2012年深圳市的食品抽检数据为基础, 从以下几个方面对食品质量安全进行分析：(1)找出深圳市这三年不合格食品中微生物、重金属、添加剂、和其它有害物质超标对于整个食品产业的影响及食品安全的变化趋势. (2)利用统计数据, 探索食品产地、抽检地点、季节因素、食品种类等与食品安全质量的关系. (3)为了加强深圳市的食品抽查有效性, 改进食品抽检的办法, 使之更科学、更有效地反映食品质量状况且不过分增加监管成本(食品抽检是需要费用的), 提出一个对抽检结果稳定且抽检频次过高的食品领域调整方案. 2 问题分析2.1 问题一的分析考虑到深圳食品安全问题主要受微生物、重金属、食品添加剂这些因素的影响, 我们统计2010至2012年深圳市的食品抽检数据后, 将数据预处理并且进行对微生物、重金属、添加剂、其它有害物质的预测. 运用食品安全综合评价的阶梯层次结构模型, 计算出2010至2012三年的食品安全综合评价指数, 最终利用综合评价指数评价这三年各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势. 2.2 问题二的分析为了反映规律性的结果, 首先对现有的统计数据进行分析, 归纳出季节、生产环节、食品包装、抽检地点、食品种类、食品产地等因素对食品质量可能存在的影响. 再用因子分析方法分析哪些因素对于食品质量有较大的影响. 最后用卡方统计量方法对不可量化的因素(如季节)进行判断、pearson积距相关系数方法对可量化因素(如食品产地可用距离量化)进行判断, 从而总结出这些因素与食品质量的关系. 2.3 问题三的分析注意到现有的抽样过程未能有效利用以前的抽检结果, 进而不会减少抽检成本. 如果能充分利用以前的抽检数据, 减少没必要的抽检, 就会大大减少人力物力的开销. 对于前两年抽检不合格的地方所抽检不合格的产品未来一年一定进行抽检, 而对于前两年合格的产品采用分层抽样的方法进行抽检. 3 模型假设(1)假设食品的安全仅与有害物质的不合格率相关. (2)假设深圳市2010至2012年深圳市的食品抽检数据客观真实, 且本文所研究的有毒有害物质仅涉及食品中微生物、重金属、添加剂及其它有害物质. (3)假设食品质量与食品种类、季节、生产地点、抽检地点等有关. (4)假设每次抽检食品费用相同. 4 符号说明及名词定义第季度(=1, 2, 3, 4)食品安全的综合评价指数一致性指标平均随机一致性指标一致性比例指标V-A层判断矩阵A的元素A-B层判断矩阵B的元素判断矩阵的最大特征值(=1, 2, 3, 4)判断矩阵的阶数判断矩阵的归一化特征相量 判断矩阵的归一化特征相量(=1, 2, 3, 4)误差平方和主成分分析法的指标变量第个评价对象的第个指标取值(=1, 2, 3, 4; =1, 2, 3, 4)指标值转化的标准化指标值食品安全的综合评价值5 模型建立与求解5.1 问题一的模型与建立5.1.1 问题一的数据预处理通过对问题分析, 我们了解到本题是一个多元评价问题. 为评价各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势. 需要将数据进行归类和统计1, 先将每一年的数据归总, 归总数据包括不合格项和合格项, 其中不合格项包括样品名称、所属环节(生产、流通、餐饮), 文字商标, 型号规格等级, 受检单位及受检单位详细地址, 生产单位及生产单位详细地址, 生产日期, 抽检日期, 样品分类(如附录7所示28类), 不合格项目及实测结果, 不合格项目所属分类(微生物、重金属、添加剂及其它)；合格项包括样品名称, 所属环节,文字商标, 型号规格等级, 受检单位及受检单位详细地址, 生产单位及生产单位详细地址, 生产日期, 抽检日期, 样品分类处理结果见附件一. 再统计出每年每季中不同指标的不合格项数和合格项数. 为了预测未来一年内的各个指标的不合格率，我们利用时间序列ARIMA序列与季节性序列进行预测2, 预测结果见图1(程序见附录1) .图1 2013年各指标不合格率预测趋势未来一年的各个微生物, 重金属, 添加剂等指标的不合格率相对于前三年有降低且稳定的趋势.5.1.2 问题一模型建立 (1) 构造食品安全综合评价的的阶梯层次结构3, 4, 5, 6. 在食品安全综合评价时, 考虑微生物、重金属、添加剂、其它有害物质等4项指标的综合效果, 将食品安全综合评价指数作为目标层, 记为, 将第一、二、三、四季度作为准则层, 将上诉四项指标作为指标层, 构造其综合评价的阶梯层次结构如图2.食品安全综合评价第一季度第二季度第三季度第四季度微生物重金属添加剂其它图2 食品安全综合评价的阶梯层次结构图(2) 建立层判断矩阵建立层判断矩阵的依据是：当认为每个季度的重要性相同时, 将每个季度都定为一级, 可得层判断矩阵. , 其中 .(3) 层判断矩阵该层的判断矩阵元素确定如下：将每个季度,的各指标用分数公式 (1)标准化, 其中是在第季度下由第个指标导致产品不合格的产品总数, 是当年所测的所有产品总数. 对于第个季度下, 取 设 (2) 其中,分别为矩阵第, 列；; . 对于任意的与相比, 其判断矩阵元素由以下确定 (3)按式(1)可得各指标观测值的标准化值, 再由式子(2), (3)得层的判断矩阵.(4) 用方根法计算出判断矩阵, 的归一化特征向量:,(=1,2,3,4).(5) 计算判断矩阵的最大特征值(=1, 2, 3, 4). 并进行一致性检验 (4)其中, n为判断矩阵的阶数计算一致性比例其中为平均一致性指标, 当时判断矩阵具有好的一致性. (6) 食品安全质量的综合评价指数7的计算 (5)其中为的分量, 为的分量(=1, 2, 3, 4; =1, 2, 3, 4).5.1.3 问题一模型求解(1) 以深圳市2011年食品安全监察所得的数据为例, 求该年食品安全的综合评价指数. 根据所处理的数据, 由于该年的数据因某些原因不存在, 则取其值为其它因子的平均值. 并按式(1)可得对应的标准化值如表1所示, 并将第四季度作为第一级, 其它三个季度作为第二级. 表1 各指标的标准化值因子微生物重金属添加剂其它A10. 03084550. 02477230. 01853920. 0073517A20. 01548090. 004611330. 02865610. 0111989A30. 02835610. 005025130. 02422830. 0129218A40. 02489410. 01146960. 02380790. 0104908由公式(2), (3)得层的判断矩阵分别为其中代表第个指标重要程度与第个指标重要程度的比值, =1, 2, 3, 4.(2) 用方根法计算出判断矩阵,(=1, 2, 3, 4)的归一化特征向量其中的4个分量分别表示层的4个因子(第一季度、第二季度、第三季度、第四季度)在综合评价中的各自的权重. 4个分量分别只表示考虑的条件下4个指标( 微生物、重金属、添加剂、其它)各自应占的权重, , , 类推. (3) 计算判断矩阵的最大特征值并进行一致性检验. 求矩阵的最大特征值, 在matlab中用dig命令能容易求出. 得计算一致性比例, 其中, 为平均一致性指标, 为判断矩阵的阶数. 查表知, 故 当时判断矩阵具有好的一致性. (4) 食品安全的综合评价指数的计算由此求得深圳市这三年食品安全综合评价指数(见表2), 求解程序见附录2,详细结果见附录3. 表2 深圳市这三年食品安全综合评价年份20100. 01276810. 00040. 01190. 02020. 000620110. 01540260. 02710. 02390. 02170. 000620120. 00588570. 02240. 00070. 03130. 0199由表2知, 2010年到2012年的食品安全综合评价指数先微升后大降. 根据食品综合评价越低越安全的准则, 可推断出这三年各主要食品领域微生物、重金属、添加剂含量等安全食品安全质量水平整体呈现出更好的趋势, 且趋势显著. 5.2 问题二的模型与求解5.2.1 问题二数据的预处理 经过统计2010到2012年深圳市的食品抽检数据, 先将数据(附件一)按季节、生产环节等的不合格率利用K-均值聚类算法聚类8. 其次, 利用因子分析法来判断食品产地与食品质量的关系、食品销售地点(即抽检地点)与食品质量的关系等. 我们利用spss软件9将各个因素进行K-均值聚类后的最终聚类中心数据表3. 表3 统计数据聚类结果聚类第一类第二类第三类季节0. 0807370. 020490. 044036生产环节0. 0071230. 0494940. 064242食品种类0. 0897820. 0159540. 143832食品产地0. 0662390. 1512350. 033724抽检地点0. 0882350. 0337240. 05934该表为每种因素k-均值聚类后的结果, 其中迭代历史记录, 聚类成员, 最终聚类中心间的距离, 方差分析(ANOVA)和每个聚类中的案例数在附件二中, 用这些最终的三个聚类结果作为每种因素聚类后的值. 5.2.2因子分析模型模型的建立(1) 进行因子分析2的指标变量有3个, 分别为, 共有5个评价对象, 对第个评价对象的第个指标取值为, ; . 将各指标值转化为标准化指标值, 有 (6)其中 , 即,为第个指标的样本均值和样本标准差. 相应地, 称 (7)为标准化指标变量. (2) 计算相关系数矩阵R计算相关系数矩阵, 有 (8) 其中, 是第个指标与第个指标的相关系数, 满足,. (3) 计算初等载荷矩阵计算相关系数矩阵的特征值, 及对应的特征向量, 其中, 初等载荷矩阵 (9)(4) 选择3个主因子根据初等荷载矩阵, 计算各个公共因子的贡献率, 并选择个主因子. 对提取的因子载荷矩阵进行旋转, 得到矩阵(其中为的前3列，为正交矩阵)构造因子模型 (10)(5)计算因子得分, 并进行综合评价用回归方法求单个因子得分函数 (11)记第个样本点对第个因子得分的估计值则有 (12) 且 (13)其中是原始数据的标准化数据矩阵; 为相关系数矩阵. 模型求解(1) 我们选取两个主因子, 第一个公因子为第一类聚类, 第二个公因子为第二类聚类. 利用Matlab程序计算得到旋转后的因子贡献及贡献率, 见表4. 表4 因子的贡献率和旋转因子分析贡献率数据旋转因子分析表因子贡献贡献率累计贡献率指标主因子1主因子211. 575452. 512952. 5129第一种聚类0. 0879550. 98900621. 030534. 351586. 8644第二种聚类-0. 907430. 0366第三种聚类0. 8626840. 225992(2) 计算得到各个因子得分函数为 (14)计算出5种因素对食品安全的综合得分见表5表5 各个因素综合排名表排名12345F10. 058080. 003041-0. 021290. 003601-0. 08004F20. 1042650. 087090. 105380. 0217810. 09561F0. 0763450. 0362790. 0288030. 010791-0. 01058因素食品种类季节抽点起点生产环节食品产地通过相关分析, 得出各个因素贡献能力与综合食品安全指数的相关系数为0. 7822, 这表明两者存在高度相关关系. 因子分析法的回归方程为, 回归方程在显著水平0. 05的情况下, 通过了假设检验. 5.2.3 卡方统计模型为了反映规律性的结果,首先需要对现有的统计信息进行定量分析. 由于问题二的求解几乎需要考虑文中提供的所有已知信息, 并且网络数据中列表所反映情况具有差异性(有的是按食品的加工, 包装, 贮存, 运输, 销售抽检. 有的按时间, 地域抽检). 这里对不同的因素进行不同的相关关系的判断. 此处涉及季节与食品合格率, 生产环节与食品合格率, 食品包装与食品合格率, 抽检地点(根据深圳有六个不同地区, 分为了六个抽检地点)与食品合格率, 食品种类(食品种类有28种,本文应用k-均值聚类方法, 将食品种类聚类成6类)与食品合格率, 食品产地(距离不同)与食品合格率, 以上六种相关关系的判定不能只由一种方法进行简单判断, 对此选择季节与食品合格率、生产环节与食品合格率、食品包装与食品合格率、抽检地点与食品合格率、食品种类与食品合格率这五项进行卡方统计量判别. 对食品产地与食品合格率这一项之间的相关关系进行pearson积距相关系数判别. 卡方统计量判别方法如下：按如下格式, 假设有两个集合, . =事件, 事件; =合格数, 不合格数.总计aba+bcdc+d总计a+cb+da+b+c+d其中, a代表事件中合格数的个数; b代表事件中不合格数的个数; c代表事件中合格数的个数; d代表事件中不合格数的个数.由独立性检验公式 (15)当都不小于5时, 可通过查阅表6来确定结论有关的可信程度.表6 卡方概率表P(K2k)0. 500. 400. 250. 150. 10k0. 4550. 7081. 3232. 0722. 706P(K2k)0. 050. 0250. 0100. 0050. 001k3. 8415. 0246. 6357. 87910. 828(1) 计算季节与食品合格率之间的卡方统计量.根据表7数据进行卡方统计量判别, 结果见表8.表7 2010到2012年各个季度的数据统计合格数不合格数合格率不合格率一季度6581. 00 543. 00 0. 9341650. 065835二季度8959. 00 379. 00 0. 9620030. 037997三季度12708. 00 726. 00 0. 9545730. 045427四季度36692. 00 1903. 00 0. 9619640. 038036表8 季节与食品合格率之间的卡方统计量一季度二季度三季度四季度一季度097. 054739. 539686. 0290二季度021. 601112. 6093三季度04. 6571四季度0分析表8数据, 可得季节与食品合格率之间有很大可能相关, 且相对第一季度到第二季度对食品合格率的影响最大, 应多加注意. (2) 生产环节与食品合格率之间的卡方统计量(以2011年的数据为例). 根据表9数据进行卡方统计量判别, 结果见表10. 表9 流通领域数据统计合格数不合格数合格率不合格率生产环节2091150. 9928770. 007123流通环节1786930. 9505060. 049494餐饮环节27531890. 9357580. 064242表10 流通领域卡方统计量生产环节流通环节餐饮环节生产环节067. 6161103. 2654流通环节04. 5284餐饮环节0由表10数据知, 生产环节与食品合格率之间也有相关性. 且由表9知餐饮环节的食品不合格率的最大. 因此相对其它环节, 应多注意餐饮环节, 其次是流通环节. (3) 食品包装与食品合格率之间的卡方统计量(以2011年的数据为例).根据表11数据进行卡方统计量判别, 可得. 得出食品包装与食品合格率之间同样有相关性. 且袋装食品的不合格率较大. 因此相对散装食品, 应多注意袋装食品. 表11 食品包装数据统计合格数不合格数合格率不合格率散装51283590. 9345730. 065427袋装83167490. 9173750. 082625(4) 抽检地点与食品合格率之间的卡方统计量. 本文按照深圳市有六个区, 抽检地点分为六个地点(以2011年的数据为例). 根据表12数据进行卡方统计量判别, 结果见表13. 表12 抽检地点数据统计合格数不合格数合格率不合格率罗湖区655440. 9370530. 062947福田区572470. 9240710. 075929南山区822440. 9491920. 050808盐田区305230. 9298780. 070122宝安区19221320. 9357350. 064265龙岗区14261380. 9117650. 088235 表13 抽检地点与食品合格率之间的卡方统计量罗湖区福田区南山区盐田区宝安区龙岗区罗湖区00. 86071. 07400. 18850. 01514. 1773福田区03. 96000. 10561. 03570. 8659南山区01. 67541. 947511. 2684盐田区00. 15951. 1426宝安区07. 3877龙岗区0由表13知抽检地点不同与食品合格率之间同样有相关性. 但相对于其它因素对食品合格率的影响, 抽检地点对食品合格率影响较不明显. 从表12可知, 福田区抽检的食品不合格率最高, 为0. 075929；南山区抽检的食品不合格率最低, 且为0. 050808. (5) 食品种类与食品合格率之间的卡方统计量(2010、2011、2012三年所有数据). 按照国家级食品分类标准将食品种类分为28种, 统计的结果见附录5. 通过k-均值聚类, 将28种食品种类聚类成6类. 具体为第1类：豆制品; 第2类：包装类, 焙烤咖啡及其制品, 炒咖啡制品, 蛋制品, 乳制品, 薯条和膨化食品, 速冻食品, 饮料;第3类：调味品, 糕点, 其他食品, 肉制品;第4类：饼干, 罐头, 冷冻饮品, 蜜饯及水果制品, 水产品;第5类：蜂产品, 酒类;第6类：茶叶, 炒货食品及坚果制品, 方便食品, 可可及焙, 粮食加工品, 蔬菜制品, 糖果制品, 特殊膳食食品;所聚类的结果见附录5. 最终统计为表14. 表14 食品种类数据统计类别合格总数不合格总数合格率不合格率17441440. 8378380. 16216226640570. 9914890. 008511369177970. 8966810. 103319428892240. 9280440. 071956510521570. 8701410. 1298596100112930. 9715640. 028436再由卡方统计量判别, 结果见表15. 表15 食品种类与食品合格率之间的卡方统计量12345610717. 515428. 301867. 31584. 3455389. 646120577. 9853307. 4469572. 633179. 91313025. 54267. 7106435. 23854036. 3214122. 223050296. 357060由表15知食品种类的变化是很显著, 说明食品种类的不同决定着有不同程度的不合格率. 且由表14可知第一类(即豆制品)的不合格率是最高的, 为0. 162162；第二类的不合格率是最低的, 为0. 008511. 由高到低依次为第一类、第五类、第三类、第四类、第六类、第二类. (6) 食品产地与食品不合格率之间使用pearson积矩相关系数判别法进行相关性的判断.此处以2011年食品运输距离长短与食品不合格率进行相关性判别, 绘表16.表16 食品产地的食品运输距离长短与不食品合格率. BJHNDWFSFJGZGXSZCQSTX(km)216080769. 414882514354301614334Y(%)13.5816. 662. 522. 893. 882. 514. 683. 725. 557. 69符号代表如下：X: 与深圳市距离 Y: 产品不合格率BJ(北京) HN(湖南) DW(东莞) FS(佛山) FJ(福建) GZ(广州) GX(广西) SZ(深圳) CQ(重庆) ST(汕头) 由pearson积矩相关系数公式 (16)可对表16求解, 得：这说明食品运输距离的长短与食品合格率相关性较明显, 运输距离越长, 越容易出现食品安全问题. 5.3 问题三的模型与求解5.3.1模型的建立为了弄清节假日和非节假日食品安全程度(不合格率)是否存在显著性差异, 对2010年到2012年的数据进行单因素方差分析(程序见附录5), 结果得到故节假日和非节假日食品安全没有显著差异.根据2010到2012年的数据中出现有严重食品安全问题的食品生产经营单位或风险较高的食品, 在未来一年中都要进行相应的产品抽查. 如某地点的产品不合格数比较多就适当加大该地点产品的抽检频次或抽检批次. 而对于一般情况的食品生产经营单位或食品采用分层抽样进行抽检 10. 设总体分为层, 以表示层的编号, 第层的单元数为, 抽检样本数为. 又(为总体)已知, , 层权也是已知的. 分别表示层总和样本的第个单位的指标值. 第层总体(样本)均值分别为 (17) (18)第层总量及样本总量分别为 (19), (20)第层方差及样本分别为 (21), (22)第层抽样比为, (23)样本均值是总体均值的无偏估计, (24)总体估计量的方差为, (25)第层的估计量方差, (26)由误差知识可得 (27)联立(25)和(27)各式可得： (28)其中d为误差，11.本文考虑简单线性费用函数,总费用为, (29) 其中是与单元抽取量无关的费用，是在第层中抽取一个单元的平均费用.各层样本量的最有分配为, (30)由于每个样本的调查费用相同，故 , (31)联立(28)和(31)可得, 总体抽样率为, (32)分层抽样率为12. (33)5.3.1问题三模型的求解 把季度抽检分为4层, 抽检地区(6个区)在前一种基础上再分为6层, 流通区域在前一种基础上在分为3层, 食品种类(六类)在前一种基础上在分为5层如图3 合格产品抽检树状图第一季度第二季度第三季度第四季度罗湖区 福田区区 南山区 盐田区宝安区 龙岗区 生产环节 流通环节 餐饮环节 豆制品 包装类 焙烤咖啡及其制品炒咖啡制品蛋制品乳制品薯类和膨化食品速冻食品饮料 调味品糕点其他食品肉制品 饼干罐头冷冻饮品蜜饯及水果制品水产品 蜂产品酒类 茶叶炒货食品及坚果制品方便食品可可及焙粮食加工品蔬菜制品糖果制品特殊膳食食品 图3 深圳市往年合格产品抽检树状图其中第四层的抽检单位分别为, 第一种聚类(豆制品), 第二种聚类(包装类、焙烤咖啡及其制品、炒咖啡制品、蛋制品、乳制品、薯条和膨化食品、速冻食品、饮料), 第三种聚类(调味品、糕点、其他食品、肉制品), 第四种聚类(饼干、罐头、冷冻饮品、蜜饯及水果制品、水产品), 第五种聚类(蜂产品、酒类), 第六种聚类(茶叶、炒货食品及坚果制品、方便食品、可可及焙、粮食加工品、蔬菜制品、糖果制品、特殊膳食食品). 再通过(程序见附录6)计算得表17表17 多层分成抽样权重季度第一季度第二季度第三季度第四季度权重41.47%26.30%26.70%5.53%地区罗湖区福田区南山区盐田区宝安区龙岗区权重15. 27%18. 42%12. 32%17. 01%15. 59%21. 40%生产环节生产环节流通环节餐饮环节权重5. 89%40. 95%53. 15%食品种类第一种聚类第二种聚类第三种聚类第四种聚类第五种聚类第六种聚类权重16. 22%0. 85%10. 33%7. 20%12. 99%2. 84%由表17可以科学的分配下一年的抽检工作(且每次抽检都是突击抽检以防作假), 很大力度上的节约人力和物力. 例如在第一季度中, 罗湖区的抽检比例占全年的15. 27%, 福田区和南山区各占18. 42%和12. 32%. 并且第一种食品聚类(豆制品)抽检占16. 22%, 第二种食品聚类(包装类、焙烤咖啡及其制品、炒咖啡制品、蛋制品、乳制品、薯条和膨化食品、速冻食品、饮料)占0. 85%等等. 由于节假日和非节假日食品安全没有显著差异, 我们提出在上述模型的基础上随机时间、随机地点进行检查. 即对前两年抽检不合格的地方所抽检不合格的产品未来一年一定进行抽检, 对前两年合格的产品采用分层抽样的方法进行抽检以达到有效性的效果；对上述需要检查的所有内容, 用随机日期和随机地点进行抽检以防止商家造假以应付抽检, 从而达到科学性的效果. 7 模型评价7.1 模型的优点(1) 通过综合评价的方法, 解决了问题中对深圳市2010到2012各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势综合评价问题. (2) 通过因子分析很好的判断各个因素对食品安全的影响大小. (3) 分层抽样的方法, 它对不同层次不同批次的产品进行抽检, 建立了样本分配率、样本方差、总体抽样率、分层抽样率等函数方程. 对任何特征分层比较明显的系统均可以采用这种方法进行抽样. 而且这种可靠性较高. (4) 模型的结构简单便于推广和改进, 对现实具有很强的指导意义. (5) 对深圳市食品抽检问题提出了更科学, 更有效切实可行的抽检方案. 7.2 模型的缺点(1) 采用 Pearson 积矩相关系数判别法, 它判断的结果是比较粗糙的, 不适于精度要求很高的问题.(2) 由于所给定的数据有限, 对个别问题的解决难免带来了误差. 参考文献:1 深圳市市场监督管理局网站http:/www. szaic. gov. cn，2013. 6. 8.2 司守奎, 孙玺菁, 数学建模算法与应用M, 北京：国防工业出版社, 2012. 3 李祚泳等, 环境质量评价原理与方法M, 北京：化学工业出版社, 2004.4 唐苏文, 马守贵, 张浩, 食品卫生安全保障体系数学模型的研究J, 数学的实践与认识, 38（14）：114-119, 2008. 5 武力, 食品安全风险评价研究J, 食品工业科技, 12（9）：45-53, 2010. 6 陈锦屏, 张志国, 关于影响食品安全因素的探讨J. 食品科学, 2005年26卷第8期. 7 李玉鑑, K均值聚类算法, http:/www. docin. com/p-436874252. html，2013. 6. 8. 8 卢纹岱, SPSS统计分析M, 北京：电子工业出版社, 2003. 9 L. Kish , 倪加勋, 抽样调查M, 北京：中国统计出版社, 1985. 附录附录1时间序列程序clc,clearload a; %把原始数据按照原来的排列格式存放在纯文本文件a.txta=a1,:;%第一列为微生物的不合格率，二列为重金属，三列为添加剂，四列为其它a=nonzeros(a); %按照原来数据的顺序去掉零元素r11=autocorr(a) %计算自相关函数r12=parcorr(a) %计算偏相关函数da=diff(a); %计算 1 阶差分r21=autocorr(da) %计算自相关函数r22=parcorr(da) %计算偏相关函数n=length(da); %计算差分后的数据个数for i=0:3for j=0:3spec= garchset(R,i,M,j,Display,off); %指定模型的结构coeffX,errorsX,LLFX = garchfit(spec,da); %拟合参数num=garchcount(coeffX); %计算拟合参数的个数%compute Akaike and Bayesian Information Criteriaaic,bic=aicbic(LLFX,num,n);fprintf(R=%d,M=%d,AIC=%f,BIC=%fn,i,j,aic,bic); %显示计算结果endendr=input(输入阶数 R);m=input(输入阶数 M);spec2= garchset(R,r,M,m,Display,off); %指定模型的结构coeffX,errorsX,LLFX = garchfit(spec2,da) %拟合参数sigmaForecast,w_Forecast = garchpred(coeffX,da,4) %计算 4步预报值x_pred=a(en