# Understanding the Polya urn model

I am having trouble understanding the plot below (taken from Edwin Chen's blog).

• What is the x-axis supposed to represent? Shouldn't color be a categorical variable? Does the x-axis have to be on the real line for the Polya urn model?
• Are we supposed to see significant changes across runs for the same alpha? • can you link the note instead of the picture? – Elvis Dec 12 '12 at 20:59
• You are allowing a metaphor to confuse you. "Color" is merely a suggestive reference to whatever attributes the balls in the urn might have. If they have a Gaussian distribution, then "colors" are real numbers. – whuber Dec 14 '12 at 19:33
• @whuber I don’t think that’s the issue; I think in the above post, one could replace "categorical" by "discrete". – Elvis Dec 14 '12 at 20:16
• That's how I understand "categorical," @Elvis, but I don't get your point. – whuber Dec 14 '12 at 20:18
• The Pólya urn process generates, even in the long run, a urn containing only a finite number of colors (almost surely). What is represented here can’t be "sample density plots of the colors in the urn"... (Hi, by the way!) – Elvis Dec 14 '12 at 20:36

As I stated in the comments, I think the Pólya urn is used to draw the centers of some normal distributions, and that the plot is a mixture of normal distributions, which seems to make sense as the text you are pointing to is revolving around models for the position of cluster centroids.

Here is a piece of R code that generates similar plots.

polya_urn_model = function(base_color_distribution, num_balls, alpha) {
balls = numeric(num_balls)
for (i in 1:num_balls)
{
balls[i] <- ifelse( runif(1) < alpha / (alpha + i-1),
base_color_distribution(), sample(balls[1:(i-1)], 1))
}
return(balls)
}

my.graph <- function(alpha)
{
N <- 1000;
x <- polya_urn_model(function() rnorm(1), N, alpha)

# the centers of the components
c <- sort(unique(x))

# their weights (proba of being in this component)
w <- as.vector(table(x)/N);

t <- seq(-5,5,length=501)
# computes the density  d(t) =\sum_i w_i f_i(t)
# where f_i = density of N( c_i, σ = 0.4)
d <- rowSums(mapply(function(ce,we) dnorm(t,mean=ce, sd=0.4)*we,c,w))
plot(t, d, type="l")
}

• I understand most of what this code is doing but, as I'm not fluent in R, there's a few lines that I'm unlcear on. Would you mind adding a few comments to explain the gist of what you're doing? Thanks. – jerad Dec 13 '12 at 16:30
• Yep, I have been a litlle cryptic. I'll do that tonight. – Elvis Dec 13 '12 at 17:07
• ✔. Note that the first function is from the blog you’ve linked, I just modified it a little bit, according to my own programming style... – Elvis Dec 13 '12 at 23:26
1. The X-axis is the support of the base distribution $G_0(x)$. In other words if $G_0(x)$ were a normal distribution then the X-axis would be the real number line.

2. Each of those plots corresponds to a random distribution, $G$, drawn from a Dirichlet Process, and smoothed with Gaussian kernels. Specifically, each of the point masses in $G$ is associated with a weighted Gaussian distribution at that location. The plots are of the resulting mixture of Gaussian density. They're all a little different even for the same $\alpha$ because they are randomly sampled distributions.

Note: I think the main source of confusion with these plots is that the X-axis is labeled "Color of ball" yet the plotted line is a continuous function. Perhaps it would have been less confusing if the plots also had colored spikes at the locations of each point mass, to symbolize ball colors. As they are, the plots conceal these centroids.

• Thanks @jerad. What I still don't understand is why the X-axis says Color of the ball and it is a continuous axis. Shouldn't color be a categorical variable? – Amelio Vazquez-Reina Dec 12 '12 at 23:10
• Oh i see. Yeah I could see why that axis label might be confusing. Technically, the X-axis is not "Ball Color" but rather the value of the data point associated with a colored ball. Remember that each ball has a color and data point assigned to it. The color of the ball corresponds to which component distribution its data point was drawn from. – jerad Dec 13 '12 at 1:13
• This is still very unclear to me. Distribution drawn from a Dirichlet process are discrete distributions, just as OP stated. The distributions plotted here are clearly continuous distributions. I suspect that these are the convolutions of such discrete distributions with a Gaussian: this makes sense as this text is centered around the way to draw centroids of clusters. If you think this is not the explanation and these plots are obtained in some other ways, please post a piece of code that would generate these plots. – Elvis Dec 13 '12 at 7:45
• @Elvis, first of all, it seems to me like your comment is in perfect agreement with my answer. In $(1)$ I tried to give a general explanation of what the X-axis is, noting that it's the support of the base distribution $G_0$. It seemed pretty obvious from the blog post that in this case $G_0$ was in fact a Gaussian, meaning that yes, the plots are of Gaussian Mixtures. – jerad Dec 14 '12 at 15:03
• Yes $G$ itself itself is always a discrete distribution, but in the context of these plots the distribution $G$ is smoothed with Gaussian kernels. But I see your objection and I'll edit my answer to clarify this point. – jerad Dec 14 '12 at 15:26