中国科学技术大学：《实用统计软件》课程课件讲义（统计计算与软件）第五讲 Monte Carlo积分和方差减少技术

Lecture5:Monte Carlo积分与方差减少技术 张伟平 Sunday 27th September,2009

Lecture 5: Monte Carlo »©Üê~E‚ ‹ï² Sunday 27th September, 2009

Contents 1 Monte Carlo Integration and Variance Reduction y 1.1 Monte Carlo Integration...............·..·· 1.1.1 Simple Monte Carlo estimator 3 1.1.2 Variance and Efficiency 12 l.2 Variance Reduction···.······-·····.···· 14 1.3 Antithetic Variables 15 1.4 Control Variates .。·。。。。。。。·。。。·24 1.4.1 Antithetic variate as control variate 30 l.4.2 Several control variates·..·.··········· 31 1.5 Importance sampling··············· 35 1.6 Stratified Sampling.············· 43 1.7 Stratified Importance Sampling ....... 47 Previous Next First Last Back Forward 1

Contents 1 Monte Carlo Integration and Variance Reduction 1 1.1 Monte Carlo Integration . . . . . . . . . . . . . . . . . . . . 3 1.1.1 Simple Monte Carlo estimator . . . . . . . . . . . . 3 1.1.2 Variance and Efficiency . . . . . . . . . . . . . . . . 12 1.2 Variance Reduction . . . . . . . . . . . . . . . . . . . . . . . 14 1.3 Antithetic Variables . . . . . . . . . . . . . . . . . . . . . . 15 1.4 Control Variates . . . . . . . . . . . . . . . . . . . . . . . . 24 1.4.1 Antithetic variate as control variate . . . . . . . . . 30 1.4.2 Several control variates . . . . . . . . . . . . . . . . 31 1.5 Importance sampling . . . . . . . . . . . . . . . . . . . . . . 35 1.6 Stratified Sampling . . . . . . . . . . . . . . . . . . . . . . . 43 1.7 Stratified Importance Sampling . . . . . . . . . . . . . . . . 47 Previous Next First Last Back Forward 1

Chapter 1 Monte Carlo Integration and Variance Reduction 蒙特卡罗(Monte Carlo)积分是一种基于随机抽样的统计方法.蒙特卡罗方法 其实也只是对一种思想的泛称，只要在解决问题时，利用产生大量随机样本， 然后对这些样本结果进行概率分析，从而来预测结果的方法，都可以称为蒙特 卡罗方法。 比如要求圆周率π的值，最著名的就是中学时学过的”割圆法”（刘徽（魏晋， 3.1416),祖冲之（南北朝，3.1415926<T<3.1415927).现在也可以使用蒙特 卡罗积分方法：由概率论知 若ruX,Yii.dU0,1),则Px2+y2≤1)= Previous Next First Last Back Forward 1

