《∏Ukeλ≤∏Uk}~Poissonλ》由会员分享，可在线阅读，更多相关《∏Ukeλ≤∏Uk}~Poissonλ（9页珍藏版）》请在金锄头文库上搜索。
1、MCMC 1 随机数的采样 1. 随机数 - 随机变量的样本 (1). U0,1随机数 0,1均匀随机数 记为 U0,1随机数 . 电脑上常用同余 法得到伪随 机数列 , 是能通 过统计中的独立性检验及均匀性检验的周期列 , 但并非真正 的独立同分布列 . 典例如 : ）周期约为14441044121041(2,1,0),2(mod =+=+ nnnnnyxyyyyy (2). F-随机数 若分布函数 F(x)严增 , 服从 U0,1, 则 ( )1F 的分布函数为 F(x). (3). 离 散随机 数 若 . 把 0,1分为 n 个不交区 间 , 使第 个区 间 的长度为 . 作 U0,1随
2、机数 U , 则 就是一个nnppxxLL11=niJiIxi1iiJip U)( 随机数 . (4) 正态随机数 121, L为 i.i.d. U0,1 随机数 , 则 )1,0(6121N+= L, 若 为 N(0,1)-随机数 , 则 + 为 随机数 . ),(2N(5) Poisson 随机数 独立的 U0,1随机数 U , 则 . L,21U+=dd:)( Rxx= )K(,0i 0)( , x 多数的 Monte Carlo 计 算 问题, 常 归 结为计算甚高维积分或庞大和数的数值计算. 而 对于甚高维积分, 或庞大 求和, 用 静 态 Monte Carlo 方法也难 以实现,
3、 因 为其计算量海大. 更适用的也是用构造 Markov 链(即 MCMC)来模拟计算. Monte Carlo 方法的长处是计算的复杂度不依赖于计算空间的维数 . 所以 对计算甚高维积分有明显的优越性 . 2. Markov 链 Monte Carlo (MCMC) 0. MCMC 可以用于 null 实现对高维分布或高维格点分布 采样, 得到 随机数; null 实现重要性采样; null 实现高维积分或项数极多的求和的数值计算; null 模拟估计最可几轨道. 例如, 如 果 模拟了 100 条轨道, 那 么就能以大概率推断, 最可几轨道就在这些轨道的领域; 在 统计量的分布未知时, 用
4、 模拟得到置信限或置信概率 null 求复杂样本空间上函数的极值(模拟退火). 设 = 是一个只取有限个值 的概率分布，(但状态 数 甚大, 例 如是一个甚高维的概率密度函数的格点近似的归一化（即除 以格点值的和）), 或者是一个超高维的 密度 1 ， 假定 i (或 x )。 我们要构造一个以 为极限分布的 Markov 链，其目的 为：得到 -随 机数；得到关于 的各种统计平均值的数值估计。 1 通过条件分布对分布 作随机采样的 Gibbs 方法 设高维总体的分布为 . Gibbs 样本生 成法, 是一 种有效可行的高维总体或复杂总体的采样法 它的想法在于构造一个互通的正常返的马氏链 ,
5、使 得这个马氏链以0: nXn为极限分布 由于当 N 充分大时, 的分 布渐近于分布NX , 就可以近似地作为分布NX 的总体的样本. (1) Gibbs 样本转移的生成 Gibbs 样本 生成法的要义, 也是把高维总体的采样, 化为一系列一维分布的采样. 在具体运作时, 对于较大的 的一个 维分布密度d d ),(1 dxxL (可能 包含一个未 有显式表示的常数因 子) , 记 ),(dxxL,1x= , ),dy,(1yyL= . 再记在第 1 至 k-1 个分量固定为且在第 k +1 至第 d 个分量 固定为 的条 件下, 第 k 个分量 取 的条件密度为 11,kyyLdkxx ,1
6、L+ ky|(ky11,kyyL, ),1 dkxxL+kdkkdkkdyxxyyxxyy+=),(),(1111LLLL),(),(11111dkkdkkxxyyCxxyyLLLL+=我们定义如下的转移概率 =dkxyp1|(ky11,kyyL, . ),1 dkxxL+则易验 证 ( 是一个转 移矩 阵, 且)xyp )(x 是它 的不变 密度, 即若取 )(x 为以 为转移矩阵的 Markov 链 的初始分布密度, 那么对于任意 n, 有分布密度)(xyp)(xnXnX . (2) Gibbs 样本生成法的具体作法 从马氏链 的样本 , 如何 求得下一步转移到 的样 本 呢 ? 由 构造
7、, 是 个一维条件分布的乘积, 再 注意到 Von Neumann 取 舍原则中的 可以与密度差一个常 数倍, 因此 在 , 固定的条件下, 要 获得服从以 为密度nXd),(1 dxxxL=11,kyyL1+nX y(xhxyp|k)nkxx ,1L+(y11,k, yLy , 的样本, 我 们可以令),1 nkxL+x ),)(dkkxyyhLL,1,(1ykx+= , 并利用 Von Neumann 取 舍 原则采样, 得 到样本 . 对于 从 1 到 合起来, 就得到 的一个样本 ky k1+nXd= (yy ,1)my,L现在我们 任取 , 按上面方法可以得到 的一个样本 ,以 此进
8、行 , )0(0yX =1X )1(y归 纳 地 我们得 到 的样本 当 充分大 时, 由于 马 氏 链 的分布近似于NXX ,2L)(,),2( NyyLNnX)(x . 于是 我们就可以认为 就是近似地服从分布)(Ny )(x 的一个样本 ,KS =)(), ijijppijij=iipijp)(ijpX PnnT1TPX,LiN ,NXLi),ppijijij1min(,1ppspijijijijij+ijpjipj )ij2 Hastings-Metropolis 采样法 (1) Metropolis sampler 设状态集合为 ,1L. Metropolis 构造了以 ,1min(
9、 , ( ij1 ) 为转移的时齐 Markov 链 , 其中 是一个待选的对称的， 互通的非周期转移矩阵 , 使用它是为了减少状态间的连接 , 以加快 Markov 链 的分布向平衡分布P =n收敛的速度 . Metropolis 提出的这种采样法, 称为 Metropolis 采样法. 在上面的条件下可以证明： (1) , 其 中 1 表示分量全为 1 的 列 向量. 因此 , 在时间充 分长以后, 可以用 Metropolis 构造的 Markov 链所处的 状态 , 作为按 分布 采集的样本. 所以它也给出了在计算机上模拟 - 随机数的一个方法. (2) 我 们 有 遍历性定理 (这时
10、并 不 需要待 选矩 阵 的非 周期 性假定 ): 无 论从什 么初始条件出发， 由转移矩 阵 P得到的 Markov 的样本 ， 以概率为 1 地成立下面 的近似关系： ,NX1当 大时， 中取到状态 的频率1X 。 (3)使用 非对称的 可能 加快收敛速度 . 这 时的典型转移有 P(,jppijij= i . (2). Hastings 推广了的 Metropolis 采样 )( ipij= , 其中 为非对 称的转移矩 阵 , 满足 : 与 要么同时为 零 , 要么同时为正 , 而 是一P (s个非负的对称矩阵 , 满足 )(,11ijppsijijijij+. 3 MCMC 的用法
11、(1). 的平均量的模拟计算 对以 为极限 的 Markov 链 , 我们可 以用它的一 段样本 , 给 出平均量 的估计NXX ,1L=KiF1iifNfFNnXn=1（ 虽然很大， 但是一般比N K 要小得多） , 它是无 偏估计. 且 的方差 F )1()nO=(FVar 。 （2）实现 Importance Sampling 作为积分 的近似计算， 我们只须考 察定义在庞大格点集 G 上的函数 格点集 G 。 xxf ),(记 =Gyyfxfxx|)(|)(|)(:)( 。 构造以 为极限分布的 Markov 链 (注意在 Metropolis Sampler 的 转移概率的表示式中只
12、出现比值ij, 所以并不涉及计算和数 , 这才能有效地构造所要的 Markov 链 ) 取其 一条样本 。 当 大时， 的分布就近似于Gyxf )(|NXNXX ,1LN 。 因此 就近似地为 的 Importance Sampling. NX )(xf(3)甚高维积分或庞大求和的模拟计算 我们沿用 （ 2） 的记号与 结论。 当 大时， 中取到状态N ,1 NXXLx 的频率 ， 因而当 大时，)(xNGy NXyf的频率,|)(|1xxf中取到状态|)(X,L， 且这个近 似式与 x 无关。 注 如果在 具体运行中, 计算结果依赖于 x, 则可以 多做几条样本, 巴对应的 结果作平均. 3
13、. 模拟退火 (simulated annealing) . 模拟退火方法的基本想法 求目标函数 的最小值时, 通常的搜索 过程, 例 如 梯度方法, 常会停止在局部极小值点. 模拟退火算法时为了克服这个困难而设计的. , 其思想是: 引进人工干扰, 使得当 算法达到局部极小值时，能有一个小概率“逃出”局部极小值的陷井。 f假定自变量只取限个值. 而加入随机干扰的算法, 就是把原来的计算过程变成随机过程变成一个马氏链， 即用 MCMC 构 造一个马氏链 , 使其 步转移矩阵的分量 0:)(nXnn)()(npij )()()(*1fifjfjfiee=)()(1ifjjfee=, 其中 是目标
14、函数， ，而是一个随机性水平的参数, 称 为 倒温度. 这 个分布称为平衡分布, 它 有一个重要的性质: 当f )(min*iffi= 时, 其极限 分布集中在目标函数 取最小的 集合 上, 即在 上. 于是人 们可 以利用 这个 马氏链 的发 展, 到达其 平衡分布, 然后令 , 以达到 取最小值的点集上. 即 f)(:*fifi =f nlimlim *)()( fXfn=而对固定的, 马氏链 的采样, 则 可以通过 Gibbs 采样法 , Metropolis 采 样法等实 现. 的过 程就是 降温 的过程 , 所 以把这 种搜 索整体 最小 的方法 , 叫 做模拟退火方法. 0:)(n
15、Xn2 用非时齐的马氏链作模拟退火 为 了 加 快 计算速度, 在理论上发 现可采取非 时齐的马氏 链, 以把马 氏链的时间 发展(t )与 降温 ( )这两步 , 归 并成一 步, 即让此 马氏 链的第 步的 转移矩 阵随nn而改变. 当 时间 趋于无穷时, 只要nn 以适当缓慢的速度趋于 , 则对无论 什么初分布0 , 此 非 时 齐 的 Markov 链 n , 在 时 刻 的分布 , 其 中 为集合上的离散均匀分布. 从而有 n()(n)*1)(min=)(: ifi )( =Pnff . 如果温度n1=T 减小过快 , 则在 直观上就相 当于在每个 温度段都还 未达到平衡 , 于是退火过程仍 旧会陷入
