3.1对偶变量法

$统计实验及R语言模拟$

 $第3章方差缩减方法$

代码如下：

#——————————————————————————3.方差缩减方法——————————————————————#
##——————————————————————3.1对偶变量法—————————————————————————##
### [例3.1]
set.seed(1)
n <- 1000
x1 <- runif(n)
mu.e <- mean(x1)
sigma2.e <- var(x1) / n
mu.e
sigma2.e
### 产生与x1完全负相关的随机数
x2 <- 1 - x1
x <- (x1 + x2) / 2
mu.avm <- mean(x)
sigma2.avm <- var(x) / n
mu.avm
sigma2.avm
### 样本均值的方差变为了0

### [例3.2]
set.seed(1)
n <- 1000
u <- runif(n)
x <- (-1 / 2) * log(1 - u)
mu.e <- mean(x)
sigma2.e <- var(x) / n
mu.e
sigma2.e
### 产生与x1负相关的指数分布随机数
x1 <- x
x2 <- (-1 / 2) * log(u)
y <- (x1 +x2) / 2
mu.avm <- mean(y)
sigma2.avm <- var(y) / n
mu.avm
sigma2.avm
sigma2.e / sigma2.avm
cov(x1, x2)
cor(x1, x2)
### 区间估计
alpha <- 0.05
Sg.e <- sqrt((cumsum(x ^ 2) - cumsum(x)) ^ 2 / (1:n)) / ((1:n) - 1)
fre.e <- cumsum(x) / (1:n)
CLL.e <- fre.e - qnorm(1 - alpha / 2) * Sg.e / sqrt(1:n)
CUL.e <- fre.e + qnorm(1 - alpha / 2) * Sg.e / sqrt(1:n)
fre.e[n]
CLL.e[n]
CUL.e[n]
CUL.e[n] - CLL.e[n]
Sg.avm <- sqrt((cumsum(y ^ 2) - (cumsum(y)) ^ 2 / (1:n)) /((1:n) - 1))
fre.avm <- cumsum(y) / (1:n)
CLL.avm <- fre.avm - qnorm(1 - alpha / 2) * Sg.avm / sqrt(1:n)
CUL.avm <- fre.avm + qnorm(1 - alpha / 2) * Sg.avm / sqrt(1:n)
fre.avm[n]
CLL.avm[n]
CUL.avm[n]
CUL.avm[n] - CLL.avm[n]
### 画成图
par(mfcol = c(1, 2))
plot(fre.e,type = 'l', ylim = c(0.2, 0.6), xlab = '模拟量', ylab = '频率', main = '一般方法')
lines(CLL.e, type = 'l', ylim = c(0.2, 0.6), col = 2)
lines(CUL.e, type = 'l', ylim = c(0.2, 0.6), col = 2)
abline(h = 1 / 2, lty = 2, col = 3)
plot(fre.avm, type = 'l', ylim = c(0.2, 0.6), xlab = '模拟量', ylab = '频率', main = '对偶变量法')
lines(CLL.e, type = 'l', ylim = c(0.2, 0.6), col = 2)
lines(CUL.e, type = 'l', ylim = c(0.2, 0.6), col = 2)
abline(h = 1 / 2, lty = 2, col = 3)

### [例3.3]
set.seed(1)
n <- 2000
x <- runif(n)
z <- exp(x)
theta.e <- mean(z)
sigma2.e <- var(z) / n
theta.e
sigma2.e
n1 <- n / 2
x1 <- x[1:n1]
y1 <- exp(x1)
y2 <- exp(1 - x1)
y <- (y1 + y2) / 2
theta.avm <- mean(y)
sigma2.avm <- var(y) / n1
theta.avm
sigma2.avm
sigma2.e / sigma2.avm
(sigma2.e - sigma2.avm) / sigma2.e
cov(y1, y2)
cor(y1, y2)