trimmed.mse <-function(n,m,k,p){ #MC est of mse for k-level trimmed mean of #contaminated normal pN(0,1)+(1-p)N(0,100) tmean <numeric(m) for (i in 1:m){ sigma <sample(c(1,10),size n, replace TRUE,prob =c(p,1-p)) x <-sort(rnorm(n,0,sigma)) tmean[i]<-sum(x[(k+1):(n-k)])/(n-2*k) mse.est <-mean(tmean2) se.mse <-sqrt(mean((tmean-mean(tmean))2))/sqrt(m) return(c(mse.est,se.mse)) ] for (k in 0:K){ mse[k+1,1:2]<-trimmed.mse(n=n,m=m,k=k,p=1.0) mse[k+1,3:4]<-trimmed.mse(n=n,m=m,k=k,p=.95) mse[k+1,5:6]<-trimmed.mse(n=n,m=m,k=k,p=.9) round(n*mse,3) ↓Ccdn Previous Next First Last Back Forward 9
trimmed.mse <- function(n, m, k, p) { #MC est of mse for k-level trimmed mean of #contaminated normal pN(0,1) + (1-p)N(0,100) tmean <- numeric(m) for (i in 1:m) { sigma <- sample(c(1, 10), size = n, replace = TRUE, prob = c(p, 1-p)) x <- sort(rnorm(n, 0, sigma)) tmean[i] <- sum(x[(k+1):(n-k)]) / (n-2*k) } mse.est <- mean(tmean^2) se.mse <- sqrt(mean((tmean-mean(tmean))^2)) / sqrt(m) return(c(mse.est, se.mse)) } for (k in 0:K) { mse[k+1, 1:2] <- trimmed.mse(n=n, m=m, k=k, p=1.0) mse[k+1, 3:4] <- trimmed.mse(n=n, m=m, k=k, p=.95) mse[k+1, 5:6] <- trimmed.mse(n=n, m=m, k=k, p=.9) } round(n*mse,3) ↓Code Previous Next First Last Back Forward 9
结果表明，均值的稳健估计（截尾均值估计）在总体分布被污染时能降低MSE. 1.2 Estimating a confidence level 在应用统计中经常遇到的一个问题是估计总体的分布.比如，许多常用的统计 推断方法和工具都是基于正态性假设下的.而在实际中，总体分布非正态是经 常的，估计量的分布可能不知道或者没有显示表示，此时，Monte Carlo方法 则可以用来进行统计推断」 假设(U,V)是未知参数的置信区间，则统计量U,V的分布都依赖于抽样分 布X的分布Fx·置信水平就是区间(U,V)能够覆盖真值的概率.因此估计置 信水平就是一个积分估计问题， 比如考虑方差的置信区间估计问题.标准的方法是对正态性假设很敏感的， 在数据（样本）偏离正态分布时，我们来使用Monte Carlo方法估计真实的置信 水平.首先看正态假定下方差的置信区间估计（标准方法）： Previous Next First Last Back Forward 10
(JL², ˛äËO( ó˛äO)3oN©Ÿ¿/ûU¸$ MSE. 1.2 Estimating a confidence level 3A^⁄O•²~ëòáØK¥OoN©Ÿ. 'X, Nı~^⁄O ̉ ê{⁄Û‰—¥ƒu5be. 3¢S•, oN©Ÿö¥² ~, O ˛©ŸåUÿ½ˆvkw´L´, dû, Monte Carlo ê{ Kå±^5?1 ⁄Ỏ. b(U, V )¥ôÎÍθò&´m, K⁄O˛U, V ©Ÿ—ù6uƒ © ŸX©ŸFX. ò&Y²“¥´m(U, V )U CXθ˝äV«. œd Oò &Y²“¥òứOØK. 'Xƒê ò&´mOØK. IOê{¥È5bÈØa, 3Í‚( ) †l©Ÿû, ·Ç5¶^Monte Carloê{O˝¢ò& Y². ƒkwb½ eê ò&´mO(IOê{): Previous Next First Last Back Forward 10
例4方差的置信区间假设X1,·,Xn,ii.d~N(4,a2),n之2，S2为 样本方差，则由 (n-1)S2 知道a2的一个100(1-a)置信上界为 (0,(n-1)s2/X2(n-1) 从而从正态总体N(0,σ2=4)中随机抽取n=20个样本，则计算g2的%95置信 上界的程序如下 FCode n<-20 alpha <-.05 x <-rnorm(n,mean=0,sd=2) UCL <-(n-1)*var(x)/qchisq(alpha,df=n-1) Code 我们可以使用图来直观上判断经验的置信水平： fccd加 Previous Next First Last Back Forward 11
~ 4: ê ò&´m bX1, · · · , Xn, i.i.d ∼ N(µ, σ2 ), n ≥ 2, S 2è ê , K d (n − 1)S 2 σ2 ∼ χ 2 n−1 σ 2òá100(1 − α)ò&˛.è (0, (n − 1)S 2 /χ2 α(n − 1)) l loNN(0, σ2 = 4)•ëŃn = 20á , KOéσ 2%95ò& ˛. ßSXe ↑Code n <- 20 alpha <- .05 x <- rnorm(n, mean=0, sd=2) UCL <- (n-1) * var(x) / qchisq(alpha, df=n-1) ↓Code ·Ç屶^„5Ü*˛‰²ò&Y²: ↑Code Previous Next First Last Back Forward 11
m<-100000 ucls<-numeric(m) for(i in 1:m){ x <-rnorm(n,mean=0,sd=2) ucls[i]<-(n-1)*var(x)/gchisq(alpha,df=n-1) ind<-ucls>4 cov.rate<-cumsum(ind)/1:m plot(2:m,cov.rate[-1],type="1") abline(h=0.95) Code 经验的置信水平是通过模拟，对理论的置信水平进行估计，其一股做法如下 Previous Next First Last Back Forward 12
m<-100000 ucls<-numeric(m) for(i in 1:m){ x <- rnorm(n, mean=0, sd=2) ucls[i] <- (n-1) * var(x) / qchisq(alpha, df=n-1) } ind<-ucls>4 cov.rate<-cumsum(ind)/1:m plot(2:m,cov.rate[-1],type="l") abline(h=0.95) ↓Code ²ò&Y²¥œL [, Ènÿò&Y²?1O. ŸòÑâ{Xe Previous Next First Last Back Forward 12
假设X~Fx,感兴趣的参数为0.则对j=1,·,m: 1.产生第个随机样本，…，日. 2.计算基于第j个样本的置信区间C 3.计算=I(0∈C) 4.计算经验的置信水平可=品∑严1 例5：置信水平的Monte Carlo估计在上一个例子中，我们使用了for循环 来实现计算a的m个置信上界.我们也可以使用replicatei函数： 下Cde n<-20 alpha <-.05 UCL <replicate(1000,expr = x <-rnorm(n,mean 0,sd =2) (n-1)*var(x)/qchisq(alpha,df n-1) ]) Previous Next First Last Back Forward 13
bX ∼ FX, a,ÎÍèθ. KÈj = 1, · · · , m: 1. )1jáëÅ x (j) 1 , · · · , x (j) n . 2. Oéƒu1já ò&´mCj 3. Oéyj = I(θ ∈ Cj ) 4. Oé²ò&Y²y¯ = 1 m Pm j=1 yj . ~ 5: ò&Y²Monte CarloO 3˛òá~f•, ·Ç¶^ for ÃÇ 5¢yOéσ 2máò&˛.. ·Çè屶^replicateºÍ: ↑Code n <- 20 alpha <- .05 UCL <- replicate(1000, expr = { x <- rnorm(n, mean = 0, sd = 2) (n-1) * var(x) / qchisq(alpha, df = n-1) } ) Previous Next First Last Back Forward 13