机器算法验证 - 在 R 中模拟狄利克雷过程 - 吾爱随笔录

在 R 中模拟狄利克雷过程

机器算法验证 r 贝叶斯马尔可夫链蒙特卡罗狄利克雷分布狄利克雷过程

2022-04-13 13:43:22

我正在阅读 LA Hannah 撰写的“广义线性模型的狄利克雷过程混合”论文。如果我想模拟以下模型

P \sim DP (c G_{0})

$\mathcal{P}\sim \text{DP}(c\mathbb{G}_0)$

θ_{i} | P \sim P

$\theta_i|\mathcal{P}\sim\mathcal{P}$

X_{i, j} | θ_{i, x} \sim N (μ_{i j}, σ_{i j}^{2}), j = 1, . . . d

$X_{i,j}|\theta_{i,x}\sim\mathcal{N}(\mu_{ij},\sigma^2_{ij}), j=1,...d$

Y_{i} | X_{i}, θ_{i, y} \sim N (β_{i 0} + \sum_{j}^{d} = β_{i j} X_{i j}, σ_{i j}^{2})

$Y_i|X_i,\theta_{i,y}\sim\mathcal{N}(\beta_{i0}+\sum^d_j=\beta_{ij}X_{ij},\sigma^2_{ij})$

在 R 中，我怎样才能得到和在 $\mathcal{P}$ $\theta_i|\mathcal{P}$

P \sim DP (c G_{0})

$\mathcal{P}\sim\text{DP}(c\mathbb{G}_0)$

θ_{i} | P \sim P

$\theta_i|\mathcal{P}\sim\mathcal{P}$

3个回答

可以肯定的是，狄利克雷过程的实现是具有可数支持的概率度量，正如D. Blackwell 在统计年鉴1（1973 年）第 1 期所证明的那样。2、356--358。您可以使用 J. Sethuraman 在Statistica Sinica , 4 , 639 (1994) 中介绍的构造性断棒表示从狄利克雷过程中采样实现。对于浓度参数且以某个分布函数为中心的 Dirichlet 过程，您必须绘制独立随机变量并计算 $c>0$ $\mathbb{G}_0$

B_{i} \sim B e t a (1, c),

$B_i\sim \mathrm{Beta}(1,c)\,,$

P_{1} = B_{1}, P_{i} = B_{i} \prod_{j = 1}^{i - 1} (1 - B_{j}), i > 1,

$P_1=B_1 \, , \qquad P_i=B_i \prod_{j=1}^{i-1}(1-B_j)\, , \qquad i>1 \, ,$ 直到对于一些你有，对于一些。然后，绘制独立的，对于，狄利克雷过程的（截断）近似实现是分布函数为了从狄利克雷过程的这种近似实现中对进行采样，请使用's给出的概率替换。

n \geq 1

$n\geq 1$

\sum_{i = 1}^{n} P_{i} \geq 1 - ϵ

$\sum_{i=1}^n P_i\geq 1-\epsilon$

0 < ϵ < 1

$0<\epsilon<1$

Y_{i} \sim G_{0}

$Y_i\sim\mathbb{G}_0$

i = 1, \dots, n

$i=1,\dots,n$

H (t) = \sum_{i = 1}^{n} P_{i} I_{[Y_{i}, \infty)} (t) .

$H(t) = \sum_{i=1}^n P_i\,I_{[Y_i,\infty)}(t) \, .$

θ_{i}

$\theta_i$ Rsample

Y_{i}

$Y_i$

P_{i}

$P_i$

现在内存太便宜了，一个更实用的截断方法是取 “足够大”。这是一个且等于分布函数的示例。 $n$ $c=2$ $\mathbb{G}_0$ $\mathrm{N}(0,10)$

c <- 2
G_0 <- function(n) rnorm(n, 0, 10)
n <- 100
b <- rbeta(n, 1, c)
p <- numeric(n)
p[1] <- b[1]
p[2:n] <- sapply(2:n, function(i) b[i] * prod(1 - b[1:(i-1)]))
y <- G_0(n)
theta <- sample(y, prob = p, replace = TRUE)

查看DPackageR 中的包。它具有许多用于从狄利克雷过程进行模拟的功能。这是文档的链接：DPackage。Zen 上面的回答也是很好的信息。

不知道为什么sample(y, prob = p, replace = TRUE)禅宗的回答是必要的。

library(tidyverse)

##concentration parameter
c <- 1000
##base distribution
G_0 <- function(n) rnorm(n, 0, 1)
##finite approximate realization of Dirichlet Process
n <- 1000
b <- rbeta(n, 1, c)
p <- numeric(n)
p[1] <- b[1]
p[2:n] <- sapply(2:n, function(i) b[i] * prod(1 - b[1:(i-1)]))
##check summation of p must be 1
sum(p)
##P(theta_i)=p_i where theta follows i.i.d G_0
theta <- G_0(n)
##plot is similar to https://en.wikipedia.org/wiki/File:Dirichlet_process_draws.svg
df1 <- data.frame(theta = theta, p = p)
df1 %>%
    ggplot(aes(x = theta , y = p)) +
    geom_col(color = "black") +
    xlim(-4,4)

其它你可能感兴趣的问题

上一篇python中加权最小二乘的多元回归？下一篇如何为正则化回归选择调整参数进行解释？