中国科学技术大学：《实用统计软件》课程课件讲义（统计计算与软件）第六讲 Monte Carlo方法在统计推断中的应用

Lecture6:统计推断中的Monte Carlo方法 张伟平 Monday 12th October,2009

Lecture 6: ⁄Ỏ•Monte Carlo ê{ ‹ï² Monday 12th October, 2009

Contents 1 Monte Carlo Methods in Inference 1 1.1 Monte Carlo Methods for Estimation............. 1.1.1 Monte Carlo Estimation and Standard Error ... 1.l.2 Estimation of MSE...................5 l.2 Estimating a confidence level··.··.·····.····· 10 l.3 Monte Carlo Methods for Hypothesis Tests·········l6 l.4 Empirical Type I error rate.··.·····.······· 17 1.4.1 Power of a Test.·.··..·.·.···...···24 l.4.2 Power Comparisons··· 30 1.5 Application:"Count Five"Test for Equal Variance ....34 Previous Next First Last Back Forward 1

Contents 1 Monte Carlo Methods in Inference 1 1.1 Monte Carlo Methods for Estimation . . . . . . . . . . . . . 1 1.1.1 Monte Carlo Estimation and Standard Error . . . . 2 1.1.2 Estimation of MSE . . . . . . . . . . . . . . . . . . . 5 1.2 Estimating a confidence level . . . . . . . . . . . . . . . . . 10 1.3 Monte Carlo Methods for Hypothesis Tests . . . . . . . . . 16 1.4 Empirical Type I error rate . . . . . . . . . . . . . . . . . . 17 1.4.1 Power of a Test . . . . . . . . . . . . . . . . . . . . . 24 1.4.2 Power Comparisons . . . . . . . . . . . . . . . . . . 30 1.5 Application: “Count Five” Test for Equal Variance . . . . . 34 Previous Next First Last Back Forward 1

Chapter 1 Monte Carlo Methods in Inference Monte Carlo方法可以指统计推断或者数值分析中任何使用随机数模拟的方 法.本章我们讨论一些统计推断中的Monte Carlo方法.Monte Carlo方法可 以用于估计统计量分布中的参数，均方误差，分位数等，也可以用于估计置信 区间覆盖真实参数的概率（置信水平），假设检验中检验的一型错误率，功效等 等 1.1 Monte Carlo Methods for Estimation 假设X1,,X为从总体X中抽取的随机样本，参数的估计量记为 a=X1,…,Xn) Previous Next First Last Back Forward

Chapter 1 Monte Carlo Methods in Inference Monte Carloê{å±ç⁄Ỏ½ˆÍ䩤•?¤¶^ëÅÍ[ê {. Ÿ ·Ç?ÿò ⁄Ỏ•Monte Carloê{. Monte Carloê{å ±^uO⁄O˛©Ÿ •ÎÍ, ˛êÿ, ©†Í, èå±^uOò& ´mCX˝¢ÎÍV«(ò&Y²), bu•uò.Üÿ«, ı . 1.1 Monte Carlo Methods for Estimation bX1, · · · , XnèloNX•ƒëÅ, ÎÍθO˛Pè θˆ = θˆ(X1, · · · , Xn) Previous Next First Last Back Forward 1

记x=(z1,·,xn)T∈R”,以及x),x②)，，为从总体X中抽取的一列独 立的随机样本观测值则有关于的性质，可以通过估计值序列()，·，x),方= 1,2,··来研究 1.1.1 Monte Carlo Estimation and Standard Error 例1：假设X1,X2i.i.d~N(0,1),估计EX1-X2 显然，0=EX1一X2的Monte Carlo估计可用通过从标准正态分布 中产生m个样本x)=(c9),x),方=1…，m.然后计算)= z》-x91,j=1,…,m.以及0的估计 m j=1 j=1 在R中很容易实现： m<-1000 Previous Next First Last Back Forward 2

Px = (x1, · · · , xn) T ∈ Rn, ±9 x (1), x(2) , · · · èloNX•ƒò’ ·ëÅ*ˇä. K k'uθˆ5ü, 屜LOäSθˆ(x (j) 1 , · · · , x (j) n ), j = 1, 2, · · · 5Ôƒ. 1.1.1 Monte Carlo Estimation and Standard Error ~ 1: bX1, X2 i.i.d ∼ N(0, 1), OE|X1 − X2|. w,, θ = E|X1 − X2|Monte CarloOå^œLlIO©Ÿ • )má x (j) = (x (j) 1 , x (j) 2 ), j = 1, · · · , m. ,￾Oé θˆ(j) = |x (j) 1 − x (j) 2 |, j = 1, · · · , m. ±9θO θˆ = 1 m Xm j=1 θˆ(j) = 1 m Xm j=1 |x (j) 1 − x (j) 2 |. 3R•ÈN¥¢y: ↑Code m <- 1000 Previous Next First Last Back Forward 2

