向量化代碼並卡住但很好

Question

下面是代碼中變量的一些示例起始值。

sd <- 2 
sdtheory <- 1.5 
meanoftheory <- 0.6 
obtained <- 0.8 
tails <- 2 

我正試圖矢量化下面的代碼。它是Bayes因子計算器的一個組成部分，最初由Dienes編寫，並由Danny Kaye改編爲R & Thom Baguley。這部分是爲了計算理論的可能性。我已經通過矢量化大大加快了速度，但是我無法匹配下面的輸出。

area <- 0 
theta <- meanoftheory - 5 * sdtheory 
incr <- sdtheory/200 
for (A in -1000:1000){ 
    theta <- theta + incr 
    dist_theta <- dnorm(theta, meanoftheory, sdtheory) 
    if(identical(tails, 1)){ 
      if (theta <= 0){ 
       dist_theta <- 0 
      } else { 
       dist_theta <- dist_theta * 2 
      } 
     } 
    height <- dist_theta * dnorm(obtained, theta, sd) 
    area <- area + height * incr 
} 
area 

下面是向量化版本。

incr <- sdtheory/200 
newLower <- meanoftheory - 5 * sdtheory + incr 
theta <- seq(newLower, by = incr, length.out = 2001) 
dist_theta <- dnorm(theta, meanoftheory, sdtheory) 
if (tails == 1){ 
    dist_theta <- dist_theta[theta > 0] * 2 
    theta <- theta[theta > 0] 
    } 
height <- dist_theta * dnorm(obtained, theta, sd) 
area <- sum(height * incr) 
area 

此代碼完全複製原始結果，如果tails <- 2。我到目前爲止所得到的所有內容都應該複製並粘貼，並提供完全相同的結果。但是，一旦tails <- 1第二個函數不再精確匹配。但就我所知，我正在做新的if聲明中的等效內容，以處理原始內容。任何幫助，將不勝感激。

（我曾嘗試創建一個更小例子，剝離下來，只是他環路和if語句和切片很少量的，我只是無法讓代碼失敗。）

Answer 1

原始計算由於浮點運算而產生錯誤;每次添加incr導致theta在它應該等於零時實際上等於7.204654e-14。所以在循環中實際上並沒有做正確的事情;當它應該是時，它不會執行<=的代碼。你的代碼是（至少，它是在我的機器上用這些起始值做的）。

您的代碼不一定能保證每次都做正確的事情;什麼seq做比反覆增加一個增量更好，但它仍然是浮點運算。你真的可能應該檢查零機器公差，也許使用all.equal或類似的東西。

Answer 2

你」重新放置觀察點theta==0。這是一個問題，因爲當theta==0時dnorm的輸出不爲零。你需要在你的輸出中看到這些觀察。

與其降低觀察值，更好的解決方案是將這些元素設置爲零。

incr <- sdtheory/200 
newLower <- meanoftheory - 5 * sdtheory + incr 
theta <- seq(newLower, by = incr, length.out = 2001) 
dist_theta <- dnorm(theta, meanoftheory, sdtheory) 
if (tails == 1){ 
    dist_theta <- ifelse(theta < 0, 0, dist_theta) * 2 
    theta[theta < 0] <- 0 
    } 
height <- dist_theta * dnorm(obtained, theta, sd) 
area <- sum(height * incr) 
area