非參貝葉斯

以前研究過一段時間的非參貝葉斯，可是對爲何叫「非參」，以及dirichlet process不是很瞭解，今天看到一篇神文，深刻淺出的娓娓道來

爲何叫「非參」：傳統的聚類在開始的時候就要設定類別的數目，而「非參」是指隨着數據的不斷增長，新的類別不斷加入

這篇神文更是簡歷了 餐廳問題和 dirichlet process 的關係

轉自：http://blog.echen.me/

Imagine you’re a budding chef. A data-curious one, of course, so you start by taking a set of foods (pizza, salad, spaghetti, etc.) and ask 10 friends how much of each they ate in the past day.

Your goal: to find natural groups of foodies, so that you can better cater to each cluster’s tastes. For example, your fratboy friends might love wings and beer, your anime friends might love soba and sushi, your hipster friends probably dig tofu, and so on.

So how can you use the data you’ve gathered to discover different kinds of groups?

One way is to use a standard clustering algorithm like k-means or Gaussian mixture modeling(see this previous post for a brief introduction). The problem is that these both assume a fixednumber of clusters, which they need to be told to find. There are a couple methods for selecting the number of clusters to learn (e.g., the gap and prediction strength statistics), but the problem is a more fundamental one: most real-world data simply doesn’t have a fixed number of clusters.

That is, suppose we’ve asked 10 of our friends what they ate in the past day, and we want to find groups of eating preferences. There’s really an infinite number of foodie types (carnivore, vegan, snacker, Italian, healthy, fast food, heavy eaters, light eaters, and so on), but with only 10 friends, we simply don’t have enough data to detect them all. (Indeed, we’re limited to 10 clusters!) So whereas k-means starts with the incorrect assumption that there’s a fixed, finite number of clusters that our points come from, no matter if we feed it more data, what we’d really like is a method positing an infinite number of hidden clusters that naturally arise as we ask more friends about their food habits. (For example, with only 2 data points, we might not be able to tell the difference between vegans and vegetarians, but with 200 data points, we probably could.)

Luckily for us, this is precisely the purview of nonparametric Bayes.*

*Nonparametric Bayes refers to a class of techniques that allow some parameters to change with the data. In our case, for example, instead of fixing the number of clusters to be discovered, we allow it to grow as more data comes in.

A Generative Story

Let’s describe a generative model for finding clusters in any set of data. We assume an infinite set of latent groups, where each group is described by some set of parameters. For example, each group could be a Gaussian with a specified mean  μi  and standard deviation  σi , and these group parameters themselves are assumed to come from some base distribution  G0 . Data is then generated in the following manner:

• Select a cluster.
• Sample from that cluster to generate a new point.

(Note the resemblance to a finite mixture model.)

For example, suppose we ask 10 friends how many calories of pizza, salad, and rice they ate yesterday. Our groups could be:

• A Gaussian centered at (pizza = 5000, salad = 100, rice = 500) (i.e., a pizza lovers group).
• A Gaussian centered at (pizza = 100, salad = 3000, rice = 1000) (maybe a vegan group).
• A Gaussian centered at (pizza = 100, salad = 100, rice = 10000) (definitely Asian).
• …

When deciding what to eat when she woke up yesterday, Alice could have thought girl, I’m in the mood for pizza and her food consumption yesterday would have been a sample from the pizza Gaussian. Similarly, Bob could have spent the day in Chinatown, thereby sampling from the Asian Gaussian for his day’s meals. And so on.

The big question, then, is: how do we assign each friend to a group?

Assigning Groups

Chinese Restaurant Process

One way to assign friends to groups is to use a Chinese Restaurant Process. This works as follows: Imagine a restaurant where all your friends went to eat yesterday…

• Initially the restaurant is empty.
• The first person to enter (Alice) sits down at a table (selects a group). She then orders food for the table (i.e., she selects parameters for the group); everyone else who joins the table will then be limited to eating from the food she ordered.
• The second person to enter (Bob) sits down at a table. Which table does he sit at? With probability  α/(1+α)  he sits down at a new table (i.e., selects a new group) and orders food for the table; with probability  1/(1+α)  he sits with Alice and eats from the food she’s already ordered (i.e., he’s in the same group as Alice).
• …
• The (n+1)-st person sits down at a new table with probability  α/(n+α) , and at table k with probability  nk/(n+α) , where  nk  is the number of people currently sitting at table k.