g <-numeric(m) for (i in 1:m){ x <-rnorm(2) g[i]<-abs(x[1]-x[2]) est <-mean(g) est Code 我们也可以计算出EX1-X2l=2/V示三1.128379以及方差Var(X1- X2)=2-4/m.因此E8=2/V元，Var(=2-4//m. 对标准差的Monte Carlo估计我们可以从一般场合出发讨论.由于样本 量为n的样本均值又的标准差为vWar(X)/m,当X的分布未知时，可以使 用”Pug-in”法估计：由 n var(x)=(4-. n i=1 Previous Next First Last Back Forward 3

g <- numeric(m) for (i in 1:m) { x <- rnorm(2) g[i] <- abs(x[1] - x[2]) } est <- mean(g) est ↓Code ·Çèå±Oé—E|X1 − X2| = 2/ √ π .= 1.128379 ±9êV ar(|X1 − X2|) = 2 − 4/π. œdEθˆ = 2/ √ π, V ar(θˆ) = [2 − 4/π]/m. ÈIOMonte Carlo O·Çå±lòÑ|‹—u?ÿ. du ˛èn˛äX¯ IOè p V ar(X)/m, X©Ÿôû, 屶 ^”Plug-in”{O: d V ar ˆ (X) = 1 m Xm i=1 (xi − x¯) 2 . Previous Next First Last Back Forward 3

因此灭标准差的估计为 s()= a[品-时]=e- 或者也可以使用无偏估计量 1 因此，前例中的标准差估计为 = 计算程序如下 sqrt(sum((g-mean(g))2))/m #sd(g)/sqrt(m)#for unbiased estimator Code Previous Next First Last Back Forward 4

œdX¯IOOè seˆ (¯x) = 1 √ m  1 m (xi − x¯) 2 1/2 = 1 m (xi − x¯) 2 1/2 . ½ˆè屶^ÆO˛ seˆ (¯x) = 1 √ m  1 m − 1 (xi − x¯) 2 1/2 . œd, c~•IOOè seˆ (θˆ) = 1 m  Xm j=1 (θˆ(j) − θˆ) 2   . OéßSXe ↑Code sqrt(sum((g-mean(g))^2))/m #sd(g)/sqrt(m) #for unbiased estimator ↓Code Previous Next First Last Back Forward 4

1.1.2 Estimation of MSE Monte Carlo方法可以用于计算一个估计量的MSE.一个估计量的MSE定 义为MSE(©=E6-.如果从总体X中产生了m个样本x(1,·,x(m), 则e的MSE的Monte Carlo估计为 1 ME=(6)-)2, m J=1 其中U)=(x). 例2：使用Monte Carlo方法估计标准正态分布的截尾均值-1的MSE. 当样本存在异常点时，截尾的样本均值常常被用来估计总体的中心.假 设X1,·,Xn为一个随机样本，X()·,X(n)为相应的次序统计量，则一 个k水平的截尾样本均值为 1 n一 -= n-2k ∑xo i=k+1 Previous Next First Last Back Forward 5

1.1.2 Estimation of MSE Monte Carlo ê{å±^uOéòáO˛MSE. òáO˛MSE½ ¬èMSE(θˆ) = E[θˆ − θ] 2 . XJloNX•) máx (1) , · · · , x(m) , KθˆMSEMonte CarloOè MSEˆ = 1 m Xm j=1 (θˆ(j) − θ) 2 , Ÿ•θˆ(j) = θˆ(x (j) ). ~ 2: ¶^Monte Carloê{OIO©Ÿó˛äX¯ [−1]MSE. 3…~:û, ó˛ä~~^5OoN•%. b X1, · · · , XnèòáëÅ, X(1), · · · , X(n)èÉAgS⁄O˛, K ò ákY²ó˛äè X¯ [−k] = 1 n − 2k nX−k i=k+1 X(i) . Previous Next First Last Back Forward 5

本例中，目标参数为9=E了=EX-1)=0.记T=X-1,则其MSE的 一个Monte Carlo估计算法如下 1.通过如下步骤产生m个重复T),j=1,…,m: 。产生总体X的样本：x,, 。从小到大排序唱≤…≤名 ·计算T0)=高∑8 2.计算MSE MSE(T)=六∑m1(T)-P 六(T6)2. 实现程序如下 TCode n<-20 m<-1000 tmean <-numeric(m) Previous Next First Last Back Forward 6

~•, 8IÎÍèθ = EX¯ = EX¯ [−1] = 0. PT = X¯ [−1], KŸMSE òá Monte Carlo Oé{Xe 1. œLXe⁄½)má­ET (j) , j = 1, · · · , m: • )oNX: x (j) 1 , · · · , x (j) n . • lå¸S x (j) (1) ≤ · · · ≤ x (j) (n) • Oé T (j) = 1 n−2 Pn−1 i=2 x (j) (i) . 2. O éMSE MSEˆ (T) = 1 m Pm j=1(T (j) − θ) 2 = 1 m (T (j) ) 2 . ¢yßSXe ↑Code n <- 20 m <- 1000 tmean <- numeric(m) Previous Next First Last Back Forward 6

