随机数产⽣,MCMC⽅法,EM算法,bootstrap⽅法
统计计算⽂献报告
摘要：本⽂主要围绕随机数的产⽣算法、MCMC算法、EM算法、Bootstrap算法四种统计计算⽅法，介绍各⽅法的算法原理；并针对每个算法给出了实例。
关键词：随机数，MCMC⽅法，EM算法，bootstrap⽅法
1．随机数的产⽣
统计模拟的最基本问题是随机数⽣成，是按统计模型⽣成数据的基础。随机数⽣成可分成两类:[0,1]区间上均匀随机数⽣成和⾮均匀随机数⽣成，是随机数⽣成问题的两个基本研究领域[1]。
1.1 [0,1]区间上均匀随机数⽣成
[0,1]上均匀随机数⽣成是随机数⽣成之基础，⾮均匀随机数是由[0,1]上均匀随机数经相应运算⽽产⽣的。[0,1]上均匀随机数的质量决定了⾮均匀随机数的质量。如何快速地⽣成[0,1]上⾼质量均匀随机数是普遍关注的问题。⽣成⽅法分为两类：物理⽅法和数学⽅法[2]。
1.1.1 物理⽅法
基本原理为：利⽤某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产⽣随机数。这些特殊设备称为随机数发⽣器。把具有随机性质的物理过程变换为随机数，使⽤物理随机随机数发⽣器，在计算机上可以得到真正的随机数，其随机性和均匀性都是很好的、⽽且是取之不尽⽤之不竭的。但是⽤物理⽅法产⽣的随机数序列⽆法重复实现，不能进⾏程序复算，给验证结果带来很⼤困难。⽽且，物理随机数发⽣器的稳定性经常需进⾏检查和维修，费⽤昂贵。因此，⼤⼤降低了这种⽅法的使⽤价值[3]。
1.1.2 数学⽅法
在计算机上⽤数学⽅法产⽣均匀随机数是指按照⼀定的计算⽅法⽽产⽣的数列，它们具有类似于均匀随机变量的独⽴抽样序列的性质，这些数既然是依照确定算法产⽣的，便不可能是真正的随机数，因此常把⽤数学⽅法产⽣的随机数称为伪随机数。⽤数学⽅法在计算机上产⽣随机数是⽬前使⽤最⼴、发展很快的⼀类⽅法。他的特地啊是占⽤内存少、速度快⼜便于复算。
1.1.3 同余法
⽬前应⽤最⼴泛的随机数产⽣的⽅法是同余法。它是由Lehmer 在1951年提出的。此⽅法利⽤数论中的同于运算来产⽣随机数，故称为同余发⽣器。它包括混合同余发⽣器和乘同于发⽣器。它是⽬前使⽤最普遍、发展迅速的产⽣随机数的数学⽅法。
基本⽅法如下：
()()()101 mod mod / 1,2,n n n
n x ax b m x ax c m u x m n -=+??=+??==? （1.1）
其中m 为模数，a 为乘⼦(乘数)，c 为增量（加数），且n x ，m ，a ，c 均为⾮负整数。12,,,k u u u 就是k 个[0,1]上的伪随机数。按（1.1）给出的为随机数⽣成⽅法称为线性同余法。当0b >时称为混合同余法，0b =时称为乘同余法，当0b =且模是素数时，称为素数模乘同余法。线性同余法产⽣的随机数序列具有周期性，其最⼤周期为m 。
著名的Coveyou 与Macpherson 混合同于发⽣器为：
15351355 1 mod(2)/2 1,2,n n n n x x u x n -?=+??==??
Kobayashi 混合同余器是
31131314159269453806245 mod(2)/2 1,2,n n n n x x u x n -?=+??==??
1.2 ⾮均匀随机数⽣成
⾮均匀随机数的⽣成是另⼀类随机数⽣成⽅法，是基于[0,1]上的均匀随机数产⽣给定概率密度函数的算法。如正态分布随机数发⽣器，指数分布随机数发⽣器等，可以利⽤分布的特殊性质，针对性极强。在⼀些实际问题或统计模拟时需要上万，⼏⼗万，甚⾄上百万个⾮均匀随机数。这时，算法的速度就显得⼗分重要。⽬前产⽣⾮聚云随机数的⼀般⽅法有：直接抽样法（反函数法）、变换抽样法、值序抽样法和舍选抽样法。下⾯针对反函数变换法和舍选法进⾏介绍。
1.2.1 反函数变换法
设随机变量X 的分布函数是()F ，且F 是连续的。那么，()U F X =是随机变量，其分布是(0,1)上均匀分布。于是，有产⽣具有给
定分布函数()F 的⾮均匀随机数X 的算法:
1=(),~(0,1)
X F U U U - 这⾥(0,1)U 是(0,1)上的均匀分布，1()F - 是F 的反函数。是⼀个通⽤算法，对任何连续分布函数的反函数变换法，简称反函数变换法。该⽅法的缺点是：并不是任何分布函数都有快速算法。对常见分布，如正态分布，F 分布等分布函数仅可以数值计算，没有解析表达式。
1.2.2 舍选法
舍选法是⾮常有⽤的随机数⽣成⽅法。其应⽤最为⼴泛，研究成果也最多。设()f 是随机变量X 的密度函数，如何⽣成有密度函数()f 的随机数是随机数⽣成中最重要的和最基本的问题。最普遍应⽤的⽅法就是舍选法。⽅法如下:
1.存在函数()M 满⾜以下条件:
()()(),(),()M x f x M x M x dx C m x C
∞
-∞≤=<∞=
（1.2） ()m 也是密度函数，⽽且密度为()m 的随机变量Y 容易⽣成。
2.⽣成随机变量~()y m 和随机数~(0,1)V U ;
3.若()()M Y V f Y <，则X Y =，否则转2;
4.输出X , X 的密度是()f
满⾜(1.2)式的函数()M 叫做函数f 的优函数。不附加任何条件的优函数是容易到的。
1.3 随机数⽣成实例
运⽤反函数法，⽤(0,1)U ⽣成exp(1)的100个随机数。
因为指数分布的分布函数为：
1 0()0 0
x e x F x x -?-≥=?
()F x --=-λ。
利⽤R 语⾔进⾏随机数⽣成。反函数法⽣成100个1=λ的指数分布： set.seed(100)
FInverse<-function(y,lambda=1){
-log(1-y)/lambda
}
randomSequence<-FInverse(runif(100,0,1))
plot(randomSequence,type="l")
为检验利⽤反函数法⽣成的随机数是否为指数分布，进⾏K-S 检验，检验结果如下：
从输出结果可以看出0.26720.05p = ，故不能拒绝反函数法⽣成的序列来⾃于exp(1)。
2． MCMC 算法
马尔可夫链蒙特卡罗(MCMC)⽅法就是其中⼀种应⽤⼴泛的模拟⽅法。MCMC ⽅法最早由蒙特卡罗(MC)⽅法发展⽽来，蒙特卡罗⽅法是⼀种根据⼤数定律重复模拟得到近似概率的过程。
MCMC 的发展最早是由统计物理学家Metropolis 等⼈将马⽒链引⼊MC ⽅法中，解开了普通Monte Carlo ⽅法不改变抽样分布的瓶颈，使之能解决许多复杂的物理系统。在这之后随着计算机技术的蓬勃发展，极⼤的推动了MCMC ⽅法的⼴泛应⽤，在统计学领域中，极⼤似然估计、⾮参数估计、贝叶斯统计学等⽅向上的很多问题得到了有效解决。
2.1 MCMC 的基本思想
MCMC ⽅法的基本思想是通过构建⼀个平稳分布为()x π的Markov 能得到服从()x π的样本, 基于此样本对()x π的分布特征进⾏各种统计推断。
2.1.1 Markov 链
若随机变量序列(0)(1)(){,,,,}n X X X 满⾜
(1)(0)(0)(1)(1)()()(1)()()P(|,,)
P(|)n n n n n n X x X x X x X x X x X x ++≤====≤=
则称此随机变量序列为Markov 链。即在已知“现在”的条件下，“将来”与“过去”是相互独⽴的，“将来”仅与“现在”有关，与“过去”⽆关。
2.1.2 Markov 链的转移核及平稳分布
设(0)(1)(2){,,}X X X 为Markov 链。其⼀步转移概率函数为：
(1)()(,)P()P(|)t t p x x x x X x X x +'''→=== （离散）
或对于任⼀可测集B ?∈B ，有
P()(,)B
x B p x x dx ''→=? (连续)
P() , ??称为该Markov 链的转移核，通常假定P() , ??与时间t ⽆关(即时齐Markov 链)。
如果分布()x π满⾜：
(,)()()p x x x dx x ππ''=?
则称()x π为转移核P() , ??的平稳分布。
2.1.3 MCMC ⽅法步骤
MCMC ⽅法可概括为如下三个基本步骤：
(1)在样本空间
x 上选择“合适”的Markov 链，使()x π为其相应的平稳分布； (2)由x 中某⼀点()()00X x =出发，利⽤(1)中的Markov 链产⽣序列() ()()12,,,n X X X ；
(3)对给定的m 和充分⼤的n ，序列()()()12,,,m m n X X X ++ 可认为是来⾃()x π的遍历样本，可基此对()x π或其数字特征进⾏估计。任⼀函数()f x 的期望估计为：
()()
11?n t t m E f f X n m π=+=-∑ MCMC ⽅法的重点是构造转移核, 使给定的概率分布()x π为相应的Markov 链的平稳分布。
2.2 MCMC 算法实例
⽤MCMC 算法的办法⽣成(2,7)Beta 。
设初始点为0.1，利⽤(0.2,0.2)U -的随机数进⾏模拟变换。
N <- 3000 #抽样个数bootstrap检验方法
x <- c()
x0 <- 0.1 # 初始值
x[1] <- x0
k <- 0 #k表⽰拒绝转移的次数
u <- runif(N) #抽取均匀分布随机数
for (i in 2:N) {
y <- x[i-1]+runif(1,-0.2,0.2)
if(0
if (u[i] <=(dbeta(y,2,7)/dbeta(x[i-1],2,7)))
x[i] <- y
else {
x[i]<- x[i-1]
k <-k + 1
}
}
else
x[i]<-x0
}
Beta的分布，需要对x进⾏采样，采得到的x是⾼度⾃相关的，要得到(2,7)
样的步长通过⾃相关系数可确定，从下图可看出⾃相关系数在延迟30阶后均在⼆倍标准差内，故确定步长为30。acf(x)
进⼀步判断x是否达到了⼀个平稳分布，做出x的样本路径。
plot(x,type="l")
可看出，样本x已是平稳的，没有漂移的情况发⽣。可以认为这个马⽒链的构造是成功的。
开始采样，从1000开始，每隔30步取⼀个数，构成新的数列y。我们对y 与x分别作K-S检验，有
y<-x[seq(1000,N,by=30)]
Beta的从输出结果可以看出0.76160.05
p ，故不能拒绝序列来⾃于(2,7)
原假设。
3．EM 算法
在统计领域⾥，主要有两⼤类计算问题，⼀类是极⼤似然估计的计算，另⼀类是Bayes计算。从计算⽅法上讲，⼆者是可以合并讨论的，因极⼤似然估计的计算类似于Bayes的后验众数的计算。Bayes计算⽅法已经有很多，⼤体上可以分为两⼤类：⼀类是直接应⽤于后验分布以得到后验均值或后验众数的估计，以及这种估计的渐进⽅差或其近似。另⼀类算法总称为数据添加算法。这是近年发展很快且应⽤很⼴的⼀种算法，它不是直接对复杂的后验分布进⾏极⼤化或进⾏模
拟，⽽是在观测数据的基础上加上⼀些“潜在数据”，从⽽简化计算并完成⼀系列简单的极⼤化或模拟，该“潜在数据”可以使“缺损数据”或未知参数。
3.1 EM 算法基本原理
设我们能观测的数据Y ,θ关于Y 的后验分布(|)p Y θ很复杂，难以直接进⾏各种统计计算，加⼊我们能假定⼀些没有能观测到的潜在数据Z 为已知（譬如，Y 为某变量的结尾观测值，⽽Z 为该变量的真值），则可能得到⼀个关于θ的简单的添加后验分布(|,)p Y Z θ，利⽤(|,)p Y Z θ的简单性我们可进⾏各种统计计算，如极⼤化，抽样等，然⽽回过头来，我们⼜可以对Z 的假定做检查和改进，如此进⾏，我们就将⼀个复杂的极⼤化或抽样问题转变为⼀系列简单的极⼤化或抽样。
3.2 EM 算法基本思想与基本步骤
EM 算法是⼀种迭代⽅法，最初是由Dempster 等提出，并主要⽤来求后验分布的众数（即极⼤似然估计），它的每⼀次迭代由两步组成：E 步（求期望）和M 步（极⼤化）。⼀般地，以(|)p Y θ表⽰θ的基于观测数据的后验分布密度函数，
成为观测后验分布，
(|,)p Z Y θ表⽰在给定θ和观测数据Y 下潜在数据Z 的条件分布密度函数。⽬的在于计算观测后验分布(|)p Y θ的众数，于
是，EM 算法如下进⾏。记()i θ为第1i +次迭代开始时后验众数的估计值，则第1i +次迭代的两步为：
E 步：将(|,)p Y Z θ或log (|,)p Y Z θ关于Z 的条件分布求期望，从⽽把Z 积掉，即
()()()(|,)E [log (|,)|,]

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。