3.2 控制变量法

最新推荐文章于 2023-10-27 17:39:55 发布

wonder1322

最新推荐文章于 2023-10-27 17:39:55 发布

阅读量2.3k

点赞数 4

分类专栏：统计实验及R语言模拟文章标签：概率论

统计实验及R语言模拟专栏收录该内容

6 篇文章 14 订阅

订阅专栏

本文详细介绍了在概率论中如何运用控制变量法来解决复杂的问题。通过实例解析，阐述了控制变量法在概率计算中的应用，帮助读者理解如何在不确定因素中找到关键变量并进行有效控制。

摘要由CSDN通过智能技术生成

在这里插入图片描述

## 【例3.4】
varh <- -0.5 * ((exp(2) - 4 * exp(1) + 3))
varf <- 1 / 12
covhf <- 0.5 * (3 - exp(1))
corhf2 <- (covhf ^ 2) / (varh * varf)
cs <- -covhf / varf
varyc <- varh - covhf ^ 2 / varf
corhf <- covhf / sqrt(varh * varf)
varh
covhf
corhf2
cs
varyc
corhf
# > varh
# [1] 0.242
# > covhf
# [1] 0.1409
# > corhf2
# [1] 0.9837
# > cs
# [1] -1.69
# > varyc
# [1] 0.00394
# > corhf
# [1] 0.9918

set.seed(1)
n <- 1000
x <- runif(n)
theta.e <- mean((exp(x)))
sigma2.e <- sum((exp(x) - theta.e) ^ 2) / (n * (n - 1))
theta.e
sigma2.e
# > theta.e
# [1] 1.717807
# > sigma2.e
# [1] 0.0002441946

set.seed(2)
## 书上对比方案样本容量变成了2000，可能会对结果有影响
m <- 2000
x1 <- runif(m)
f1 <- x1
h1 <- exp(x1)
covhf.m <- cov(h1, f1)
varf.m <- var(f1)
cs.e <- -covhf.m / varf.m
cs.e
# > cs.e
# [1] -1.687069

mu <- 1/2
f <- x
h <- exp(x)
y <- h + cs.e * (f - mu)
theta.cv <- mean(y)
sigma2.cv <- sum((y - theta.cv) ^ 2) / (n * (n - 1))
theta.cv
sigma2.cv
(sigma2.e - sigma2.cv) / sigma2.e
cor(h, f)
cor(h, f) ^ 2
# > theta.cv
# [1] 1.71833
# > sigma2.cv
# [1] 4.00642e-06
# > (sigma2.e - sigma2.cv) / sigma2.e
# [1] 0.9835933
# > cor(h, f)
# [1] 0.9917653
# > cor(h, f) ^ 2
# [1] 0.9835983

set.seed(2)
## 将对比方案样本容量变成了1000
m <- 1000
x1 <- runif(m)
f1 <- x1
h1 <- exp(x1)
covhf.m <- cov(h1, f1)
varf.m <- var(f1)
cs.e <- -covhf.m / varf.m
cs.e
# > cs.e
# [1] -1.690381

mu <- 1/2
f <- x
h <- exp(x)
y <- h + cs.e * (f - mu)
theta.cv <- mean(y)
sigma2.cv <- sum((y - theta.cv) ^ 2) / (n * (n - 1))
theta.cv
sigma2.cv
(sigma2.e - sigma2.cv) / sigma2.e
cor(h, f)
cor(h, f) ^ 2
# > theta.cv
# [1] 1.718328
# > sigma2.cv
# [1] 4.011995e-06
# > (sigma2.e - sigma2.cv) / sigma2.e
# [1] 0.9835705
# > cor(h, f)
# [1] 0.9917653
# > cor(h, f) ^ 2
# [1] 0.9835983
## 样本容量为1000时，与2000的差别很小
## 但实际上样本容量不一致，对比的意义不大，都设为1000为好

在这里插入图片描述


## 【例3.5】
n <- 1000
set.seed(1)
U <- runif(n)
set.seed(2)
V <- runif(n)
cor(U, V)
h <- function(U, V){
  exp((U + V) ^ 2)
}
theta.e <- mean(h(U, V))
sigma2.e <- sum((h(U, V) - theta.e) ^ 2) / (n * (n - 1))
theta.e
sigma2.e
# > theta.e
# [1] 4.978106
# > sigma2.e
# [1] 0.039346

