> 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