Note a couple things:

• The more people (data points) there are at a table (cluster), the more likely it is that people (new data points) will join it. In other words, our groups satisfy a rich get richer property.
• There’s always a small probability that someone joins an entirely new table (i.e., a new group is formed).
• The probability of a new group depends on  α . So we can think of  α  as a dispersion parameter that affects the dispersion of our datapoints. The lower alpha is, the more tightly clustered our data points; the higher it is, the more clusters we have in any finite set of points.

(Also notice the resemblance between table selection probabilities and a Dirichlet distribution…)

Just to summarize, given n data points, the Chinese Restaurant Process specifies a distribution over partitions (table assignments) of these points. We can also generate parameters for each partition/table from a base distribution  G0  (for example, each table could represent a Gaussian whose mean and standard deviation are sampled from  G0 ), though to be clear, this is not part of the CRP itself.

Code

Since code makes everything better, here’s some Ruby to simulate a CRP:

Chinese Restaurant Process chinese_restaurant_process.rb

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

# Generate table assignments for `num_customers` customers, according to # a Chinese Restaurant Process with dispersion parameter `alpha`. # # returns an array of integer table assignments def chinese_restaurant_process(num_customers, alpha) return [] if num_customers <= 0 table_assignments = [1] # first customer sits at table 1 next_open_table = 2 # index of the next empty table # Now generate table assignments for the rest of the customers. 1.upto(num_customers - 1) do |i| if rand < alpha.to_f / (alpha + i) # Customer sits at new table. table_assignments << next_open_table next_open_table += 1 else # Customer sits at an existing table. # He chooses which table to sit at by giving equal weight to each # customer already sitting at a table. which_table = table_assignments[rand(table_assignments.size)] table_assignments << which_table end end table_assignments end

And here’s some sample output:

Chinese Restaurant Process chinese_restaurant_process.rb

1 2 3 4 5 6 7 8 9 10 11 12 13 14

> chinese_restaurant_process(num_customers = 10, alpha = 1) 1, 2, 3, 4, 3, 3, 2, 1, 4, 3 # table assignments from run 1 1, 1, 1, 1, 1, 1, 2, 2, 1, 3 # table assignments from run 2 1, 2, 2, 1, 3, 3, 2, 1, 3, 4 # table assignments from run 3 > chinese_restaurant_process(num_customers = 10, alpha = 3) 1, 2, 1, 1, 3, 1, 2, 3, 4, 5 1, 2, 3, 3, 4, 3, 4, 4, 5, 5 1, 1, 2, 3, 1, 4, 4, 3, 1, 1 > chinese_restaurant_process(num_customers = 10, alpha = 5) 1, 2, 1, 3, 4, 5, 6, 7, 1, 8 1, 2, 3, 3, 4, 5, 6, 5, 6, 7 1, 2, 3, 4, 5, 6, 2, 7, 2, 1

Notice that as we increase  α , so too does the number of distinct tables increase.

Polya Urn Model

Another method for assigning friends to groups is to follow the Polya Urn Model. This is basically the same model as the Chinese Restaurant Process, just with a different metaphor.

• We start with an urn containing  αG0(x)  balls of 「color」  x , for each possible value of  x . ( G0  is our base distribution, and  G0(x)  is the probability of sampling  x  from  G0 ). Note that these are possibly fractional balls.
• At each time step, draw a ball from the urn, note its color, and then drop both the original ball plus a new ball of the same color back into the urn.

Note the connection between this process and the CRP: balls correspond to people (i.e., data points), colors correspond to table assignments (i.e., clusters), alpha is again a dispersion parameter (put differently, a prior), colors satisfy a rich-get-richer property (since colors with many balls are more likely to get drawn), and so on. (Again, there’s also a connection between this urn model and the urn model for the (finite) Dirichlet distribution…)

To be precise, the difference between the CRP and the Polya Urn Model is that the CRP specifies only a distribution over partitions (i.e., table assignments), but doesn’t assign parameters to each group, whereas the Polya Urn Model does both.

