食品安全综合评价及预测.doc

食品安全综合评价及预测摘 要本文主要研究两个问题，其一是如何系统地对近年我国食品安全进行评价；其二是如何准确有效地预测未来几年我国食品安全情况。针对问题一：本文采用2000-2005年间食品化学残留合格率、微生物污染合格率、食品卫生总体合格率、食品安全认知度、食品出口产值、中毒人数、死亡人数以及食品安全法规数作为研究对象。由于影响食品安全的因子甚多，所以采用主成分分析法将多种因子进行降维处理，可得第一、第二、第三主成分的贡献率分别为56.16%、24.94%、11.97%，累计贡献率为93.07%，超过85%，故只需求出第一、第二、第三主成分z1、z2、z3、即可。其中化学残留合格率、微生物污染合格率、食品出口产值以及食品安全法规数占第一主成分的荷载分别为0.9903、0.9254、0.9926、0.9352，与第一主成分呈现明显的正相关，故此4类因子可归结为第一主成分。相似得，食品安全认知度归结为第二主成分，食品卫生总体合格率与死亡人数归结为第三主成分。接下来便是确定各类指标因子的权重，由前数据可知，第一主成分的4个因子主成分荷载类似，故权重系数可取56.16%4=0.1404；第二主成分只有一个因子，故权重等于贡献率，为0.2494；第三主成分有两个指标因子，根据主成分荷载分配贡献率可得权重分别为0.0731、0.0356。最后，便是根据各类因子数值及相应权重确定食品安全评价综合指数，其应为不同指标加权之和，F=0731X1+0.1404X2+0.1404X3+0.2494X4-0.0356X6+0.1404X7+0.1404X8，得出2000-2005年的食品安全评价指数分别为51.43、52.59、56.12、55.68、59.52、61.24。 针对问题二：本文采用2000-2005年已得出的每年食品安全评价指数，利用灰色系统理论，预测未来3年2011-2013年的评价指数值为65.4912、65.5376、66.0012，并进行检验，90%，精度为一级，故可以用此模型进行预测。由预测值可见，虽然未来几年食品安全指数逐年增加，但增加极为缓慢，食品安全仍然存在较大风险。关键词：食品安全 主成分分析 评价指数 灰色系统预测1.问题重述近年来，地沟油、硫磺姜、回炉面包、三鹿牛奶等食品安全事件频发，食品安全问题似乎是当下中国最热门的一幕丑剧。今年，瘦肉精事件尚未平息，4月，上海染色馒头又像一重磅炸弹引爆全国；随后接棒的则有温州毒馒头、广东墨汁粉条、辽宁毒豆芽等。常年风起云涌的食品安全事件，像一双双粗重的钢琴手，一次次肆虐地撩拨人们敏感的神经。为此，人大常委会建议：把食品安全作为国家安全组成部分，其重要性不亚于金融安全、粮食安全、能源安全、生态安全；加大对违反食品安全行为的惩处力度，以法治的方式维护公民的“食品安全权”；推进食品，尤其是农产品的规模化、产业化生产经验，从根本上保证食品质量；参照科学依据，结合中国国情，制定出全国统一的食品安全检验和认定标准；充分发挥社会组织作用，强化食品生产行业协会的作用，加强行业自律；加强宣传力度，将涉及到食品安全方面的法律法规汇编成册，针对从业者流动性强等特点，加大宣传力度并从严考核；充分发挥社会监督的作用，鼓励全社会共同参与食品安全监督等。请你们选择感兴趣的关于食品安全问题的某个侧面(比如：食品中污染物的分布，食品安全事件预警，您也可以选择其它方面研究)，建立数学模型，利用互联网数据，定量评估食品安全问题。 2.问题分析本课题的内容主要包括对食品安全问题的分析、评价及预测。由于食品安全是由多种因素共同作用产生的，存在着较大的模糊性和不确定性，因此要选取较为合理、代表性强、易于量化的评价指标来建立评价体系。本文采取食品卫生总体合格率、化学残留检测合格率、我国食品安全法规及标准健全程度等等为评价指标，确定评价指标的权重系数，由此而建立综合评价体系。而影响食品安全的因子甚多，鉴于这种情况，我们可以用主成分分析法把多种复杂因素限制为主要的几种代表因素，使之针对性较强，更易于研究对策。更为重要的一个问题是如何准确的预测未来食品安全情况，基于此系统为信息不完全明确的灰色系统，可以采用GM（1，1）模型进行比较准确地预测。3.模型假设1.假设所收集数据真实有效。2.假设食品安全评价主要由我们所选的指标确定，其他因素影响甚微，可以忽略。3.假设食品安全指数具有可预测性。4.符号说明F食品安全综合评价指数各评价指标各主成分第i个主成分的特征值特征向量主成分荷载灰色关联度5.模型的建立与求解5.1食品安全指标主成分分析5.1.1 主成分分析概述 利用主成分分析原理，先对指标值进行标准化处理，然后求出两两指标间的相关系数，得到相关系数矩阵R；求矩阵R的特征值、特征向量和方差贡献率，确定提取的主成分个数。5.1.2 食品安全综合评价指标 本文采用数据为2000-2005年统计数据。a 化学残留合格率与微生物污染合格率表一 我国2000-2005年食品化学残留及微生物检验合格率（%）年度食物化学残留检验合格率食品微生物合格率200087.1274.40200187.2374.40200288.9179.05200390.4984.97200491.0185.59200593.6286.58b 我国食品卫生总体合格率表二 我国2000-2005年食品卫生总体检验合格率（%）年度食品卫生监测总体合格率（%）200088.87200188.10200289.50200390.50200488.50200587.50c 我国2000-2005年食品出口产值表三 我国2000-2005年食品出口产值年度食品出口总产值备注2000122.8食物及活微生物2001127.82002146.22003175.32004188.72005242.0d.我国2000-2005年食品安全法规统计表表四 我国2000-2005年食品安全法规及标准统计表年份食品安全法标准200041602001416020024160200314360020041437292005405529e我国2000-2005安全管理认知度统计表表五 我国2000-2005年食品安全管理评价及食品安全认知度统计表年份食品安全认知度200039.6200139.6200339.6200338.5200439.6200539.6f.我国2000-2005年食物中毒统计表表六 我国2000-2005年食物中毒统计年份中毒人数死亡人数20006273150200113515130200213515130200310526330200414000283200590212355.1.3 主成分分析求解 （1）数据标准化利用所给数据，得到原始数据矩阵a 88.7 87.12 74.4 39.4 6273 151 122.8 488.1 87.23 74.4 38.7 13515 132 127.8 489.5 88.91 79.05 39.2 12547 137 146.2 490.5 90.49 84.97 38.5 10526 330 175.3 1488.9 91.01 85.59 39.6 14000 284 188.7 1487.5 93.62 86.58 39.8 9021 236 242 40 ，由于数据的不统一性，可以先进行标准化处理，数据标准化处理得到新矩阵方法：用matlab计算可得，（程序见附录）表七 矩阵标准化指标因子X1X2X3X4X5X6X7X8Z1-0.1581-1、0468-1.14160.3922-1.5832-0.7207-0.9870-0.6689Z2-0.7275-1.0027-1.1416-0.98060.8525-0.9464-0.8757-0.6689Z30.6010-0.3289-0.31620.00000.5269-0.8870-0.4660-0.6689Z41.54980.30480.7345-1.3728-0.15281.40580.18180.0478Z50.03160.51340.84460.78451.01560.85930.48010.0478Z6-0.12681.56021.02031.1767-0.65900.28911.66681.9113（2）计算相关系数（i，j=1，2，p）为原变量与的相关系数， =，其计算公式为：代入数据可得，例如： 用matlab计算可得，（程序见附录）表八相关系数矩阵Z1Z2Z3Z4Z5Z6Z7Z8Z11.0000-0.16750.1326-0.59550.11930.4201-0.2703-0.4425Z21.00000.93860.47860.01520.66320.99340.9224Z31.00000.30530.14810.84740.90270.7529Z41.0000-0.2452-0.05410.51200.5229Z51.00000.0668-0.0136-0.2164Z61.00000.62070.4739Z71.00000.9550Z81.0000（3）计算特征值与特征向量解特征方程 ，常用雅可比法求出特征值，并使其按大小顺序排列 ； 分别求出对应于特征值的特征向量 (i=1,2,p),要求 =1，即=1,其中表示向量的第j个分量。用matlab计算得，（程序见附录）表九 主成分与特征值主成分Z1Z2Z3Z4Z5Z6Z7Z8特征值4.49271.99480.95750.40760.1474000（4）计算主成分贡献率及累积贡献率贡献率：累计贡献率：注：一般取累计贡献率达8595%的特征值所对应的第一、第二、第m（mp）个主成分。 代入数据可求得， 特征值贡献率累计贡献率Z14.492756.1656.16Z21.994824.9481.10Z30.957511.9793.07可见，第一、第二、第三主成分累积贡献率已达到93.07%，超过85%，所以取3个主成分。5.1.4 指标因子权重计算主要成分荷载：计算可得如下主成分Z1Z2Z3X1-0.2183-0.87610.2971X20.9903-0.0523-0.0488X30.9254-0.3610-0.0426X40.53220.6688-0.0435X5-0.0500-0.4198-0.9025X60.6804-0.64290.1895X70.99260.0342-0.0598X80.93520.23790.0652其中化学残留合格率、微生物污染合格率、食品出口产值以及食品安全法规数占第一主成分Z1的荷载分别为0.9903、0.9254、0.9926、0.9352，与第一主成分呈现明显的正相关，故此4类因子可归结为第一主成分。根据相同算法，可得食品安全认知度归结为第二主成分，食品卫生总体合格率与死亡人数归结为第三主成分。然后便确定各类指标因子的权重，由前数据可知，第一主成分的4个因子主成分荷载类似，故权重系数可取56.16%4=0.1404；第二主成分只有一个因子，故权重等于贡献率，为0.2494；第三主成分有两个指标因子，根据主成分荷载分配贡献率可得权重分别为0.0731、0.0356。最后，便是根据各类因子数值及相应权重确定食品安全评价综合指数，其应为不同指标加权之和，F=0731X1+0.1404X2+0.1404X3+0.2494X4-0.0356X6+0.1404X7+0.1404X8，代入原始数据得出2000-2005年的食品安全评价指数如下表所示，表十 2000-2005综合指数值年份200020012002200320042005综合指数值51.4352.5956.1255.6859.5261.24整理成图，可见，2000-2005年间食品安全综合指数整体上呈上升趋势，但在2003年少许下降。5.2灰色系统理论预测与检测5.2.1 灰色系统模型建立灰色系统理论的微分方程成为Gm模型，G表示gray（灰色），m表示model（模型），Gm（1，1）表示1阶的、1个变量的微分方程模型。Gm（1，1）建模过程和机理如下：记原始数据序列为非负序列其中，其相应的生成数据序列为其中，为的紧邻均值生成序列其中，称为Gm(1,1)模型，其中，b是需要通过建模求解的参数，若为参数列，且， 则求微分方程的最小二乘估计系数列，满足 称为灰微分方程，的白化方程，也叫影子方程。如上所述，则有1.白化方程的解或称时间响应函数为2.Gm(1,1)灰微分方程的时间响应序列为3.取，则4.还原值首先运用主成份分析法对2000年到2005年的食品卫生监测合格率、不同食品化检测合格率、不同食品生产和消费总量、食品化学残留检测合格率、食品安全管理评价及食品安全认知度统计表、食物中毒统计、食品安全法规及标准统计表、食品出口总值八个指标的综合分析评估，分别计算出我国2000年到2005年的食品安全综合指数估算值。见表十，代入数据可得，= =（51.43，104.38，159.29，216.24，275.28，336.52）还原出的模拟值，由， 得 =（51.43，52.96，54.92，56.95，59.05，61.24）5.2.2灰色系统模型的检验=37.81=37.78=0.03 =0.990.9精度为一级，可以用 进行预测5.2.3灰色系统模型预测2011-2013年预测值为 =（65.49，65.57，66.00）可见，虽然未来3年食品安全指数呈上升趋势，但是上升速度缓慢，我国食品安全仍然面临着较大的风险。6.模型的优缺点评价6.1优点：1. 运用主成份分析法对涉及到食品安全的各种因素进行了综合考虑筛选，并且引入了权重，使得各种因素对食品安全的综合评价影响充分体现出来。2. 运用matlab对数据进行合理的处理计算，是计算数据和计算过程更加清晰全面地展现出来。3. 用灰色系统预测模型对食品安全问题未来几年的预测，同时对模型的精度进行检核，是预测结果更具有可信性。6.2缺点：1. 权重的引入具有一定的主观因素，可能是模型计算结果受到一定影响。2. 由于对数据的掌握有限，所以使得对食品安全指数的综合评估可能受到影响。7.参考文献【1】中国卫生检验杂志，2006，16（1）：117.【2】中国卫生部公告2000-2005.【3】中国统计年鉴2000-2005.【4】刘思峰.灰色系统理论及其应用（第三版） 社会科学出版社 2004年.8.附录程序1 a=88.7 87.12 74.4 39.4 6273 151 122.8 488.1 87.23 74.4 38.7 13515 132 127.8 489.5 88.91 79.05 39.2 12547 137 146.2 490.5 90.49 84.97 38.5 10526 330 175.3 1488.9 91.01 85.59 39.6 14000 284 188.7 1487.5 93.62 86.58 39.8 9021 236 242 40%原始数据矩阵aa = 1.0e+004 * 0.0089 0.0087 0.0074 0.0039 0.6273 0.0151 0.0123 0.0004 0.0088 0.0087 0.0074 0.0039 1.3515 0.0132 0.0128 0.0004 0.0089 0.0089 0.0079 0.0039 1.2547 0.0137 0.0146 0.0004 0.0091 0.0090 0.0085 0.0039 1.0526 0.0330 0.0175 0.0014 0.0089 0.0091 0.0086 0.0040 1.4000 0.0284 0.0189 0.00140.0088 0.0094 0.0087 0.0040 0.9021 0.0236 0.0242 0.0040 b=zscore(a)%标准化后矩阵bb = -0.1581 -1.0468 -1.1416 0.3922 -1.5832 -0.7207 -0.9870 -0.6689 -0.7275 -1.0027 -1.1416 -0.9806 0.8525 -0.9464 -0.8757 -0.6689 0.6010 -0.3289 -0.3162 0.0000 0.5269 -0.8870 -0.4660 -0.6689 1.5498 0.3048 0.7345 -1.3728 -0.1528 1.4058 0.1818 0.0478 0.0316 0.5134 0.8446 0.7845 1.0156 0.8593 0.4801 0.0478 -1.2968 1.5602 1.0203 1.1767 -0.6590 0.2891 1.6668 1.9113 r,p = corrcoef(b)%相关系数矩阵rr = 1.0000 -0.1675 0.1326 -0.5955 0.1193 0.4201 -0.2703 -0.4425 -0.1675 1.0000 0.9386 0.4786 0.0152 0.6632 0.9934 0.9224 0.1326 0.9386 1.0000 0.3053 0.1481 0.8474 0.9027 0.7529 -0.5955 0.4786 0.3053 1.0000 -0.2452 -0.0541 0.5120 0.5229 0.1193 0.0152 0.1481 -0.2452 1.0000 0.0668 -0.0136 -0.2164 0.4201 0.6632 0.8474 -0.0541 0.0668 1.0000 0.6207 0.4739 -0.2703 0.9934 0.9027 0.5120 -0.0136 0.6207 1.0000 0.9550 -0.4425 0.9224 0.7529 0.5229 -0.2164 0.4739 0.9550 1.0000p = 1.0000 0.7512 0.8023 0.2124 0.8218 0.4069 0.6044 0.3796 0.7512 1.0000 0.0055 0.3370 0.9772 0.1511 0.0001 0.0088 0.8023 0.0055 1.0000 0.5563 0.7794 0.0332 0.0137 0.0841 0.2124 0.3370 0.5563 1.0000 0.6396 0.9190 0.2991 0.2871 0.8218 0.9772 0.7794 0.6396 1.0000 0.8999 0.9796 0.6805 0.4069 0.1511 0.0332 0.9190 0.8999 1.0000 0.1885 0.3424 0.6044 0.0001 0.0137 0.2991 0.9796 0.1885 1.0000 0.00300.3796 0.0088 0.0841 0.2871 0.6805 0.3424 0.0030 1.0000 v,d=eig(r)%特征值向量矩阵v,特征值dv = -0.1182 0.1211 0.3051 0.4909 0.3869 0.3036 -0.6203 -0.1030 0.1248 0.6291 -0.5214 0.3060 -0.0321 -0.0499 -0.0373 0.4672 0.5141 -0.6374 -0.1998 0.0725 0.1626 -0.0435 -0.2556 0.4366 -0.0274 0.0965 0.1923 -0.1077 0.8075 -0.0445 0.4735 0.2511