在R中创建新的概率分布（依赖于先前的r.v。）

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：当前r.v。，X_(t-1)：前一个r.v.，其中t = 1,2，...，T，t = 0时的初始值已给出。最后，a和b是两个已知概率。

Answer 1

不确定我是否完全理解，但是您可以执行以下操作。

如果我们说'y'是先前的实现，'x'是当前的实现，那么我们有：

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

因此，我们可以在[0,1]中生成统一变量U，如果y = 0，则在U <= 1-a时将x = 0，否则x = 1；如果y = 1，则如果U <= b，则将x = 0，否则将x = 1

以下函数可以解决问题，其中x0是x的初始值：

rhany <- function(n, a, b, x0 = 0) {
    sim <- c(x0, runif(n-1))
    for (i in 2:n){
        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)
    }
    sim
}

因此，如果您运行该功能：

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