Code

Again, here’s some code for simulating a Polya Urn Model:

Polya Urn Model polya_urn_model.rb

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

# Draw `num_balls` colored balls according to a Polya Urn Model # with a specified base color distribution and dispersion parameter # `alpha`. # # returns an array of ball colors def polya_urn_model(base_color_distribution, num_balls, alpha) return [] if num_balls <= 0 balls_in_urn = [] 0.upto(num_balls - 1) do |i| if rand < alpha.to_f / (alpha + balls_in_urn.size) # Draw a new color, put a ball of this color in the urn. new_color = base_color_distribution.call balls_in_urn << new_color else # Draw a ball from the urn, add another ball of the same color. ball = balls_in_urn[rand(balls_in_urn.size)] balls_in_urn << ball end end balls_in_urn end

And here’s some sample output, using a uniform distribution over the unit interval as the color distribution to sample from:

Polya Urn Model polya_urn_model.rb

1 2 3 4 5 6

> unit_uniform = lambda { (rand * 100).to_i / 100.0 } > polya_urn_model(unit_uniform, num_balls = 10, alpha = 1) 0.27, 0.89, 0.89, 0.89, 0.73, 0.98, 0.43, 0.98, 0.89, 0.53 # colors in the urn from run 1 0.26, 0.26, 0.46, 0.26, 0.26, 0.26, 0.26, 0.26, 0.26, 0.85 # colors in the urn from run 2 0.96, 0.87, 0.96, 0.87, 0.96, 0.96, 0.87, 0.96, 0.96, 0.96 # colors in the urn from run 3

Code, Take 2

Here’s the same code for a Polya Urn Model, but in R:

Polya Urn Model polya_urn_model.R

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

# Return a vector of `num_balls` ball colors according to a Polya Urn Model # with dispersion `alpha`, sampling from a specified base color distribution. polya_urn_model = function(base_color_distribution, num_balls, alpha) { balls = c() for (i in 1:num_balls) { if (runif(1) < alpha / (alpha + length(balls))) { # Add a new ball color. new_color = base_color_distribution() balls = c(balls, new_color) } else { # Pick out a ball from the urn, and add back a # ball of the same color. ball = balls[sample(1:length(balls), 1)] balls = c(balls, ball) } } balls }

Here are some sample density plots of the colors in the urn, when using a unit normal as the base color distribution:

Notice that as alpha increases (i.e., we sample more new ball colors from our base; i.e., as we place more weight on our prior), the colors in the urn tend to a unit normal (our base color distribution).

And here are some sample plots of points generated by the urn, for varying values of alpha:

• Each color in the urn is sampled from a uniform distribution over [0,10]x[0,10] (i.e., a [0, 10] square).
• Each group is a Gaussian with standard deviation 0.1 and mean equal to its associated color, and these Gaussian groups generate points.

Notice that the points clump together in fewer clusters for low values of alpha, but become more dispersed as alpha increases.

Stick-Breaking Process

Imagine running either the Chinese Restaurant Process or the Polya Urn Model without stop. For each group  i , this gives a proportion  wi  of points that fall into group  i .

So instead of running the CRP or Polya Urn model to figure out these proportions, can we simply generate them directly?

This is exactly what the Stick-Breaking Process does:

• Start with a stick of length one.
• Generate a random variable  β1∼Beta(1,α) . By the definition of the Beta distribution, this will be a real number between 0 and 1, with expected value  1/(1+α) . Break off the stick at  β1 ;  w1  is then the length of the stick on the left.
• Now take the stick to the right, and generate  β2∼Beta(1,α) . Break off the stick  β2  into the stick. Again,  w2  is the length of the stick to the left, i.e.,  w2=(1−β1)β2 .
• And so on.

Thus, the Stick-Breaking process is simply the CRP or Polya Urn Model from a different point of view. For example, assigning customers to table 1 according to the Chinese Restaurant Process is equivalent to assigning customers to table 1 with probability  w1 .

Code

Here’s some R code for simulating a Stick-Breaking process:

Stick-Breaking Process stick_breaking_process.R

1 2 3 4 5 6 7 8 9 10 11 12 13