Chapter 1 Monte Carlo Integration and Variance Reduction ÑAk¤(Monte Carlo)»©¥ò´ƒuëŃ⁄Oê{. ÑAk¤ê{ Ÿ¢ èê¥Èò´géç°, êá3)˚ØKû, |^)å˛ëÅß ,￾È ˘ (J?1V«©¤, l 5˝ˇ(Jê{, —å±°èÑA k¤ê{. 'Xᶱ«πä, ÅÕ¶“¥•ÆûÆL”{”(4†(üA, 3.1416), y¿É(Hä, 3.1415926< π <3.1415927)). y3è屶^ÑA k¤»©ê{: dV«ÿ e r.v. X, Y i.i.d U(0, 1), K P(X2 + Y 2 ≤ 1) = π 4 Previous Next First Last Back Forward 1

因此，开=4P(X2+Y2≤1)≈4#{x2+2≤1}/m. ·n=1000,元=3.168 ●n=100.000,元=3.14312. ●n=107,元=3.141356. Previous Next First Last Back Forward 2

œd, π = 4P(X2 + Y 2 ≤ 1) ≈ 4#{x 2 + y 2 ≤ 1}/n. • n = 1000, ˆπ = 3.168. • n = 100, 000, ˆπ = 3.14312. • n = 107 , ˆπ = 3.141356. Previous Next First Last Back Forward 2

1.1 Monte Carlo Integration 假设g是一个可积函数，我们要计算∫公g(x)dc.回忆在概率论中，若随机变 量X的密度为f(x),则随机变量Y=9(X)的期望为 )- g(x)f(x)dz 如果从X的分布中产生了一些随机数，则Eg(X)的无偏估计就是相应的样本平 均值. 1.1.1 Simple Monte Carlo estimator 考虑估计6=g(c)dz.若X1,·,Xm为均匀分布U(0,1)总抽取的样本，则 由强大数律知 =9m=1∑g(X) m i=1 以概率1收敛到期望Eg(X).因此6g(r)dc的简单的Monte Carlo估计量 为gm(X) Previous Next First Last Back Forward 3

1.1 Monte Carlo Integration bg¥òá建Í, ·ÇáOé R b a g(x)dx. ££3V«ÿ•, eëÅC ˛Xó›èf(x), KëÅC˛Y = g(X)œ"è Eg(X) = Z ∞ −∞ g(x)f(x)dx XJlX©Ÿ•) ò ëÅÍ, KEg(X)ÆO“¥ÉA² ˛ä. 1.1.1 Simple Monte Carlo estimator ƒOθ = R 1 0 g(x)dx. eX1, · · · , Xmè˛!©ŸU(0, 1)oƒ, K dråÍÆ θˆ = gm(X) = 1 m Xm i=1 g(Xi) ±V«1¬Òœ"Eg(X). œd R 1 0 g(x)dx{¸Monte Carlo O˛ ègm(X). Previous Next First Last Back Forward 3

例1：（简单的Monte Carlo积分）计算 -dr 的简单onte Carlo估计以及与积分值相比较. TExample m<-10000 x <-runif(m) theta.hat <-mean(exp(-x)) print(theta.hat) print(1 -exp(-1)) #[1]0.6355289 #[1]0.6321206 ↓Example 估计为0÷0.6355，而积分值为0=1-e-1÷0.6321. 若要计算9(x)d,此处a<b为有限数.则作一积分变量代换使得积分限 从0到1.即作变换y=(x-a)/(b-a),因此 9a= g(y(b-a)+a)(b-a)dy Previous Next First Last Back Forward 4

~1: ({¸Monte Carlo »©) Oé θ = Z 1 0 e −x dx {¸Monte CarloO±9Ü»©äÉ'. ↑Example m <- 10000 x <- runif(m) theta.hat <- mean(exp(-x)) print(theta.hat) print(1 - exp(-1)) #[1] 0.6355289 #[1] 0.6321206 ↓Example Oèθˆ .= 0.6355, »©äèθ = 1 − e−1 .= 0.6321. eáOé R b a g(x)dx, d?a < bèkÅÍ. Käò»©C˛ìܶ»©Å l01. =äCÜy = (x − a)/(b − a), œd Z b a g(x)dx = Z 1 0 g(y(b − a) + a)(b − a)dy Previous Next First Last Back Forward 4

另外一种做法就是利用均匀分布U(α，b),即 例2：简单Monte Carlo积分（续）计算 o-l erds 的Monte Carlo估计并和积分值相比较. TExample m<-10000 x <-runif(m,min=2,max=4) theta.hat <mean(exp(-x))*2 print(theta.hat) print (exp(-2)-exp(-4)) #[1]0.1172158 #[1]0.1170196 Previous Next First Last Back Forward 5

, ò´â{“¥|^˛!©ŸU(a, b), = Z b a g(x)dx = (b − a) Z b a g(x) 1 b − a dx ~2: {¸Monte Carlo »©(Y) Oé θ = Z 4 2 e −x dx Monte CarloOø⁄»©äÉ'. ↑Example m <- 10000 x <- runif(m, min=2, max=4) theta.hat <- mean(exp(-x)) * 2 print(theta.hat) print(exp(-2) - exp(-4)) #[1] 0.1172158 #[1] 0.1170196 Previous Next First Last Back Forward 5

即，估计的值为8÷0.1172，而真值为0=e-2-e-4÷0.1170. 总结一下，积分∫2g(x)dc的简单Monte Carlo估计方法为 1.从均匀分布U(a,b)中产生i.i.d样本X1,·,Xm 2.计算gm(X=品9(X) 3.a=(b-a)(gm(X)- 例3：Monte Carlo积分，无穷积分比如使用Monte Carlo积分方法计算标 准正态的分布函数 Φ()= =e-号t /2 首先我们不能直接使用以前的方法（因为积分区间无界），但是我们可以将 此问题分为两种情形：￡>0和x≤0.若x>0,注意到对标准正态分布，积分 Previous Next First Last Back Forward 6

↓Example =, Oäèθˆ .= 0.1172, ˝äèθ = e−2 − e−4 .= 0.1170. o(òe, »© R b a g(x)dx{¸Monte CarloOê{è 1. l˛!©ŸU(a, b)•)i.i.dX1, · · · , Xm; 2. Oégm(X) = 1 m g(Xi) 3. θˆ = (b − a)(gm(X)). ~3: Monte Carlo»©, 𻩠'X¶^Monte Carlo»©ê{OéI O©ŸºÍ Φ(x) = Z x −∞ 1 √ 2π e − t 2 2 dt. ƒk·ÇÿUܶ^±cê{(œè»©´mÃ.), ¥·Çå±Ú dØK©è¸´ú/: x > 0⁄x ≤ 0. e x > 0, 5øÈIO©Ÿ, »© Previous Next First Last Back Forward 6

区间可以分为(-0,0)和(0，)，因此只需要计算积分8=片e号d止即可.故 可以使用之前的方法.但是需要从均匀分布U(O,x)中抽取样本，若x发生变化， 则均匀分布也就变化了.若要求从U(0,1)中抽样，则可以作变换y=t/x,则 0=f se-(wdy 因此，0=Eyze-(xY)],其中Y~U(0,1).从而产生U(0,1)的ii.d随机 数u1,,,um,则的估计为 =gm=∑e-m2 n i1 此时对x>0,Φ(x)的估计为0.5+/√2；对x≤0，计算Φ(x)=1-(-x) TExample x<-seq(.1,2.5,1 ength=10) m<-10000 u <runif(m) Previous Next First Last Back Forward 7

´m屩è(−∞, 0)⁄(0, x),œdêIáO黩θ = R x 0 e − t 2 2 dt =å.  屶^Écê{. ¥Iál˛!©ŸU(0, x)•ƒ, exu)Cz, K˛!©Ÿè“Cz . eá¶lU(0, 1)•ƒ, Kå±äCÜy = t/x, K θ = Z 1 0 xe−(xy) 2 dy œd, θ = EY [xe−(xY ) 2 ], Ÿ•Y ∼ U(0, 1). l )U(0, 1)i.i.dëÅ Íu1, · · · , um, K θOè θˆ = gm(u) = 1 m Xm i=1 xe−(xu) 2 . dûÈx > 0,Φ(x)Oè0.5 + θ/ˆ √ 2π; Èx ≤ 0, OéΦ(x) = 1 − Φ(−x). ↑Example x <- seq(.1, 2.5, length = 10) m <- 10000 u <- runif(m) Previous Next First Last Back Forward 7

cdf <-numeric(length(x)) for (i in 1:length(x)){ g<-x[1]*exp(-(u*x[1])2/2) cdf [i]<-mean(g)/sqrt(2 pi)+0.5 Phi <-pnorm(x) print(round(rbind(x,cdf,Phi),3)) ↓Example 例4：Monte Carlo积分，无穷积分（续）对上例，我们可以使用另外一种方 式(hit-or-miss).记Z~N(0,1),则对任何实数x有 Φ(x)=P(Z≤x)=EI(Z≤x): 从而从标准正态中产生随机样本z1,·,之m后，即可得到Φ(x)的估计为 =e≤ 三1 Previous Next First Last Back Forward 8

cdf <- numeric(length(x)) for (i in 1:length(x)) { g <- x[i] * exp(-(u * x[i])^2 / 2) cdf[i] <- mean(g) / sqrt(2 * pi) + 0.5 } Phi <- pnorm(x) print(round(rbind(x, cdf, Phi), 3)) ↓Example ~4: Monte Carlo»©, ð»©(Y) È˛~, ·Ç屶^, ò´ê ™(hit-or-miss). PZ ∼ N(0, 1), KÈ?¤¢Íxk Φ(x) = P(Z ≤ x) = EI(Z ≤ x). l lIO•)ëÅz1, · · · , zm￾, =åΦ(x)Oè Φ( ˆ x) = 1 m Xm i=1 I(zi ≤ x). Previous Next First Last Back Forward 8