simulate dirichlet process in R

Question

I am reading the paper of "Dirichlet Process Mixtures of Generalized Linear Models" authored by L. A. Hannah. If I would like to simulate the following model

$$\mathcal{P}\sim \text{DP}(c\mathbb{G}_0)$$ $$\theta_i|\mathcal{P}\sim\mathcal{P}$$ $$X_{i,j}|\theta_{i,x}\sim\mathcal{N}(\mu_{ij},\sigma^2_{ij}), j=1,...d$$ $$Y_i|X_i,\theta_{i,y}\sim\mathcal{N}(\beta_{i0}+\sum^d_j=\beta_{ij}X_{ij},\sigma^2_{ij})$$

In R, how can I get $\mathcal{P}$ and $\theta_i|\mathcal{P}$ in $$\mathcal{P}\sim\text{DP}(c\mathbb{G}_0)$$ $$\theta_i|\mathcal{P}\sim\mathcal{P}$$

Answer 1

With certainty, realizations of a Dirichlet Process are probability measures with countable support, as proved by D. Blackwell, The Annals of Statistics 1 (1973), no. 2, 356--358. You can sample realizations from a Dirichlet Process using the constructive stick-breaking representation introduced by J. Sethuraman, Statistica Sinica, 4, 639 (1994). For a Dirichlet process with concentration parameter $c>0$ and centered at some distribution function $\mathbb{G}_0$, you must draw independent random variables $$ B_i\sim \mathrm{Beta}(1,c)\,, $$ and compute $$ P_1=B_1 \, , \qquad P_i=B_i \prod_{j=1}^{i-1}(1-B_j)\, , \qquad i>1 \, , $$ until for some $n\geq 1$ you have $\sum_{i=1}^n P_i\geq 1-\epsilon$, for some $0<\epsilon<1$. Then, drawing independent $Y_i\sim\mathbb{G}_0$, for $i=1,\dots,n$, the (truncated) approximate realization of the Dirichlet Process is the distribution function $$ H(t) = \sum_{i=1}^n P_i\,I_{[Y_i,\infty)}(t) \, . $$ To sample the $\theta_i$'s from this approximate realization of the Dirichlet Process use R's sample with replacement from the $Y_i$'s with probabilities given by the $P_i$'s.

Memory is so cheap nowadays that a more practical way to do the truncation is by taking $n$ "big enough". Here is an example with $c=2$ and $\mathbb{G}_0$ equal to the $\mathrm{N}(0,10)$ distribution function.

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)

Answer 2

Check out the package DPackage in R. It has a lot of functionality for simulating from the Dirichlet Process. Here is a link to the documentation: DPackage. Zen's answer above is pretty good info as well.

Answer 3

Not sure why sample(y, prob = p, replace = TRUE) from zen's answer is necessary.

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)