Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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

share|improve this question
1  
I've entered the formulas in your images as LATEX. Please double-check that I have not inadvertently introduced any errors. –  user777 Apr 25 at 13:45
    
@user777, Thanks! –  user785099 Apr 25 at 14:10

2 Answers 2

up vote 5 down vote accepted

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)
share|improve this answer
    
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. –  user785099 Apr 25 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 at 14:39
    
Thank you very much! –  user785099 Apr 25 at 15:19
    
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$ –  user785099 Apr 25 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 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.