Создайте новое распределение вероятностей (опираясь на предыдущий rv) в R

Question

Я хочу написать этот pdf и сгенерировать случайные числа, используя его.

Предположим, что X - это случайная переменная, которая принимает значения только 0 или 1, pdf выглядит следующим образом:

P(X_t) = a^[(1-X_(t-1))*X_t] * (1-a)^[(1-X_(t-1))*(1-X_t)] * b^[X_(t-1)*(1-X_t)] * (1-b)^[X_(t-1)*X_t]

X_t: текущий rv, X_(t-1): предыдущий rv, где t = 1,2, ..., T и задано начальное значение при t = 0. Наконец, a и b - две известные вероятности.

Answer 1

Не уверен, что я полностью понял, но вы могли бы поступить следующим образом.

Если мы говорим «у» - это предыдущая реализация, а «х» - текущая, то мы имеем:

P(x=0|y=0) = 1-a
P(x=1|y=0) = a
P(x=0|y=1) = b
P(x=1|y=1) = 1-b

Итак, мы можем сгенерировать равномерные переменные U в [0,1], и если y = 0, тогда установить x = 0, если U <= 1-a, иначе x = 1;и если y = 1, тогда установите x = 0, если U <= b, иначе установите x = 1 </p>

. Следующая функция может помочь, где x0 - начальное значение x:

* 1009. *

Итак, если вы запустите функцию:

rhany(10, 0.1, 0.7)
[1] 0 1 0 1 1 1 0 1 1 0

По общему признанию, цикл for замедляет функцию;На моей машине требуется около 9 секунд, чтобы сгенерировать 1e7 вариаций. Вы можете переопределить, используя пакет Rcpp:

library(Rcpp)

cppFunction('NumericVector rhanya(double a, double b, NumericVector zs) {
    int n = zs.size();
    NumericVector sim = zs;
    for (int i = 1; i < n; i++) {
        sim[i] = (sim[i-1] == 0) * ((sim[i] <= 1-a) * 0 + (sim[i] > a) * 1) + (sim[i-1] == 1) * ((sim[i] <= b) * 0 + (sim[i] > 1-b) * 1);
    }
    return(sim);
}')

rhany1 <- function(n, a, b, x0 = 0) {
    sim <- c(x0, runif(n-1))
    rhanya(a, b, sim)
}

Эта функция rhany1 занимает менее 0,5 секунды для генерации 1e7 вариаций.

Вы можете проверить, что две функции rhany и rhany1 дают одинаковые результаты, когда установлено одинаковое начальное число:

set.seed(123); bc <- rhany(1e7, 0.1, 0.7)
set.seed(123); ac <- rhany1(1e7, 0.1, 0.7)
all.equal(ac, bc)
[1] TRUE