食品卫生安全保障体系数学模型的研究.pdf

第38卷第14期 2008年7月 数学的实践与认识 MA THEMA T ICS I N PRACT ICE AND THEORY Vol138 No114 July 2008 食品卫生安全保障体系数学模型的研究 唐苏文 马守贵 张 浩 东南大学 移动通信国家重点实验室 江苏 南京 210096 摘要 主要针对食品安全保障体系中的食物摄入量模型 污染物分布模型和风险评估模型做了以下研究 提出一种基于多层次划分法的抽样方案设计抽样 针对食品的分类问题 提出一种基于贡献率的食品分类 法 提出利用广义g h分布估计污染物的偏态分布的方法 进一步基于Bootstrap法提出了一种利用少量数 据再抽样估计偏态分布特性的方法 提出一种尾部放大法 巧妙地将概率搜索和理论分析相结合 实现高精 度地估计污染物摄入量的右分位点 此外 对相关的几个理论问题进行研究 并基于Monte2Carlo模拟给出相 应的解决方案 关键词 多层次划分法 基于贡献率的食品分类法 广义g h分布 Bootstrap法 尾部放大法 Monte2 Carlo模拟 1 引 言 收稿日期 2008204201 请参见2007年第四届全国研究生数学建模竞赛A题 建立食品卫生安全保障体系数学 模型及改进模型的若干理论问题 2 模型假设及参数说明 经过分析 我们发现这是一真刀真枪的实际问题 为了使问题得到很好的解决 特作一 些假设和符号说明 限于篇幅 这里不再给出 下文用到时再作说明 3 人群食物摄入量模型 人群食物摄入量模型是用于估计不同地区 不同性别 不同年龄 不同季节 不同劳动强 度 不同经济收入的人群各类食品的一天摄入量 不同地区 不同性别 不同年龄 不同季节 不同劳动强度等对人群食物摄入量有不同程 度的影响 设计抽样方案显得异常重要 这里我们提出多层次划分法设计抽样 多层次划分 法是指根据不同的因素对调查人群进行划分 第一层为目标层 即所有要调查研究的对象 如全体国民 第二层则在第一层的基础上考虑某一因素比如年龄对第一层进行划分 分成几 类不同属性的人群 第三层在第二层基础上考虑另外一个因素比如地区对第二层的几类人 群进行划分 如此进行下去 直到满足一定的条件结束划分 划分过程如图1所示 人群食物摄入量模型要调查的对象是全国31个省 自治区和直辖市抽中样本住户的常住人口 当然多层次划分抽样调查的第一层为全国31个省 自治区和直辖市抽中样本住户的常住人口 接下来考虑第二层 将全国各地区按照经济发展情况 1 分成 811 812 816等六类地 区 其中 811包括北京 上海 天津 重庆 哈尔滨 沈阳 大连 济南 青岛 宁波 南京 广州 图1 调查抽样设计方案 西安 武汉 厦门 其余的城市为 812 813主要分布在长江三角洲 环渤海以及南部沿海农村 经济区 814主要分布在华北平原 四川盆地 东南丘陵以及豫皖鄂等长江中游农村经济区 815主要分布汾渭谷地 太行山 大别山农村经济区 其余地区为 816 根据1992年全国营养 调查的样本资料 2 进行抽样实验后 估计660000例样本就可以满足95 以上精度和准确 度要求 应答率按90 计 则本模型需要的调查样本量为740000个 第二层按地区经济发展水平划分好后 可进行第三层划分 第三层划分的因素可以选择 年龄 季节等 由于中国居民消费的食品种类比其他国家居民消费的食品种类复杂得多 包括 主食 肉类 蔬菜 水果 水 饮料 各种调味剂和经过加工的食品 细分将达数千种以上 在实际调 查过程中进行如此详细地分类 其调查工作量太大 而如果随意粗糙进行分类 则将影响调 查的精度 因此需要根据污染物分布模型的数据合理设计抽样调查中食物的分类办法 为 此 我们选用基于贡献率的食品分类法 以铅为例 铅的污染主要集中在蛋类 罐头蔬菜中 尤其是皮蛋 由于其特殊的加工工艺 已成为重金属污染食品的一个重要问题 根据污染物 分布模型 找出累计贡献率达到90 以上的食品作为监测对象 将最主要的污染食品皮蛋单 独分为一类进行监测 因此可将食品品种分为 粮食 蔬菜 水果 肉 鱼 蛋 皮蛋 奶粉 鲜 奶 豆类 小食品 罐头 饮料 茶叶 在进行抽样调查时 将重点监测的几种污染物 铅 镉 有 机磷 有机氯 每个污染物选取一中贡献率最大的食品单独最为一类 如铅对应的皮蛋 镉对 应的茶叶 有机氯对应的畜肉等 2 由上述方法确定了样本点后 可采用称重法 回顾法 食物频率法来作膳食调查 3 4 污染物分布模型 我国建立食品卫生安全保障体系的时间还不长 根据国际上的热点和我国的国情 据初 步估计 我国现阶段可能会集中力量对众多污染物中少数几种危害面广 后果严重的污染 物 如铅 镉 有机磷 有机氯等实行监控 其他污染物的监控工作则待时机成熟后再推广 由于已知化学污染物的浓度这一随机变量服从偏态分布 且抽样数据分布极不均匀 右 分位点以后的数据比较多 右分位点以前的数据只有总数据的2 如何充分利用已知数据 估计概率密度函数并用右分位点以后的数据作修正是概率密度函数估计的关键 可以先使 用分位点左侧的数据和右侧的部分数据进行概率密度函数估计 而用剩下的数据对估计结 果进行误差分析 评价估计的效果 这样复用数据可以达到提高精度的目的 由于这一实际问题我们能得到的有效数据极少 一般的假设检验统计方法恐怕难以胜 1514期唐苏文 等 食品卫生安全保障体系数学模型的研究 任 广义g h分布可以很好的估计偏态分布 需要的样本点少 估计精度高 Bootstrap估计 法是一种新的增广样本统计方法 这种方法可以由较少的原始样本数据再抽样 获得较多的 再生样本 可以大大提高估计精度 本文拟采用这两种方法估计污染物的概率分布函数 4 1 广义g h分布估计法 在实际应用中 经常使用正态分布来描述观测数据的概率分布规律 这要求观测数据的 分布是单峰对称的 但是本模型的污染物分布是单峰非对称的 与正态分布相比有两种变 化 一是分布有偏 二是尾部的概率有变化 这里我们采用Turkey等提出的一类分布模型 g h分布或广义g h分布 4 Z表示具有正态分布的随机变量 则随机变量 X X Z A B g h A B gZ 1 g e hZ2 2 1 服从g h分布 其中A B 0 g和h都是常数 这里约定 X Z A B 0 h li m g X Z A B g h A B Ze hZ2 2 2 假设g和h是Z的函数g Z 和h Z 即 X X Z A B A B e g Z Z 1 g Z e h Z Z2 3 称为广义g h分布模型 不难验证 g ZP g Z1 P g ZP 即g 是偶函数 故可取g Z 为Z 2 的多项 式 类似的 因为X0 5 XP X1 P X0 5 e g ZP ZP 所以h Z 也可取为Z 2 的多项式 最简 单的情况是 g ZP g0 g2Z 2 P h ZP h0 h2Z 2 P 4 其中g0 g2 h0 h2都是常数 这时有 g0 g2Z 2 P 1 ZP ln X1 P X0 5 X0 5 XP lnB h0 2 Z 4 P ln g0 g2Z 2 P X1 P X0 5 e g0ZP g2Z2 P 1 5 设有观测数据x1 x2 xn 欲用g h或广义g h分布模型来描述这组数据 将广义 g h分布建立观测数据的分布模型的方法和步骤参见文献 4 根据文献 5 所给六省样品中DON污染物浓度数据 下面利用广义g h分布进行概率 分布函数估计 因为污染物含量总是非负的 显然污染物含量越高则概率密度就越小 近似服从偏态分 布 可以用广义g h分布来估计其分布函数 按照上面的方法估计概率分布函数的过程如 表1所示 估计的概率分布函数如图2所示 由图2可以看出 玉米和小麦中DON含量的分 布函数是偏态分布的 用广义g h分布来估计其分布函数是合理的 表1 DON污染物浓度密度估计表 粮食种类 A g Z h Z B 玉米25 880 551Z2 0 06840 2194Z2 3 132622 934 小麦50 040 8694Z2 1 4690 4258Z2 2 04697 7424 25数 学 的 实 践 与 认 识38卷 图2 g h模拟方法获得污染物浓度概率分布图 4 2 Bootstrap估计法 1977年 美国Stanford大学统计学教授Eforn在前人研究成果的基础上提出了一种新 的增广样本统计方法 Bootstrap法 6 为解决小样本试验评估问题提供了新的方向 该 方法直接利用验前样本试验数据确定验前分布 并能充分利用子样本本身的信息 对于总体 分布不需要作假定 Bootstrap法的数学描述为 设随机样本X x1 x2 xn 来自未知的总体分布F R X F 为某个预先选定的随机变量 是X和F的函数 先要求根据观测样本X x1 x2 xn 估计R X F 的分布特征 例如设 F 是总体F的某个参数 Fn是观测样本X 的经验分布函数 Fn 是 的估计 记估计误差为 R X F Fn F Tn 6 现要求有观测样本X x1 x2 xn 估计R X F 的分布特征 实质就是一个再抽 样过程 传统的Bootstrap方法存在缺陷 随机样本的产生受原样本序列的影响较大 尤其是 大值与小值均受到限制 而本问题中需要求出99 999 的右分位点 对大值精度要求较高 因此传统方法不合适 采用指数函数拟合修正经验函数 7 能够使仿真生成的随机样本不受 原始样本序列范围的限制 得到了较好的评估结果 因此本文中也采用这个方法 其具体做 法如下 在n m个样本之前逐段构造线性经验分布函数 其中 n是样本总数 m是为尾部样本 数 其后的样本用指数分布拟合 并应使修正的经验分布均值与原始样本均值相等 通常m 取5以下的整数即可使右尾部拟合方差较小 修正的样本经验分布函数为 Fn x 0 x x 1 i n x x i n x i 1 x i x i x 0 0 x 0 其中 0 x 1e xdx 又设P x xlm in p1 P x xlmax p2 p1 p2可由尾部放大 法的前几步得到 又 p1 1 xlm in 0 x 1e xdx p2 1 xlmax 0 x 1e xdx 8 由式 8 可以计算出 和 值 问题转化为求xP使得 1 xP 0 x 1e xdx 99 999 查 分布表可求出99 999 右分位点xP 图4 尾部放大法的放大过程 6 若干理论问题的深入研究 以上我们分别建立了人群食物摄入量模型 污染物分布模型和风险评估模型进行食品 暴露评估 除此之外 还有一些理论问题需要研究 下面将研究几个典型问题 6 1 高精度地估计概率密度函数右分位点的方法 把每个人每天某种污染物摄入量看成是一个随机变量 我们关心的不仅是它的均值 更 关心的是它的99 999 右分位点 下面我们将重点讨论如何利用有限的抽样数据估计随机 变量的99 999 右分位点 本文在第五部分已经提出用尾部放大法可以较好地完成估计 当 然还有其他方法也可以胜任 这里再用一种方法进行尝试 本文在 4 2 节提出用Bootstrap法估计污染物的概率分布函数 实际上Bootstrap法的 再抽样思想完全可以估计99 999 右分位点 估计的步骤如下 5514期唐苏文 等 食品卫生安全保障体系数学模型的研究 图5 分布计算偏态分布分位点示意图 1 观测样本X x1 x2 xn 的值为有限的总体样本 由它可以构造样本的经验分 布函数Fn 2 根据经验累积分布函数Fn用仿真法随机抽取大样本 得到随机抽样集合X X 1 X 2 X N N 表示抽样次数 X k 为一个很大的数 被称为Bootstrap增广再生子样 表 示如下 X k x k 1 x k 2 x k n k 1 2 N 3 用符号 F 表示累积分布函数的99 999 右分位点 直接由Bootstrap增广再生子 样得到 F 的估计 F 6 2 分类标准不同时数据处理的方法 在本文的食品摄入量模型中 食品分类采用的是基于贡献率的食品分类法 当然也可以 采用其他分类法 当采用不同的食品分类法时 数据如何转化是必须要解决的 否则会造成 大量数据的浪费 下面我们从理论层面分析该问题 设有两种分类标准 一种为U u1 u2 uN T 另一种为V v1 v2 vM T 其中U是 随机变量u1 u2 uN构成的总体 V是由随机变量v1 v2 vM构成的总体 由于U和V都是由 基本的食品按照不同的组合形成不同的分类方法 因此一定存在一N M阶常数矩阵A 使得 U A V 那么 由于V的分布在可以统计得到 处理不同分类标准数据的问题就转化为已知V的概 率分布 求U的概率分布的问题 为此设ui ai1v1 ai2v2 ai MvM 其中v1 v2 vM均为偏态分布 若v1 v2 vM之间有相关性 直接求ui的分布是困难的 出于以下两方面考虑我们设v1 v2 vM之间 相互独立 其一是在某些场合这种假设是合理的 其二是V的各个分量经过适当的线性变换 可以转化为独立形式 有了v1 v2 vM之间相互独立的条件 推导ui的分布就容易多了 我们拟利用概率特 征函数推导ui的分布 设vj的密度函数为fvj v 定义其概率特征函数为fvj v 的Fourier 逆变换 由相关定理 8 可知 ui的概率密度函数为 fvi v 1 2 5 v1 ai1w 5v2 ai2w 5vM ai Mw e jw vdw 9 其中 5vj w f vj x e jw xdx 65数 学 的 实 践 与 认 识38卷 我们关心的并不是ui的概率密度函数 而是总体U的概率密度函数 这时 已知V的概 率密度函数要求其线性变换U的概率密度函数 需要用到Jacobin行列式 这就要求矩阵A 是可逆矩阵 定理2 设u g v A v 其中A为N N的可逆矩阵 则U的概率密度函数为 fU u fV A 1u det A 1u 10 6 3 利用地方调查数据估计全国水平的方法 当需要利用某些省 市的日常监测数据来估计全国的情况时也面临着两者的概率分布 函数可能并不相同的问题 如何利用部分省市的调查数据估计全国的分布是很有研究价值 的 本文提出线性转化法可以很好的解决该问题 设反映全国整体水平的随机变量为U 反映各个省市水平的随机变量为U1 U2 UN 其中U所反映的全国水平是U1 U2 UN所反映省市水平的总体 随机变量U1 U2 UN 的分布相似 反映部分省水平的随机变量U1 U2 UM的分布是已知的 问题的核心是用 U1 U2 UM的分布估计U的分布 既然U所反映的全国水平是U1 U2 UM所反映省市水平的总体 设U与U1 U2 UM的关系为U U1 U2 UM 下面本文拟用解析的办法确定U的分布 设随机变量函数U U1 U2 UM 的定义域为高维立方体 a1 b1 a2 b2 aM bM U1 U2 UM 的联合概率密度函数为fU u1 u2 uM 则 FU u P U u C arg U1 U2 UM u u 1 u2 uM fU u1 u2 uM M i 1 dui 11 这样就可以求出U的分布 6 4 衔接不配套抽样调查数据的方法 在本文的风险评估模型中 食品摄入量的数据要与污染物模型中的数据配套 当这两批 数据不配套时 就要设法处理这两批数据 否则会对数据造成浪费 设某个地区或全国某种食物的摄入量为随机变量X 某种污染物在该类食品中的含量 为随机变量Y 则该污染物在该地区的风险即摄入污染物的大小Z X Y也为随机变量 估 计出随机变量Z的分布特性就可以进行风险评估 设X的概率密度函数为fX x 定义域为 a1 b1 Y的概率密度函数为fY y 定义域为 a2 b2 下面在X与Y独立的情况下求Z的概率密度函数 由于X与Y是独立的 X与Y的 联合密度函数为 fX Y x y fX x fY y x y a1 b1 a2 b2 0else 12 则FZ z P Z z x y a1 b1 a2 b2 xy z fX Y x y d xdy 13 这样就可以求出Z的概率分布函数 求出概率分布函数后利用导数关系就可求出Z的概率 密度函数 7 结束语 本文主要针对食品安全保障体系中的食物摄入量模型 污染物分布模型和风险评估模 型做了以下研究 针对设计最优调查方案问题 提出一种基于多层次划分法的抽样方 7514期唐苏文 等 食品卫生安全保障体系数学模型的研究 案 此外 针对食品的分类问题 进一步提出一种基于贡献率的食品分类法 定义了概率 密度函数的估计误差 提出利用广义g h分布估计污染物的偏态分布的方法 进一步基于 Bootstrap法提出了一种利用少量数据再抽样估计偏态随机变量分布特性的方法 实例计算 表明上述两种方法在小样本的情况下有比较理想的效果 提出一种尾部放大法 巧妙地 将概率搜索和理论分析相结合 可以高精度地估计出污染物摄入量的右分位点 基于该方 法 本文进一步建立了点估计 单一分布评估和概率评估模型 基于以上三个模型 本文进一步深入地讨论了相关的几个理论问题的解决方案 收到较 好的效果 参考文献 1 中国财经网 http www 90531 com htm l 2555 htm l 2 王茂起等 中国2000年食品污染状况监测与分析 J 中国食品卫生杂志 2002 14 2 328 3 江国红等 膳食调查方法 J 职业与健康 2002 18 3 66267 4 侯紫燕等 一种适用于拟合单峰分布数据的模型 J 郑州大学学报 2001 33 2 26229 5 李晓云等 2005年中国居民膳食DON污染调查及暴露评估 J 长治医学院学报 2007 21 2 1012103 6 黄玮等 极小子样试验的虚拟增广样本评估方法 J 西北工业大学学报 2005 23 3 3842387 7 肖刚等 系统可靠性分析中的蒙特卡罗方法 M 科学出版社 2003 8 陈明 信息与通信工程中的随机过程 M 科学出版社 2005 6 9 盛骤等 概率论与数理统计 M 高等教育出版社 2001 12 10 朱道元 多元统计分析与软件SAS M 东南大学出版社 1999 8 Research on MathematicalM odel of Food Health and Security System TAN G Su2w en MA Shou2gui ZHAN G Hao N ationalM obile Communication Research L aboratory of Southeast U niversity N anjing 210096 China Abstract The paper mainly conducts the follow ing research on food intake model pollutant distribution model and risk asses