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
Apr 25 '14 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 <- b
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. Apr 25 '14 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
Apr 25 '14 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$ Apr 25 '14 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
Apr 25 '14 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.