為什么我的回圈需要永遠完成？

我正在運行一個演算法來從 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 分鐘。我做錯了什么會導致長時間運行嗎？另外，任何人都可以建議一種方法來矢量化上述代碼，以便減少我的處理時間嗎？任何幫助表示贊賞。

編輯：在下面發布整個問題以獲得更好的解釋。

Using the rejection method, we can generate samples from a standard Gaussian distribution N(0,1) using samples from an exponential distribution, Exp(1). The algorithm is as follows:

1. Generate Y1, Y2, independent samples from Exp(1).
2. If Y2 - (Y1 -1)^2 / 2 >= 0, set X = Y1. Otherwise, go back to step 1.
3. Generate a sample U from a uniform U(0,1). If U > 0.5, set Z = X. Otherwise, set Z = -X. The variable Z follows a Gaussian distribution.

a) Implement the algorithm in R. Provide your code in a separate file, rejection.R.

b) Generate 10000 samples of a Gaussian distribution using your code and report the sample mean and standard deviation.

c) Modify your code to count how many samples of the exponential and uniform distribution you need in order to obtain 10000 samples of a standard Gaussian. Run the code 10 times and report the average number of samples required.

d) When we accept the sample Y1 in step 2, the random variable Y = Y2 - (Y1 -1)2 / 2 also follows an exponential distribution with rate 1, and it is independent of Y1. Modify your code to recycle Y as one of the samples required in step 1. Submit your code in a separate file, rejection2.R.

e) Count how many samples is your code using now. Run the code 10 times and report the average number of samples required.

uj5u.com熱心網友回復：

我想你對措辭感到困惑。關鍵是你應該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 提供的更多資訊，正確的演算法是簡單地從指數分布中采樣 Y2（如果Y被拒絕）：

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

^{Created on 2022-03-27 by the reprex package (v2.0.1)}

標籤：r performance loops rstudio vectorization

上一篇：有人可以為我解釋這個回收器視圖配接器，類頭等

下一篇：我被困在匯編中撰寫一個回圈來讀取5個字符并以相反的順序列印它們