在 R 中创建一个新的概率分布（依赖于之前的 rv）答案

【问题标题】：Create a new probability distribution (relying on the previous r.v.) in R在 R 中创建一个新的概率分布（依赖于之前的 rv）
【发布时间】：2020-01-30 21:25:15
【问题描述】：

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

【问题讨论】：

目前尚不清楚 X(t) 如何依赖于 X(t) - 根据您的公式。应该是 X(t) 是 X(t-1) 和 X(t-2) 的函数吗？
不，它是像伯努利这样的 pdf，即 P(x)= p^x 。 (1-p)^(1-x)。所以通常 pdf 依赖于 x，但这个 pdf 中的内容是它也依赖于 x 和之前的 x。
我只是写了 x_t 以使它与前面的 X_(t-1) 不同
我在之前的回答中犯了一个很大的错误，但我现在已经更正了。请看一下。

标签： r function statistics probability-density

【解决方案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 则设置 x = 0 如果 U

以下函数可以解决问题，其中 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循环会减慢函数的速度；在我的机器上生成 1e7 变量大约需要 9 秒。您可以使用 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

【讨论】：