为什么我的循环要花很长时间才能完成？

Question

我是运行一个从 2 个指数生成正态分布的算法，如下所示：

set.seed(69420)
j = 1 
Z = c()

# Algorithm
while(j <= 10000){ Y1 <- rexp(1); Y2 <- rexp(1)
    if(Y2 - (Y1 -1)^2/2 >= 0 ){
            X = Y1
            U <- runif(1,0,1)
            if(U > 0.5){
                    Z[j] = X
            } else {
                    Z[j] = -X
            }
            j = j+1
    }
}

然后要求我修改代码如下：“当我们接受样本Y1（如X=Y1）时，随机变量Y=Y2-(Y1-1)2/2也服从指数分布比率为 1，并且它独立于 Y1。修改您的代码以将 Y 作为所需样本之一进行回收。“

我编写了以下代码来完成上述任务。

set.seed(69420)
j = 1
Z = c()
count = 0
Y2 <- rexp(1)
#Algorithm
while(j <= 50){
    Y1 <- rexp(1)
    Y = Y2 - (Y1 -1)^2/2
    if(Y >= 0){
        X = Y1
        Y2 = Y
        U <- runif(1,0,1)
        if(U > 0.5){
            Z[j] = X
        }else {
            Z[j] = -X
        }
        j = j+1
    }
    else if(Y < 0 & j == 1){
       Y2 <- rexp(1)
    }
}

但是，我的循环会一直保持运行，甚至生成 50 次迭代也需要 5 分钟。有什么我做错的导致长时间运行的吗？另外，任何人都可以建议一种对上述代码进行矢量化的方法，以便减少我的处理时间吗？感谢任何帮助。

编辑：在下面发布整个问题以获得更好的解释。

使用拒绝方法，我们可以从标准高斯生成样本使用来自指数分布 Exp(1) 的样本的分布 N(0,1)。这算法如下：

从 Exp(1) 生成 Y1、Y2、独立样本。
如果 Y2 - (Y1 -1)^2 / 2 >= 0，则设置 X = Y1。否则，返回步骤 1。
从统一的 U(0,1) 生成样本 U。如果 U > 0.5，则设置 Z = X。否则，设置 Z = -X。变量 Z 服从高斯分布。

a) 在 R 中实现算法。在单独的文件中提供您的代码，rejection.R。

b) 使用您的代码生成 10000 个高斯分布样本并报告样本均值和标准差。

c) 修改您的代码以计算指数和均匀样本的数量您需要分布以获得 10000 个标准高斯样本。运行代码 10 次并报告所需样本的平均数量。

d) 当我们在步骤2中接受样本Y1时，随机变量Y = Y2 - (Y1 -1)2 / 2 也服从速率为1的指数分布，与Y1无关。修改您的代码以回收 Y 作为步骤 1 中所需的示例之一。在单独的文件中提交您的代码，rejection2.R.

e) 计算您的代码现在使用了多少样本。运行代码 10 次和报告所需的平均样本数。

Answer 1

我认为您对措辞感到困惑。关键是你应该在每次迭代的循环中收获 Y，它应该是一个指数分布的变量，平均值大约为 1:

set.seed(69420)
j = 1 
Z = c()
Y = c()

# Algorithm
while(j <= 10000){ Y1 <- rexp(1); Y2 <- rexp(1)
    if(Y2 - (Y1 -1)^2/2 >= 0 ){
            X = Y1
            Y[j] <- Y2 - (Y1 -1)^2 / 2
            U <- runif(1,0,1)
            if(U > 0.5){
                    Z[j] = X
            } else {
                    Z[j] = -X
            }
            j = j+1
    }
}

所以 Z 应该服从正态分布：

hist(Z)

Y 呈指数分布

hist(Y)

Y 的均值应该接近 1:

mean(Y)
#> [1] 0.9870445

编辑

根据来自 OP 的更多信息，如果 Y 被拒绝，正确的算法只是从指数分布中采样 Y2：

set.seed(69420)
j = 1
Z = c()
count = 0
Y2 <- rexp(1)
#Algorithm
while(j <= 10000){
    Y1 <- rexp(1)
    Y = Y2 - (Y1 -1)^2/2
    if(Y >= 0){
        X = Y1
        Y2 = Y
        U <- runif(1,0,1)
        if(U > 0.5){
            Z[j] = X
        }else {
            Z[j] = -X
        }
        j = j+1
    }
    else {
       Y2 <- rexp(1)
    }
}

mean(Z)
#> [1] -0.00591165

sd(Z)
#> [1] 0.9961794

^{由 reprex package (v2.0.1)}

于 2022-03-27 创建

为什么我的循环要花很长时间才能完成？

Why is my loop taking forever to complete?

performance

loops

r

vectorization