# Return a vector of weights drawn from a stick-breaking process # with dispersion `alpha`. # # Recall that the kth weight is # \beta_k = (1 - \beta_1) * (1 - \beta_2) * ... * (1 - \beta_{k-1}) * beta_k # where each $\\beta\_i$ is drawn from a Beta distribution # \beta_i ~ Beta(1, \alpha) stick_breaking_process = function(num_weights, alpha) { betas = rbeta(num_weights, 1, alpha) remaining_stick_lengths = c(1, cumprod(1 - betas))[1:num_weights] weights = remaining_stick_lengths * betas weights }

And here’s some sample output:

Notice that for low values of alpha, the stick weights are concentrated on the first few weights (meaning our data points are concentrated on a few clusters), while the weights become more evenly dispersed as we increase alpha (meaning we posit more clusters in our data points).

Dirichlet Process

Suppose we run a Polya Urn Model several times, where we sample colors from a base distribution  G0 . Each run produces a distribution of colors in the urn (say, 5% blue balls, 3% red balls, 2% pink balls, etc.), and the distribution will be different each time (for example, 5% blue balls in run 1, but 1% blue balls in run 2).

For example, let’s look again at the plots from above, where I generated samples from a Polya Urn Model with the standard unit normal as the base distribution:

Each run of the Polya Urn Model produces a slighly different distribution, though each is 「centered」 in some fashion around the standard Gaussian I used as base. In other words, the Polya Urn Model gives us a distribution over distributions (we get a distribution of ball colors, and this distribution of colors changes each time) – and so we finally get to the Dirichlet Process.

Formally, given a base distribution  G0  and a dispersion parameter  α , a sample from the Dirichlet Process  DP(G0,α)  is a distribution  G∼DP(G0,α) . This sample  G  can be thought of as a distribution of colors in a single simulation of the Polya Urn Model; sampling from  G  gives us the balls in the urn.

So here’s the connection between the Chinese Restaurant Process, the Polya Urn Model, the Stick-Breaking Process, and the Dirichlet Process:

• Dirichlet Process: Suppose we want samples  xi∼G , where  G  is a distribution sampled from the Dirichlet Process  G∼DP(G0,α) .
• Polya Urn Model: One way to generate these values  xi  would be to take a Polya Urn Model with color distribution  G0  and dispersion  α . ( xi  would be the color of the ith ball in the urn.)
• Chinese Restaurant Process: Another way to generate  xi  would be to first assign tables to customers according to a Chinese Restaurant Process with dispersion  α . Every customer at the nth table would then be given the same value (color) sampled from  G0 . ( xi  would be the value given to the ith customer;  xi  can also be thought of as the food at table  i , or as the parameters of table  i .)
• Stick-Breaking Process: Finally, we could generate weights  wk  according to a Stick-Breaking Process with dispersion  α . Next, we would give each weight  wk  a value (or color)  vk  sampled from  G0 . Finally, we would assign  xi  to value (color)  vk  with probability  wk .

Recap

Let’s summarize what we’ve discussed so far.

We have a bunch of data points  pi  that we want to cluster, and we’ve described four essentially equivalent generative models that allow us to describe how each cluster and point could have arisen.

In the Chinese Restaurant Process:

• We generate table assignments  g1,…,gn∼CRP(α)  according to a Chinese Restaurant Process. ( gi  is the table assigned to datapoint  i .)
• We generate table parameters  ϕ1,…,ϕn∼G0  according to the base distribution  G0 , where  ϕk  is the parameter for the kth distinct group.
• Given table assignments and table parameters, we generate each datapoint  pi∼F(ϕgi)  from a distribution  F  with the specified table parameters. (For example,  F  could be a Gaussian, and ϕi  could be a parameter vector specifying the mean and standard deviation).

In the Polya Urn Model:

• We generate colors  ϕ1,…,ϕn∼Polya(G0,α)  according to a Polya Urn Model. ( ϕi  is the color of the ith ball.)
• Given ball colors, we generate each datapoint  pi∼F(ϕi) .

In the Stick-Breaking Process:

• We generate group probabilities (stick lengths)  w1,…,w∞∼Stick(α)  according to a Stick-Breaking process.
• We generate group parameters  ϕ