рекурсивно генерировать экспоненциальные случайные величины

Question

Я новичок в рекурсии в R. Я пытаюсь сгенерировать экспоненциальные случайные величины, которые удовлетворяют определенному условию в R. Вот простой пример, показывающий мою попытку генерировать две независимые экспоненциальные случайные величины.

КОД

#### Function to generate two independent exponential random variable that meet two criteria 
gen.exp<-function(q1,q2){ 
    a=rexp(1,q1) # generate exponential random variable 
    b=rexp(1,q2) # generate exponential random variable 
    if((a>10 & a<10.2) & (a+b>15)){ # criteria the random variables must meet 
    return(c(a,b)) 
    }else{ 
    return(gen.exp(q1,q2)) #if above criteria is not met, repeat the process again 
    } 

} 

Пример: q1 = 0,25, q2 = 2

 gen.exp(.25,2) 

Когда я бегу выше код, я получаю следующие ошибки:

1. Ошибка: оценка вложен слишком глубоко: бесконечная рекурсия/параметры (выражения =)?
2. Ошибка во время обертывания: оценка вложен слишком глубоко: бесконечная рекурсия/параметры (выражения =)?

Я уже пробовал модифицировать options(expressions=10000) и дать R увеличить количество итераций. Это, похоже, не помогает в моем случае (возможно, я не использую этот вариант правильно). Я понимаю, что генерирование непрерывных распределений с жесткими критериями, как и выше, может быть, проблема. Сказав это, есть ли способ избежать ошибок? Или, по крайней мере, повторять рекурсию всякий раз, когда возникает ошибка? Рекурсия здесь убита? Существуют ли более простые/лучшие способы генерации желаемых случайных величин? Я ценю любые идеи.

Answer 1

Вы используете простой отказ пробник с два условия

• a > 10 & a < 10.2
• a+ b > 15

Однако, возможность согласования обоих из них является низким, то есть очень медленно. Однако, поскольку вас интересуют экспоненциальные случайные числа, мы можем избежать моделирования чисел, которые мы бы отклонили.

Для генерации экспоненциального случайного числа, мы используем формулу

-rate * log(U) 

, где U представляет собой U(0,1) случайное число.Таким образом, чтобы генерировать значение от экспоненциального распределения больше, чем 10 (скажу), мы просто делаем

-log(U(0, exp(-10*rate))/rate 

или R код

-log(runif(1, 0, exp(-10*rate)))/rate 

Мы можем использовать подобный трюк для верхних границ.

Использование @ функции Роланда сверху, это дает

gen.exp = function(q1, q2, maxiter = 1e3){ 
    i = 0 
    repeat { 
    i = i + 1 
    upper = exp(-10*q1) 
    lower = exp(-10.2*q1) 
    a = -log(runif(1, lower, upper))/q1 
    b = -log(runif(1, 0, exp(-4.8*q2)))/q2 

    if((a>10 & a<10.2) & (a+b>15)) {message(i); return(c(a,b)) 
    if (i > maxiter) stop(paste("Conditions not fulfilled after", maxiter, "tries.")) 
    } 
} 

Обратите внимание, я также распечатаны, сколько итераций потребовалось. Для ваших параметров мне нужно около 2 итераций.

Answer 2

Не использовать рекурсию для этого:

gen.exp<-function(q1, q2, maxiter = 1e3){ 
    i <- 0 
    repeat { 
    i <- i + 1 
    a=rexp(1,q1) # generate exponential random variable 
    b=rexp(1,q2) # generate exponential random variable 
    if((a>10 & a<10.2) & (a+b>15)) return(c(a,b)) # criteria the random variables must meet 
    if (i > maxiter) stop(paste("Conditions not fulfilled after", maxiter, "tries.")) 
    } 
} 

set.seed(42) 
gen.exp(.25, 2) 
#Error in gen.exp(0.25, 2) : Conditions not fulfilled after 1000 tries. 

Вероятность b быть больше, чем 4,8 является:

pexp(4.8, 2, lower.tail = FALSE) 
#[1] 6.772874e-05 

Давайте попробуем с большим количеством итераций:

gen.exp(.25, 2, maxiter = 1e7) 
#[1] 10.08664 5.55414 

Конечно, этот РНГ настолько медленный, что он почти бесполезен. Было бы лучше производить более крупные партии a и b сразу.