m <- 2000
set.seed(3)
U1 <- runif(m)
set.seed(4)
V1 <- runif(m)
cor(U1, V1)
f1 <- function(U, V){
  U + V
}
f2 <- function(U, V){
  (U + V) ^ 2
}
f3 <- function(U, V){
  exp(U + V)
}
covhf1 <- cov(h(U1, V1), f1(U1, V1))
covhf2 <- cov(h(U1, V1), f2(U1, V1))
covhf3 <- cov(h(U1, V1), f3(U1, V1))
varf1 <- var(f1(U1, V1))
varf2 <- var(f2(U1, V1))
varf3 <- var(f3(U1, V1))
cs1 <- -covhf1 / varf1
cs2 <- -covhf2 / varf2
cs3 <- -covhf3 / varf3
c('cs2' = cs1, 'cs2' = cs2, 'cs3' = cs3)
# > c('cs2' = cs1, 'cs2' = cs2, 'cs3' = cs3)
# cs2        cs2        cs3 
# -10.970274  -6.096927  -4.240745

mu1 <- 1
mu2 <- 7/6
mu3 <- (exp(1) - 1) ^ 2
y1 <- h(U, V) + cs1 * (f1(U, V) - mu1)
y2 <- h(U, V) + cs2 * (f2(U, V) - mu2)
y3 <- h(U, V) + cs3 * (f3(U, V) - mu3)
theta1.cv <- mean(y1)
theta2.cv <- mean(y2)
theta3.cv <- mean(y3)
sigma2.cv1 <- sum((y1 - theta1.cv) ^ 2) / (n * (n - 1))
sigma2.cv2 <- sum((y2 - theta2.cv) ^ 2) / (n * (n - 1))
sigma2.cv3 <- sum((y3 - theta3.cv) ^ 2) / (n * (n - 1))
c('theta1.cv' = theta1.cv, 'theta2.cv' = theta2.cv, 'theta3.cs' = theta3.cv)
# theta1.cv  theta2.cv  theta3.cs 
# 5.002781   4.994224   4.988433 

c('sigma2.cv1' = sigma2.cv1, 'sigma2.cv2' = sigma2.cv2, 'sigma2.cv3' = sigma2.cv3)
# sigma2.cv1  sigma2.cv2  sigma2.cv3 
# 0.016383028 0.009319562 0.008609911 

vrr1 <- (sigma2.e - sigma2.cv1) / sigma2.e
vrr2 <- (sigma2.e - sigma2.cv2) / sigma2.e
vrr3 <- (sigma2.e - sigma2.cv3) / sigma2.e
c('vrr1' = vrr1, 'vrr2' = vrr2, 'vrr3' = vrr3)
# vrr1      vrr2      vrr3 
# 0.5836164 0.7631383 0.7811744 

cor1 <- cor(h(U, V), f1(U, V))
cor2 <- cor(h(U, V), f2(U, V))
cor3 <- cor(h(U, V), f3(U, V))
c('cor1' = cor1, 'cor2' = cor2, 'cor3' = cor3)
# cor1      cor2      cor3 
# 0.7653984 0.8748768 0.8851338 

alpha <- 0.05
CLL.h <- theta.e - qnorm(1 - alpha / 2) * sqrt(sigma2.e)
CUL.h <- theta.e + qnorm(1 - alpha / 2) * sqrt(sigma2.e)
CLL.f1 <- theta1.cv - qnorm(1 - alpha / 2) * sqrt(sigma2.cv1)
CUL.f1 <- theta1.cv + qnorm(1 - alpha / 2) * sqrt(sigma2.cv1)
CLL.f2 <- theta2.cv - qnorm(1 - alpha / 2) * sqrt(sigma2.cv2)
CUL.f2 <- theta2.cv + qnorm(1 - alpha / 2) * sqrt(sigma2.cv2)
CLL.f3 <- theta3.cv - qnorm(1 - alpha / 2) * sqrt(sigma2.cv3)
CUL.f3 <- theta3.cv + qnorm(1 - alpha / 2) * sqrt(sigma2.cv3)
UL <- CUL.h - CLL.h
UL1 <- CUL.f1 - CLL.f1
UL2 <- CUL.f2 - CLL.f2
UL3 <- CUL.f3 - CLL.f3
Re <- data.frame('置信下限' = c(CLL.h, CLL.f1, CLL.f2, CLL.f3),
                 '置信上限' = c(CUL.h, CUL.f1, CUL.f2, CUL.f3),
                 '区间宽带' = c(UL, UL1, UL2, UL3))
rownames(Re) <- c('h', 'f1', 'f2', 'f3')
options(digits = 4)
Re
#      置信下限 置信上限 区间宽带
# h     4.589    5.367   0.7776
# f1    4.752    5.254   0.5017
# f2    4.805    5.183   0.3784
# f3    4.807    5.170   0.3637