# simulate dirichlet process in R

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}$$

• I've entered the formulas in your images as LATEX. Please double-check that I have not inadvertently introduced any errors.
– Sycorax
Commented Apr 25, 2014 at 13:45

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)

• Hi Zen, thanks for the answer. By the way, I can see how does the parameter $c$ get involved in your code. But I cannot see how does $\mathbb{G}_0$ get involved in your code? This is how I get confused. Commented Apr 25, 2014 at 14:15
• Their $\mathbb{G}_0$ is my $F$. Sorry for changing the notation. I'll fix it. In the R code, I took $\mathbb{G}_0$ as $\mathrm{N}(0,10)$. It can be anything you want. If you still have any doubts, just ask. I wrote a paper with a very short appendix that outlines the main DP results. It may help you to understand the roles of $c$ and $\mathbb{G}_0$. It's open acess. Take a look: scirp.org/journal/PaperDownload.aspx?paperID=44043
– Zen
Commented Apr 25, 2014 at 14:39
• Hi Zen, I found that R also provides a dirichlet function as rdirichlet. Can it be used to get the same result as your stick-breaking code, with $c=2$ and $\mathbb{G}_0$ Commented Apr 25, 2014 at 18:04
• @user785099: In general, no. The Dirichlet Process is a different, more general, animal whose definition depends on the usual Dirichlet distribution. They become the same only when the sampling space is finite (which isn't the case of your problem). You should read carefully the appendix of the paper linked above.
– Zen
Commented Apr 25, 2014 at 18:17

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.

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)