生物数学模型参数估计.pdf

生物数学模型参数估计生物数学模型参数估计 陕西师范大学陕西师范大学 数学与信息科学学院数学与信息科学学院 唐三一唐三一 2015年年7月月16日日 2014年广东省年广东省Dengue疫情数据疫情数据 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 13579111315171921232527293133353739414345474951535557 广州新增病例 广州累计病例 广东新增病例 广州累计病例 新增住院 现住院病例 累计住院病例 累计治愈出院 疫情数据与实际问题的解决疫情数据与实际问题的解决 基本再生数是多大 基本再生数是多大 Basic reproduction number 疫情的峰值及其达峰时间是什么 疫情的峰值及其达峰时间是什么 Peak size and peak time 到底有多少人会染病 到底有多少人会染病 Final Size 如何评估控制措施的有效性 如何评估控制措施的有效性 什么样的控制措施是最有效的 什么样的控制措施是最有效的 疫情数据与实际问题的解决疫情数据与实际问题的解决 Time Number of pests pest population time to treat time to treat EIL ET Time1 Time2 Economic Injury Level EIL lowest population density that will cause economic damage Economic Threshold ET population density at which control measures should be invoked to prevent an increasing pest population from reaching the economic injury level 综合害虫治理与最优经济阈值估计综合害虫治理与最优经济阈值估计 小区小区 1 2 3 4 5 6 7 阈值 头 百株 1 2 5 10 20 40 60 0 10 20 30 40 50 60 70 80 90 南1 南2 南3 南4 南5 南6 南7 2012年南边七个小区的盲蝽种群动态图年南边七个小区的盲蝽种群动态图 0 20 40 0 20 40 0 20 40 0 20 40 0 20 40 0 20 40 60 80 May 28 Jun 2 Jun 7 Jun 12 Jun 17 Jun 22 Jun 27 Jul 2 Jul 7 Jul 12 Jul 17 Jul 22 Jul 27 Aug 1 Aug 6 Aug 11 Aug 16 Aug 21 Aug 26 Aug 31 Sep 5 0 20 40 60 80 1 2012 Middle 2012North 2012South Action threshold a b 2 c 5 d 10 e 20 g f 4060 在防治虫量为在防治虫量为 1 a 2 b 5 c 10 d 20 e 40 f 60 g 头头 百株下百株下 小区内盲蝽的种小区内盲蝽的种 群动态 注 黑群动态 注 黑 虚线为防治虫虚线为防治虫 量 超过了该虫量 超过了该虫 量即打药防虫 量即打药防虫 0102030405060 60 70 80 90 100 110 0102030405060 0 2 4 6 8 10 12 024681012 60 70 80 90 100 110 B Yield kg Treatment No of insects A Number of spray Treatment No of insects C 2012 Middle 2013 south 2013 Middle 2013 North Yield kg Number of Spray 设定的防控虫量 棉田小区产量及打药次数设定的防控虫量 棉田小区产量及打药次数 之间的关系之间的关系 012345 0 10 20 30 40 50 60 70 80 90 100 110 120 0510036904812160102005 10 15 20 25 30 7月12日 棉枝顶端受害花蕾总数 p 0 584 7月17日 p 0 305 p 0 112 7月22日 小区实际调查的虫量 p 0 005 7月27日 p 0 000 8月2日 p 0 002 8月7日 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 100 204060800 10203040506070102030048121602468 10 12 p 0 0157 8月11日 棉枝顶端花蕾受害总数 p 0 011 8月16日 p 0 015 8月21日 实际调查的虫量 p 0 0025 8月26日 p 0 042 8月31日 p 0 000 9月5日 2012年每小区棉株顶端花蕾受害总数与小区内实际调查到的虫年每小区棉株顶端花蕾受害总数与小区内实际调查到的虫 量之间的关系 量之间的关系 21个小区图中共个小区图中共21个点 个点 常用数学软件常用数学软件 常用参数估计方法及其简单程序实现常用参数估计方法及其简单程序实现 绪论绪论 Bayes公式和线性回归模型公式和线性回归模型 最小二乘法最小二乘法 极大似然估计极大似然估计 Gibbs抽样技术抽样技术 Metropolis Hastings算法算法 MCMC 单种群生物动力系统单种群生物动力系统 科学科学出版社 出版社 20082008 552552页 页 6767万万字字 唐三一 肖燕妮唐三一 肖燕妮 生物数学进阶生物数学进阶 西安交通大学出版社西安交通大学出版社 肖燕妮 周义仓 唐三一肖燕妮 周义仓 唐三一 参考书目参考书目 参考书目参考书目 绪论绪论 参数估计参数估计 利用已知的实验观测数据并结合数学模型来估计模型利用已知的实验观测数据并结合数学模型来估计模型 中的未知参数中的未知参数 是统计推断的基本问题之一是统计推断的基本问题之一 参数估计的两种方法参数估计的两种方法 点估计点估计 是依据样本估计总体分布中所含的未知参数或未知参是依据样本估计总体分布中所含的未知参数或未知参 数的函数数的函数 其常见方法有其常见方法有 最小二乘法 极大似然估计法 贝叶斯估计法最小二乘法 极大似然估计法 贝叶斯估计法 区间估计区间估计 是依据抽取的样本是依据抽取的样本 根据一定的正确度与精确度根据一定的正确度与精确度 的要求的要求 构造出适当的区间构造出适当的区间 作为总体分布的未知参数或参作为总体分布的未知参数或参 数的函数的真值所在范围的估计数的函数的真值所在范围的估计 其常见方法有其常见方法有 利用利用已知的抽样分布 利用区间估计与假设检验的联系已知的抽样分布 利用区间估计与假设检验的联系 利用利用大样本理论大样本理论 1 Bayes公式和线性回归模型公式和线性回归模型 1 1 Bayes公式和似然函数公式和似然函数 1 2 一元线性回归模型一元线性回归模型 1 3 多元线性回归模型多元线性回归模型 1 4 Bayes线性模型线性模型 1 1 Bayes公式和公式和似然函数似然函数 对某个物种具有 年的统计观测数据 记为 1 2 我 们希望利用这组统计数据估计该种群的内禀增长率和环境容纳 量等 个未知参数 记作 1 2 即在给定的条件 下如 何确定 的问题 也就是在 发生的条件下 中各个参数值可能性 大小的问题 记作 在参数估计中 通常知道未知参数服从某 个分布比如正态分布或Gamma分布等 它是在不知道事件 是否 发生的情况下对未知参数 可能性大小的认识 记作 这样利 用Bayes公式得到 1 1 p Yp pY p Y 1 1 Bayes公式和公式和似然函数似然函数 在公式 1 1 中 称作先验概率 prior distribution 称为后 验概率 posterior distribution 边缘密度函数 公式 1 1 综合 了先验信息与实验提供的新信息 新的观测统计数据 并反映了 先验分布向后验分布的转化 通常情况下参数 是参数空间上的一个随机变量 这样样本 的分布族应理解为条件分布族 即不论是连续随机变量还是离散 随机变量 都可表示为依赖参数 的条件密度函数 其中 1 2 为依赖 的统计观测值 根据公式 1 1 可以看 出 统计观测值 的产生要分成两步 第一步第一步 设想从先验分布 产生一个样本参数 这一步是不 知道的 只能从人们对参数分布 的了解和经验判断得到 1 1 Bayes公式和公式和似然函数似然函数 第二步第二步 根据取样的参数值 从总体分布p 产生一个样 本 这一步是具体的 人们能够推断的 可理解为对于给定的 参数 求解模型 模型的解在给定时刻上的值即为样本 根据联合密度概率公式 此样本 发生的概率为 公式 1 2 综合了总体信息和样本信息 常称为似然函数 下面给出一个关于似然函数的例子 1 1 2 n i i p Yp y 1 3 L Yp Y 1 1 Bayes公式和似然函数公式和似然函数 例例1 1 考虑对一个物种密度或数量 的 次统计 假设每一次统 计误差 服从正态分布 0 2 其中统计精度 2是未知的参数 此时 1 2 样本分布的似然函数 样本参数 是设想出来并由先验分布 随机取样的 因 此不仅要考虑 还需要结合先验信息 即对 的一切可能取样加 以考虑 此时样本 和参数 的联合分布为 22 2 1 1 exp 1 4 2 n n i i pYL Ypyp 1 5 p Yp 1 1 Bayes公式和似然函数公式和似然函数 在给定先验分布和样本分布以后 可以根据 Bayes公式 1 1 计算 后验分布 由于对于所有可能的 值 边缘密度函数 与 没有关系 即相对于 来说 是一个常数 因此在计算后 验分布 时 通常忽略因子 而把 Bayes 公式 1 1 改写为 如下等价的形式 当样本分布为似然函数时 1 6 pYp Yp 1 7 pYL Yp 后验分布似然函数 先验分布 1 1 Bayes公式和似然函数公式和似然函数 其中符号 表示两边仅相差一个不依赖参数 的常数因子 比如 根据例 1 1 给出的似然函数 后验分布可以等价地改写为 22 2 1 1 exp 2 n n i i pYL Ypyp 1 2 一一元线性回归模型元线性回归模型 一元线性回归模型在生物数学模型中有非常广泛的应用 生 物统计学中的相关与回归分析既有联系也有区别 相关是研究变 量间是否有关系以及关系的程度如何 当通过相关分析 确定了 变量间存在着相关关系后 然后寻求一种函数关系式来构建它们 之间的关系 这就是回归分析 所构建的函数关系式就称为回归 方程 最简单的回归方程是一元线性回归方程 它是用于分析一个 自变量 与一个因变量 之间线性关系的数学模型 一般形式为 其中 0 1称作回归参数 它们通常是未知的 需要结合实际数据 并通过一定的统计方法进行估计 01 1 8 YX 1 2 一元线性回归模型一元线性回归模型 在实际应用中 当我们获得具有相关关系的两个变量 和 的多组观测值后 如果能通过适当的方法确定两个未知参数 0 1 这样就确定了正确的回归模型 1 8 对未知参数 0 1进行 估计时常用的方法就是最小二乘法和 Bayes方法 上面提及到一元线性回归模型在生物数学模型中具有非常 广泛的应用 这是因为实际应用中很多生物模型的参数估计问题 都可以直接或间接地利用一元回归模型的参数估计方法 如下面 给出的几个例子 例例1 2 设一个种群的数量 符合Malthus增长规律 1 9 dN t rN t dt 1 2 一元线性回归模型一元线性回归模型 初始为 0 0的解为 0 由此可以看出种群在任何 时刻的数量完全由参数初始 0和内禀增长率 确定 关于这两个 参数的估计可以转化为如下的线性回归方程 这样估计初始 0和内禀增长率 就转化为确定回归方程 1 10 的 回归系数 0 1 例例1 3 当考虑 Logistic模型的最大持续收获时 其在系统平衡态处 的持续产量为 01001 ln ln 1 10 YtYNNr 其中 1 1 11 qE qENqKE r 产量 1 2 一元线性回归模型一元线性回归模型 如果记 则方程 1 11 可以转化为如下形式的回归方程 这样估计参数 和 就转化为确定方程 1 12 的回归系数 0 1 例例1 4 对于离散的Ricker模型 先变形然后在两边取对数得 因此 如果记 则得到相同的线性回归方程 2 01 q KYield qKY rE 01 1 12 YE 1 exp1 1 13 t tt N NNr K 1 ln t t t Nr rN NK 1 01 ln t t Nr Yr NK 01 t YN 1 2 一元线性回归模型一元线性回归模型 实际上 在进行数据统计或实验的过程中 不可避免地要受到统 计误差 实验误差和实验条件等因素的随机干扰 当在模型 1 8 中考虑随机因素对因变量 的影响时 回归方程 1 8 变为如下的 一元回归模型 其中 为不可观察的随机变量 也称为随机误差或残差 一般来说 服从均值为零 方差为 2的正态分布 即 0 2 由于随 机变量的期望值为零 因此有 此时 可以看成一个随机变量 方差 2的大小反映了该随机变量 的分散程度 方差越大随机变量取值越分散 即波动越大 方差越 小变量取值越集中 即波动越小 01 1 14 YX 01 E YE Y XX 1 2 一元线性回归模型一元线性回归模型 一元回归模型在拟合具有噪音 随机干扰 的观测数据中的应用是 非常广泛的 下面利用实例加以说明 例例1 5 对白鼠从出生后第3天起 每隔3天称一次体重 一直称到第 18天 数据如下表所示 试计算日龄 与体重 的回归关系 表表1 1 白鼠日龄和体重的关系白鼠日龄和体重的关系 若把表1 1中的每对数据 对应于横轴为 和纵轴为 的坐标 系中的一个点 并把所有的点都在坐标系中标出 就得到图1 1 序号序号 1 2 3 4 5 6 日龄 3 6 9 12 15 18 体重 6 12 17 22 25 29 1 2 一元线性回归模型一元线性回归模型 图图1 1 白白鼠日龄和体重的关系图鼠日龄和体重的关系图 其中圆圈表示其中圆圈表示表表1 1中中的实际的实际数据数据 实点为实点为最最 小二小二乘法拟合得到的数值乘法拟合得到的数值 1 2 一元线性回归模型一元线性回归模型 由于任何一条直线的方程都可以写成 的形式 基中 表 示直线的斜率 表示直线与 轴的截距 如果 给定了 那么这条 直线就给定了 但是现在我们所考虑的是具有回归关系的量 是 随机变量 因此可以假设随机变量 与自变量 有关系 其中 是两个参数 是随机误差 要求 0 2 对于样本观测值 就有 我们的任务就是设法确定参数 使得直线 能较好地反映点 之间的变化关系 这条直线我们将称为回 归直线 yabxe 2 1 2 0 iiii yabxeineN 其中 yabx 1 3 多元多元线性回归模型线性回归模型 影响生物个体的因素是多样多样的 比如身高 体重等 当因变 量受多个自变量的影响时 就需要利用多元回归模型 多元线性 回归模型与一元线性回归模型形式上基本类似 只不过自变量由 一个增加到两个以上 并且因变量 与 个自变量 1 2 之 间存在线性关系 即 其中 0 1 为 1个未知参数 为随机误差项 公 式 1 15 称为多元线性回归模型 由于随机误差项的期望值为零 模型 1 15 两边求期望便能得到相应的多元回归方程 当观测的数据为 组 即 yi 1 2 1 2 时 多元线性回归模型 1 15 变为如下方程组 01122 1 15 kk YXXX 1 3 多元线性回归模型多元线性回归模型 其矩阵形式为 与一元线性回归模型一样 对于多元线性回归模型 我们仍旧希 望根据观测样本估计模型中的各个参数 即 0 1 有关的参数估计方法包括最小二乘法 极大似然估计和 Bayes 方法 所以在介绍上述三种参数估计方法之前 首先介绍下面的 Bayes线性模型 它是Bayes统计推断的基础 01 122 1 16 iiikkii yxxx 11121100 21222211 12 1 1 1 k k nnnknkk yxxx yxxx yxxx 1 4 Bayes线性模型线性模型 由模型 1 16 如果记 0 1 1 2 则多元回归模型可以 描述为 或 其中 1 2 为因变量样本观测值向量 1 2 为随机误差项向量 0 1 2 1 17 k ijjii j yxin 1 18 YB 11211 12222 12 1 1 1 k k nnkn xxx xxx B xxx 1 4 Bayes线性模型线性模型 为自变量观测值矩阵或设计矩阵 如何根据已有的观测值和 Bayes方法估计参数向量 根据Bayes公式和关系 1 7 需要按 以下三个步骤进行 为所有未知参数选取先验分布 即确定 其中 2 给出给定参数下观测数据的似然函数 或 利用 Bayes 公式 1 1 确定后验分布 2 最小二乘法最小二乘法 最小二乘法是根据观测数据以误差的平方和最小为准则 估计线性或非线性模型中未知参数的一种基本参数估计方法 它首先在1794年被德国数学家高斯解决行星轨道预测问题时提 出 其基本思路是选择估计量使模型输出与实测数据之差的平方 和达到最小 线性最小二乘法是应用最广泛的参数估计方法 它在 理论研究和数据拟合中都具有重要的作用 比如可以通过最小化 误差的平方和来寻找数据的最佳函数拟合 2 1 线性模型的最小二乘法线性模型的最小二乘法 2 2 非非线性模型线性模型的的最小二乘法最小二乘法 2 1 线性模型线性模型的最小二乘法的最小二乘法 对于例1 5中给出的白鼠日龄和体重的数量关系 我们希望 知道什么样的直线能够较好地反映点 之间的变化关系 进而可以利用这一关系来预测后期白鼠体重随时间的变化 不妨假设点 满足直线关系 那 将表示观测值与回归线的误差 这样 对于所观测的 个点 来说 如果我们给 以适当的值 使得离差 的平方和 达到最小 那么就认为由 所决定的直线是最好的 利用这种办 法来确定位置参数的方法称为最小二乘法 那怎样选择参数 使得离差的平方和达到最小呢 这实际上是一个求最小值的问题 利用函数求极值的方法我们分下面三步来实现这一目标 2 2 11 nn ii iii Qeyabx 2 1 线性模型的最小二乘法线性模型的最小二乘法 将 对 求偏导数 令它们各自等于零有 整理上面两式 得到关于 的线性方程组 关于 求解上面的方程组 得到 的估计值分别为 1 1 2 0 2 0 n ii i n iii i Q yabx a Q yabx x b 11 2 111 nn ii ii nnn iii iiii naxby xaxbx y 111 1 2 2 2 1 11 1 1 nnn n iiii ii iii i n nn i i i iii aybx x yxyxxyy n b xx xx n 2 1 线性模型的最小二乘法线性模型的最小二乘法 现在就可以利用表1 1中的队列数据 并结合 的计算公 式得到参数 的估计 当统计数据很多 上面的计算还是相对比 较复杂的 此时我们可以借助数学软件中常用的曲线拟合 多项 式拟合命令得到相应的参数估计 对于形如一元线性回归或多元 线性回归模型的最小二乘法估计问题 数学软件Matlab中有一个 命令 可以直接得到参数 的估计 比如当 分别为表表1 1中的时间 向量和白鼠体重 利用命令 1 得到两个估计值为 由此得到拟合数据如图1 1中的实点 Par polyfit x y 1 2 6 1 5143 ab 图图1 1 白鼠日龄和体重的关系图白鼠日龄和体重的关系图 其中圆圈表示其中圆圈表示表表1 1中的实际数据中的实际数据 实点为最实点为最 小二乘法拟合得到的数值小二乘法拟合得到的数值 拟合结果拟合结果 Polyfittsy m 2 2 非非线性模型线性模型的最小二乘法的最小二乘法 关于非线性模型的最小二乘法 希望通过下面的实例分 析介绍如何利用数值方法寻求参数估计和数据拟合问题 利用模型 来拟合英国传染病监测中心在1978年发布的有关流感病人的统 计数据 该数据统计了两个周内每天英格兰北部一所男孩寄宿学 校流感爆发和流行的情况 总共有763个学生 从发现一个染病者 开始统计染病学生的人数 总共统计了15天 每一天统计得到的染 病人数分别为 2 1 dS t S t I t dt dI t S t I tI t dt 1 3 7 25 72 222 282 256 233 189 123 70 25 11 4 2 2 非非线性模型的最小二乘法线性模型的最小二乘法 这样模型 2 1 中两个变量的初始值分别为 0 762 0 1 利用模型 2 1 中的 来拟合上述数据 关键是确定模型中 的两个未知参数 即参数向量 记 表示模型 2 1 在 给定参数 下对应的数值解 其中 2 对应于第二个方程的数 值解 并用向量 表示15个实际统计数据 不妨把第一天的时 间记为 0 0 则 0 1 2 14为对应于15个统计数据 的时间 因此为了寻求最佳的数据拟合 必须找到最佳的参数组 和 使得 达到最小 14 2 2 0 2 2 kk k LSf tID 2 2 非非线性模型的最小二乘法线性模型的最小二乘法 由于模型 2 1 没有解析解 2 的形式非常复杂 此时使得离 差的平方和 达到最小就是一个非线性的最小二乘法问题 可以 借助数学软件 中 和 等命令进行数值 求解 由于参数 和 具有明显的生物背景 所以它们的取值应该是 非负的 这是一种先验信息先验信息 所以在寻求这两个参数的估计时 我 们可以利用这一信息而对参数 和 加以约束 求解这一具有约束 条件的优化问题 我们利用命令 来实现 我们用我们用 来实现来实现 Parest SI mParest SI m 程序中的 最终输出参数 和 的估计值分别为 0 0022 0 4529 利用这两个参数重新运行微分方程得到统计数 据的最佳拟合曲线 结论显示模拟的感染者数量变化曲线与实际 数据非常吻合 如下图2 1所示 2 2 非非线性模型的最小二乘法线性模型的最小二乘法 图图2 1 模型模型 2 1 的模拟曲线与实际统计数据的比较的模拟曲线与实际统计数据的比较 圆圈表示统计数据圆圈表示统计数据 参数参数 为为 0 0022 0 4529 初始值为初始值为 0 762 0 1 3 非非线性模型线性模型的的极极大似然估计大似然估计 极大似然估计方法是求参数估计的另一种简单有效的方法 首先由德国数学家高斯在1821年提出 后来英国的统计学家Fisher 在1922年再次提出了这个思想 并给出了极大似然估计这一名称 以及这种方法具备的一些性质 极大似然估计的基本思想是 如果 已知某组参数能使样本出现的概率最大 就不会再去选择其它小 概率的样本 所以干脆就把这组参数作为估计的真实值 当然极大 似然估计只是一种粗略的数学期望 要知道它的误差大小还要做 区间估计 也就是说 选择参数 使已知数据 在某种意义最有可 能出现 某种意义是指使得似然函数 达到最大 这里 是数据 的概率分布函数 与与最小二乘法不同的是最小二乘法不同的是 极大似然法极大似然法需需 要已知要已知这个概率分布函数这个概率分布函数 这也使得极大似然估计方法在 应用中受到限制 3 非非线性模型线性模型的的极大似然估计极大似然估计 尽管如此 极大似然估计方法的应用还是非常广泛的 为了解 决这一问题 一般可假设 是正态分布函数 这时极大似然估 计与最小二乘法估计就一样了 一般来说 设总体 的分布密度为 其中 是需要估计 的未知参数向量 1 2 是总体的一个随机样本 它的联合 密度是 对于已知的 来说 函数 1 2 给出了随机变量 1 2 取值 1 2 的相对可能性 使得 1 2 取得极大值的样本观测值 1 2 就是 最 有可能出现 的数值 12 1 n ni i L y yyp y 3 非非线性模型线性模型的的极大似然估计极大似然估计 另外 如果在这个联合密度中 未知 而我们得到样本的观测 值是 1 2 这时这组观测值 极有可能 来自哪个密度函数 呢 这个密度函数的参数 应该是使得 1 2 作为 的 函数 取得极大的那个 值 于是我们就以 1 2 作为参数 的估计值 以这种思想得到的估计值称为极大似然估 计 为了记号简单我们记 如果 1 2 是使 达到极大的 值 那么称 1 2 为 的极大似然估计量 由于ln 是 的单调 函数 所以 和ln 在同一个 值达到它们的极大值 12 n LL y yy 3 非非线性模型线性模型的的极大似然估计极大似然估计 然而求ln 的极大值有时更容易 因此求 的极大似然估计时 只要求解下列的似然方程 即可 一般来说 当参数向量含有 个参数时的极大似然估计量 就是下面 个方程的解 下面给出求极大似然估计估计的一般步骤 12 ln 0 n dL y yy d 12 k 12 ln 0 1 2 k i L ik 3 非非线性模型线性模型的的极大似然估计极大似然估计 根据分布写出似然函数 对似然函数取对数 对ln 1 2 求偏导 建立似然方程 解出 1 2 得到参数 1 2 的极大似然估计 例例3 1 设 1 2 是取自正态总体 2 的随机样本 2是 未知参数 试求 2的极大似然估计 12 1 n ki i Lp y 12 ln k L 12 ln 0 k i L 3 非非线性模型线性模型的的极大似然估计极大似然估计 解 根据密度函数 得到似然函数 相应的似然方程为 由这两个方程分别求解 2得到极大似然估计量 它们刚好是样本的一阶原点矩和二阶原点中心矩 22 2 11 exp 22 p yy 22 2 1 11 exp 2 2 n i n i Ly 2 1 2 242 1 ln1 0 ln1 0 2 n i i n i i L y Ln y 222 11 11 nn iin ii yyyys nn 3 非非线性模型线性模型的的极大似然估计极大似然估计 在生物统计的许多问题中 极大似然估计是较为容易得到的 不 仅如此 在数学上还可以证明 在一定条件下 只要样本容量足够 大 极大似然估计量可以和参数的真值任意接近 而且在一定意义 上来说 极大似然估计是最好的估计 因此极大似然法因此极大似然法可以可以是是在在 理论上比较优良理论上比较优良 适用范围较广的一个估计方法适用范围较广的一个估计方法 下面向大家介绍如何一种利用统计观测数据和极大似然估 计法估计某种传染病比如A H1N1在一个地区的基本再生数的方 法 从发现病例的第一天开始 每天报告的新发病例数为 直到第 天结束 如下表3 1所示 这样报告的数据记为 01 3 1 T NNNN 3 非非线性模型线性模型的的极大似然估计极大似然估计 表表3 1 某地区日报告新发病某地区日报告新发病人数人数 记 为第 1 2 天续代时间续代时间的概率函数 其中续续代代 时间时间表示一个病人感染疾病出现症状到被这个病人传染的表示一个病人感染疾病出现症状到被这个病人传染的病人病人 也也出现症状的时间出现症状的时间 这里要求 是为了报告数据时间足够长 使得二代病例甚至是三代病例出现 以便正确估计基本再生数 在上述假设基础上 和 在2008年推导出基本再生数 0和概率函数p 1 2 由如下的似然函数确定 日期日期 3 4 5 6 7 8 9 10 11 12 人数人数 3 22 0 16 4 10 4 13 8 13 日期日期 13 14 15 16 17 18 19 20 21 人数人数 3 10 20 18 19 37 29 35 28 3 非非线性模型线性模型的的极大似然估计极大似然估计 其中 0 1 最大续代时间 对于A H1N1 6 和 是gamma函数 如果存在从外地的输入病例 记作 1 2 则似 然函数变为 有证据显示续代时间服从均值为 方差为 的gamma 分布 由于 对于某些疾病 续代时间已经研究得比较清楚 因此我们可以假设 均值 和方差 是已知的 然后利用极大似然估计就可以估计基本 再生数了 0 1 exp 3 2 1 t N T tt t t L Rp N N 1 exp 3 3 1 tt NY T tt c t tt L Rp N Y NY 3 非非线性模型线性模型的的极大似然估计极大似然估计 等式 3 2 两边取对数得到对数似然函数为 关于 0求导数并解似然方程得到似然估计值为 从公式 3 5 可以看出 如果我们能够给出 的计算公式 那 0就可 以利用表 3 1 中提供的数据 如下图3 1所示 估计参数 0了 由于 续代时间 服从均值为 方差为 的gamma分布 所以 可以由 0 1 ln lnln 3 4 T tttt t L Rp NNN 1 0min 11 3 5 T t t k tT jtj tj N R p N 1 1 3 6 j x j j pxedx 3 非非线性模型线性模型的的极大似然估计极大似然估计 计算 其中 是形状参数 是率参数 这样如果假设 4 1 6 0的极大似然估计值 0 1 663 图图3 1 某地某地区日报告新发病例区日报告新发病例数数 4 Gibbs抽样技术抽样技术 Bayes统计推断 的基本方法是将关于未知参数的先验信息与样本 信息综合 再根据贝叶斯定理 得出后验信息 然后根据后验信息 去推断未知参数 MCMC方法是最近发展起来的一种简单且行之 有效的 Bayes计算方法 主要由 Metropolis Hastings 算法和 Gibbs 抽样组成 基本思想是通过重复抽样通过重复抽样 建立建立一个平稳分布为一个平稳分布为所所求求 后验分布后验分布的的 Markov 链链 从而得到后验分布的样本从而得到后验分布的样本 然后然后基于这基于这 些些样本样本再作各种统计推断再作各种统计推断 比如比如估计参数的均值估计参数的均值 方差和自方差和自相相 关性等关性等 4 1 Gibbs抽样技术的具体抽样技术的具体步骤步骤 4 2 Gibbs抽样抽样方法方法的的应用举例应用举例 4 1 Gibbs抽样技术抽样技术的具体步骤的具体步骤 从公式 1 7 可以看出 联合后验分布密度函数通常由参数向量 所满足的概率分布决定 是一个条件分布 如果该条件后验分布是 已知的分布函数 那么就可以利用随机取样方法直接从该分布中 抽样待估参数的估计值 这种抽样方法称为 Gibbs 抽样方法 Gibbs 抽样方法可以描述为 设未知参数个数为 即未知参数向量 1 2 其中联 合后验分布为 p 1 2 通常来说这个联合分布的精确形 式是不能知道的 如果能够知道每一个参数 在给定其它参数条 件下的完全条件分布 111 iiil pY 4 1 Gibbs抽样技术抽样技术的具体步骤的具体步骤 就能利用Gibbs抽样方法获得联合分布的一组抽样值 过程如下 给定初始值向量 给定初始时刻 0 充分大的正整数 和 时开始每一步存储参数抽样 1 当 时结束循环 0 0 0 0 0 123 l 4 1 Gibbs抽样技术抽样技术的具体步骤的具体步骤 上述过程中的正整数 称为 期 它是为了避免未知参 数初始值和早期抽样值远离真实值而舍弃的前期抽样 这样当 充分大时得到一个关于未知参数的抽样 1 或 随机数列 它可以近似地看成联合分布 1 2 的一组 抽样值 利用此 随机数列就可以做相应的统计推断了 比如 可以求参数 的估计值 同时这一序列的标准差还可以用来做有关的假设检验 下面给出 一个具体的实例讨论如何根据似然函数 先验信息以及得到未 知参数的完全条件分布 然后根据 抽样方法的上述步骤实 现参数的估计 并给出相应的 程序 供大家参考 i 1 1 j M i j m Mm 4 2 Gibbs抽样抽样方法方法的的应用举例应用举例 为了说明 抽样技术的操作过程或 分析的基本思想 我 们考虑一组从均值为 标准差为 2的正态分布抽样得到的数据 1 2 这样第 次观测数据 的似然函数为 所以对于 个数据的联合似然函数为 2 2 2 2 1 exp 4 1 2 2 i i y Ly 2 2 2 2 1 2 2 2 2 1 1 exp 2 2 11 exp2 4 2 2 2 n i i n ii y Ly ynyn 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 具体目标是如何根据这组数据估计正态分布的两个参数即均值 和方差 为此适当选取这两个参数的先验分布 然后联合似然函数 并结合 Bayes公式 分别给出这两个未知参数的完全条件分布或 后验分布 从而利用Gibbs抽样方法估计这两个参数 为了实现上 述目标 我们分三种情况逐步加以考虑 方差已知 均值未知 假设方差 2是已知的 而均值 是未知的 即在给定方差 2 的前提 下如何得到均值 的完全条件分布 根据Bayes分析 我们还需知道 均值 的先验分布p 不妨假设其服从均值为 0方差为 0 2的正 态分布 即 N 0 0 2 所以 2 0 2 2 0 0 1 exp 4 3 2 2 p 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 其中参数 0和 0 2称为超参数 这样根据Bayes公式得到参数 的后 验分布为 其中 2 2 2 2 0 22 22 1 0 0 2 2 2 0 22 1 0 22 0 2222 00 2 0 2 0 111 exp2exp 22 22 1 exp2 22 exp 22 exp 2 n ii n ii pyLy p ynyn ynyn nyn 2 222 0 4 4 2 nyn 1 22 0 2222 00 1 nny 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 在 4 4 中最后等式中凑平方并忽略与 无关的常数因子 得到未知 参数 的完全条件后验分布服从均值为 和方差为 2的正态分布 即 此时我们注意到参数 的先验分布和完全条件后验分布具有完全 相同的分布类型正态分布 只是参数不同而已 因此如果一个参数 的先验分布与后验分布具有相同类型的分布 则称该先验分布是 其相应后验分布的共轭先验分布 均值已知 方差未知 如果假设正态分布的均值 是已知的 如何寻求方差 2的共轭先 验分布 2 2 2 exp 4 5 2 py 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 为此先介绍 Gamma 分布和逆Gamma 分布函数的有关性质 Gamma分布由两个参数 和 决定 其中 0为形状参数 0 为尺度参数 相应的密度函数为 或简记为随机变量 服从 Gamma 分布 即 设 1 通过概率运算得到 的密度函数为 1 exp 4 6 a a b p x a bxbx a Ga Xa b 1 exp 4 7 a a b p z a bzb z a 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 即逆 Gamma 分布 或简记为随机变量 服从逆 Gamma 分布 即 在正态分布均值 是已知的前提条件下 方差 2的共轭先验分布 通常取为逆 Gamma 分布 即 其中参数 和 称为超参数 实际上 当均值 为常数时 似然函数 4 2 与逆 Gamma 分布 4 8 的联合密度 即 2的后验分布为 IG Za b 2212 exp 4 8 a a b pb a 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 其中 从 4 9 不难看出 2的后验分布仍是参数为 的逆Gamma 分布 即IG 所以说 2 是当正态分布均值 已知时似然 函数 2 的共轭先验分布 222 2212 2 2 1 2 22212 2 1 12 2 11 exp exp 2 2 1 exp exp 2 exp 4 9 n a n a i i n na i i a pyL yp b yb a yb b 2 1 1 2 2 n i i aanbyb 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 均值未知 方差未知 当均值 和方差 2都未知时 可以结合前面两种情况得到后验分 布 这样只需把先验分布 2 看做为方差的条件分布 即 然后结合上面的两个步骤中的计算结果得到参数 和 2在给定观 测数据 条件下的后验分布满足 公式 4 11 给出了参数 和 2的完全条件分布 即参数 服从 参数 2服从 222 4 10 ppp 222 IG 4 11 pypy pyNa b 22 4 12 pyN 正态分布 2 IG Gamma 4 13 pya b 逆分布 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 利用 4 12 和 4 13 给出的两个条件分布进行抽样 就可获得一 个 Gibbs 随机数列 从而获得 和 2的抽样值 为了确定正态分布两个未知参数的完全条件分布 首先需要 确定这两个未知参数先验分布中的四个未知超参数 这也是一个 比较困难的问题 这里就不介绍了 下面只考虑最简单情形 即给 定四个超参数 0 0 0 2 1 2 1 为了说明Gibbs抽样 技术的有效性 我们通过模拟数据加以验证 为此从正态分布 1 42 中随机产生长度为 30 的样本观测值 即样本观测值的长 度 30 这样就可以根据上面介绍的步骤实施 Gibbs 抽样 下面介绍数值实现步骤和Matlab程序 任意设定未知参数的 初始取值为 0 0 2 下面的 Matlab 程序 MCMC Gibbs1 m 实 现了一维正态分布两个参数mu 和 sigma的Gibbs 抽样 得到了 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 相应的随机数列 图4 1给出了 9000 次burn in 期后的 Gibbs 随机 数列和相应的直方图 由此可以看出 Gibbs随机数列的均值与真 实值 1 4 非常接近 图图4 14 1 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 剩下的一个问题是我们希望知道先验分布是如何影响估计 值的精确性的 即是否得到的参数估计值依赖超参数的选取 为了说明这一点 我们在图4 2中给出了正态分布两个参数 先验分布和后验分布的数值实现 比较图4 2和图4 1可以看出 尽 管未知参数先验信息的均值与真实值相差很远 但最终得到的后 验分布却与Gibbs 随机数列确定的参数分布非常接近 这就说明 了当未知参数的共轭先验分布已知 就能得到该参数的完全条件 分布 从而实施Gibbs抽样技术得到未知参数较为精确的估计值 Gibbs抽样技术具有过程简单 容易数值实现等优点 4 2 Gibbs抽样方法抽样方法的的应用举例应用举例 图图4 24 2 5 Metropolis Hastings算法算法 Gibbs抽样技术在完全条件分布已知或易于求出的情况下是一种 行之有效的MCMC参数实现方法 在实际应用过程中 由于未知参 数比较多 要想得到所有参数的完全完全条件后验分布条件后验分布几乎是不可能 的 这种情形下Gibbs抽样方法就失效了 此时就需要应用MCMC 中的另外一种常用的方法即Metropolis Hastings 算法 简称MH 算法 MH 算法是比 Gibbs 抽样更一般化的MCMC方法 最初由 Metropolis 等在 1953 年提出 后来经过 Hastings 在 1970 年扩展 完善 现如今该方法在生物数学各领域的模型确定 参数估计和 数据分析等领域具有非常广泛的应用 这里我们就其最简单的情 形和基本思想向大家做个介绍 5 Metropolis Hastings算法算法 MH 算法的基本思想是从稳定的Markov链中抽样以模拟服 从一定概率分布的观察值 比如假设 MCMC 链的当前状态为 1 则Markov 链的下一状态的取值由如下两个步 骤产生 第一第一步步 从一个建议分布 中抽样侯选值 其中 第第二二步步 利用一个接受概率 来判断这个取样值 是否 可以作为 的下一个取样值 如果接受 则 1 否则拒绝 该次取样并且 1 即该链在这一步没有移动 所以说 MH 算法也称为接受 拒绝方法 1 1 111 tttt iin i 5 Metropolis Hastings算法算法 当未知参数个数比较多时 如果从建议分布 一次取样 一个参数并更新该参数 就会使得过程非常复杂或计算量非常大 所以我们有时可以从同一个建议分布中一次取样多个参数 称为 联合更新 这样在构造Markov 链时就有两种不同的 MH 取样方法 单参数更新和联合更新 单参数更新是指一次只更新一个参数 联合更新是指一次更新多个或所有待估计的未知参数 下面只介 绍联合更新方法的具体步骤 给定初始值向量 给定初始时刻 0 循环 0 0 0 0 0 123 n 5 Metropolis Hastings算法算法 从建议分布 中抽样 从均匀分布 0 1 中抽取随机数 如果 则 1 接受侯选值 否则 1 拒绝侯选值 1 经过一个 期后每 次循环后存储 1 当 充分大后结束循环 MCMC链收敛为止 5 Metropolis Hastings算法算法 MCMC链收敛的判别标准很多 可以根据观察 数学方法比如 Geweke收敛性判别等 从MH取样的具体步骤我们知道 实现MH 方法的关键是接受概率 的确定和建议分布q 的选取 下面分别加以简单说明 5 1 接受概率的确定接受概率的确定 5 2 建议分布的选取建议分布的选取 5 3 应用举例应用举例 5 1 接受概率的接受概率的确定确定 假设 是Markov链的当前取样值 从建议分布q 取样 如果 建议分布q 满足关系 则称q 具有可逆性 或者说选取的建议分布q 是一个正 确的Markov链的转移核 然而等式 5 1 在大多数情况下并不成立 则下面的情形之一一定发生 或 此时为了保证建议分布q 具有可逆性 需要在不等式 5 2 或 5 3 中大的一边乘以 5 3 qpYqpY 5 1 接受概率的确定接受概率的确定 此时q 可近似地视为转移核 在这种情形下接受取 样 的概率为 同样对于不等式 5 3 如果在右边乘以概率 则有等式 此时q 视为转移核 这样接受取样 的概率为1 综上所述 定 义接受取样 的接受概率为 5 4 qpYqpY 5 5 qpY qpY 5 6 qpYqpY min 1 5 7 qpY qpY 5 1 建议建议分布的分布的选取选取 利用相同的方法 可以得到一次更新一个参数或几个参数的接受 概率的计算公式 比如如果我们采用单参数更新方法一次更新一 个参数 那么更新参数被接受的接受概率为 min 1 5 8 iiii iii iiii qpY qpY ii i ii 5 2 接受概率的确定接受概率的确定 从接受概率的计算公式可以看出 更新参数的接受概率的大小取 决于建议分布的选取 建议建议分布分布的选取在一定程度上决定了参数 更新的快慢和速度 所以说选择合适的建议分布在实现MH算法 时就显得非常重要了 那在实际应用中如何选取建议分布呢 总 的原则就是使得更新参数的接受概率能够维持在一定水平上并 且计算起来比较方便 常见的选取方法有下面的几种 Metropolis取样取样或或对称分布对称分布 Metropolis在实现其提出的算法的同