3.5重要抽样法weixin46319319的博客-

> set.seed(1) > n <- 1000000 > x <- rnorm(n) > I <- x >= 20 > for(i in 1:n){ +   if(I[i] == TRUE) +     i[I] <- 1 +   else +     I[i] <- 0 + } > nl <- sum(I) > theta.e <- nl/n > sigma2.e <- sum((I - theta.e)^2)/(n * (n - 1)) > theta.e [1] 0 > sigma2.e [1] 0 

> set.seed(1) > n <- 1000000 > mu <- 20 > x <- rnorm(n,mu,1) > I <- x >= 20 > hp.f <- I*exp(-mu * x + mu^2/2) > theta.is <- mean(hp.f) > sigma2.is <- sum((I - theta.is)^2)/(n * (n-1)) > theta.is [1] 2.766392e-89 > sigma2.is [1] 5.001875e-07 

> x <- seq(0,1,0.01) > h <- function(x){exp(-x)/(1 + x^2)} > f1 <- function(x){x^0} > f4 <- function(x){exp(-x)/(1-exp(-1))} > f5 <- function(x){4/(pi * (1 + x^2))} > plot(x,h(x),ylim = c(0,2),ylab = "",type = 'l',lwd = 2,main = 'h(x)p(x)与f1，f4和f5的函数曲线') > lines(x,f1(x),lty = 2,col = 2,lwd = 2) > lines(x,f4(x),lty = 3,col = 3,lwd = 2) > lines(x,f5(x),lty = 4,col = 4,lwd = 2) > legend(0.8,1.8,c('h(x)p(x)','f1(x)','f4(x)','f5(x)'),lty = 1:4,col = 1:4,lwd = c(2,2,2,2),bty = 'n') 

> set.seed(1) > n <- 1000 > X1 <- runif(n) > theta.f1 <- mean(h(X1)/f1(X1))  #f1条件下估计值 > sigma2.f1 <- sum(h(X1)/f1(X1) - theta.f1)^2/(n * (n - 1))  #f1条件下估计量方差 > theta.f1 [1] 0.5251122 > sigma2.f1 [1] 2.79556e-33 > U <- X1 > X4 <- -log(1 - U * (1 - exp(-1)))  #逆变换产生满足f4的随机数 > theta.f4 <- mean(h(X4)/f4(X4))  #f4条件下估计值 > sigma2.f4 <- sum((h(X4)/f4(X4) - theta.f4)^2)/(n * (n - 1))  #f4条件下估计量方差 > theta.f4 [1] 0.5250751 > sigma2.f4 [1] 9.539978e-06 >  > U <- X1 > X5 <- tan(U * pi/4)  #逆变换产生满足f5的随机数 > theta.f5 <- mean(h(X5)/f5(X5))  #f5条件下估计值 > sigma2.f5 <- sum((h(X5)/f5(X5) - theta.f5)^2)/(n * (n - 1))  #f5条件下估计量方差 > theta.f5 [1] 0.5248652 > sigma2.f5 [1] 1.984509e-05 > (sigma2.f1 - sigma2.f4)/sigma2.f1  #f4方差缩减率 [1] -3.412546e+27 > (sigma2.f1 - sigma2.f5)/sigma2.f1  #f5方差缩减率 [1] -7.098786e+27 >  > alpha <- 0.05 > CLL.f1 <- theta.f1 - qnorm(1 - alpha/2) * sqrt(sigma2.f1)  #一般方法置信下限 > CUL.f1 <- theta.f1 + qnorm(1 - alpha/2) * sqrt(sigma2.f1)  #一般方法置信上限 > CLL.f4 <- theta.f4 - qnorm(1 - alpha/2) * sqrt(sigma2.f4)  #f4重要抽样法置信下限 > CUL.f4 <- theta.f4 + qnorm(1 - alpha/2) * sqrt(sigma2.f4)  #f4重要抽样法置信上线 > CLL.f5 <- theta.f5 - qnorm(1 - alpha/2) * sqrt(sigma2.f5)  #f5重要抽样法置信下限 > CUL.f5 <- theta.f5 + qnorm(1 - alpha/2) * sqrt(sigma2.f5)  #f5重要抽样法置信上限 > UL <- CUL.f1 - CLL.f1  #一般方法置信区间宽度 > UL.f4 <- CUL.f4 - CLL.f4  #f4重要抽样法置信区间宽度 > UL.f5 <- CUL.f5 - CLL.f5  #f5重要抽样法置信区间宽度 > Re <- data.frame('置信下限' = c(CLL.f1,CLL.f4,CLL.f5),'置信上限' = c(CUL.f1,CUL.f4,CUL.f5),'区间宽度' = c(UL,UL.f4,UL.f5)) > rownames(Re) = c('一般方法','重要函数f4','重要函数f5') > Re             置信下限  置信上限     区间宽度 一般方法   0.5251122 0.5251122 2.220446e-16 重要函数f4 0.5190214 0.5311288 1.210742e-02 重要函数f5 0.5161340 0.5335964 1.746243e-02