for (i in 1:m){ x <-sort(rnorm(n)) tmean[1]<-sum(x[2:(n-1)])/(n-2) mse <-mean(tmean"2) mse sqrt(sum((tmean mean(tmean))2))/m #s日 Code 截尾均值的MSE的估计为0.0504531(se÷0.007).样本均值的MSE为Var(X)/n= 1/20=0.05.另一方面，中位数本质上也是一种截尾均值，其截掉了除中间的 一个或两个点外的所有点.样本中位数的MSE估计如下 TCode n<-20 m<-1000 tmean <-numeric(m) for (i in 1:m){ x <-sort(rnorm(n)) tmean[i]<-median(x) ] Previous Next First Last Back Forward 7

for (i in 1:m) { x <- sort(rnorm(n)) tmean[i] <- sum(x[2:(n-1)]) / (n-2) } mse <- mean(tmean^2) mse sqrt(sum((tmean - mean(tmean))^2)) / m #se ↓Code ó˛äMSEOè0.0504531( ˆse .= 0.007). ˛ä MSEèV ar(X)/n = 1/20 = 0.05. ,òê°, •†Íü˛è¥ò´ó˛ä, ŸK ÿ•m òὸá: §k:. •†ÍMSEOXe ↑Code n <- 20 m <- 1000 tmean <- numeric(m) for (i in 1:m) { x <- sort(rnorm(n)) tmean[i] <- median(x) } Previous Next First Last Back Forward 7

mse <mean(tmean2) mse sqrt(sum((tmean mean(tmean))2))/m #se Code 从而样本中位数的MSE估计为0.075(s2=0.0086). 例3：比较标准正态分布与如下混合(“污染”)正态分布下的k水平戴尾均值 估计的MSE. pN(0,a2=1)+(1-p)N(0,a2=100) 我们写一个函数来对不同的k和P计算截尾均值X-的MSE.从混合正 态中产生样本时，要根据P(a=1)=p:P(a=10)=1-p来随机选 择a.注意正态随机数产生函数rnorm可以使用参数向量作为标准偏差.考 虑p=1.0,0.95,0.9以及k=0,1,·,n/2.因此程序如下 fccd加 n<-20;K<-n/2-1;m<-1000; mse <-matrix(0,n/2,6) Previous Next First Last Back Forward 8

mse <- mean(tmean^2) mse sqrt(sum((tmean - mean(tmean))^2)) / m #se ↓Code l •†ÍMSEOè0.075( ˆse = 0.0086). ~ 3: 'IO©ŸÜXe·‹(“¿/”)©ŸekY²ó˛ä OMSE. pN(0, σ2 = 1) + (1 − p)N(0, σ2 = 100) ·ÇòáºÍ5Èÿ”k⁄pOéó˛äX¯ [−k]MSE. l·‹ •) û, áä‚P(σ = 1) = p; P(σ = 10) = 1 − p5ëÅ¿ Jσ. 5øëÅÍ) ºÍrnorm屶^ÎÍï˛äèIO†.  ƒp = 1.0, 0.95, 0.9 ±9k = 0, 1, · · · , n/2. œdßSXe ↑Code n <- 20; K <- n/2 - 1; m <- 1000; mse <- matrix(0, n/2, 6) Previous Next First Last Back Forward 8