Векторизация кода и застрял, но хорошо

Question

Вот несколько примеров начальных значений для переменных в приведенном ниже коде.

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

Я пытаюсь векторизовать следующий код.Это компонент калькулятора байесовского фактора, который был изначально написан Динесом и адаптирован для R Дэнни Кейем и Томом Багули.Эта часть предназначена для расчета вероятности для теории.У меня эта вещь значительно ускорилась за счет векторизации, но я не могу сопоставить вывод с битом ниже.

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

И ниже - векторизованная версия.точно копирует результаты оригинала, если tails <- 2.Все, что у меня здесь есть, должно просто скопировать и вставить и дать точно такие же результаты.Однако после tails <- 1 вторая функция больше не соответствует точно.Но насколько я могу судить, я делаю эквивалент в новом операторе if тому, что происходит в оригинале.Буду признателен за любую помощь.

(Я попытался создать более минимальный пример, урезая его до цикла «он» и «если операторы и небольшое количество фрагментов», и я просто не смог заставить код выйти из строя.)

Answer 1

Вы отбрасываете наблюдения, где theta==0. Это проблема, потому что вывод dnorm не равен нулю, когда theta==0. Вам нужны эти наблюдения в вашем выводе.

Вместо того, чтобы отбрасывать наблюдения, лучшим решением было бы установить эти элементы на ноль.

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

Answer 2

Исходный расчет имеет ошибку из-за арифметики с плавающей запятой; добавление incr каждый раз приводит к тому, что theta фактически равно 7.204654e-14, когда оно должно равняться нулю. Так что на самом деле это не правильно делает на этом проходе через цикл; он не делает код <=, когда это должно быть. Ваш код (по крайней мере, с этими начальными значениями на моем компьютере).

Ваш код не обязательно гарантирует, что каждый раз поступает правильно; то, что делает seq лучше, чем добавление приращения снова и снова, но это все еще арифметика с плавающей запятой. Вы действительно должны, вероятно, проверять в пределах допуска машины ноль, возможно, используя all.equal или что-то подобное.