We will discuss the Dirichlet distribution, and revisit the Bayesian bootstrap comparing it with the classical bootstrap.

How to simulate from the Dirichlet distribution?¶

Consider the Dirichlet distribution with parameters $a_1, \dots, a_k$. Let us assume that $a_1, \dots, a_k$ are all strictly positive. We will occasionally take some of the $a_j$-s to be small to mimic the effect of them being exactly zero.

The Dirichlet distribution is a distribution on probability vectors $(p_1, \dots, p_k)$ (probability vector means that each $p_j \geq 0$ and $\sum_{j=1}^k p_j = 1$). In other words, if we simulate from the Dirichlet distribution, we obtain probability vectors.

When $k = 3$, we can visualize probability vectors as points in the triangle with vertices $(1, 0, 0)$ (corresponding to $p_1 = 1$), $(0, 1, 0)$ (corresponding to $p_2 = 1$) and $(0, 0, 1)$ (corresponding to $p_3 = 1$). This is done below.

From the above plot (with varying $a_1, a_2, a_3$), the following observations can be made:

1. When $a_i = 1$ for all $i$, the generated probability vector samples are uniformly distributed over the triangle (when $k \geq 3$ the triangle is referred to as the probability simplex).

2. When $a_1 = 1 = a_3$ and $a_2 = 5$ is large, then the probability vectors are piled up near the corner corresponding to $p_2 = 1$.

3. When $a_i = \epsilon$ for a small $\epsilon$ (e.g., $\epsilon = 0.0002$), the points are basically superimposed on the corners of the triangle.

4. When $a_1 = \epsilon$ and $a_2 = a_3 = 1$, then $p_1$ will always be very nearly zero. So the points will be on the line joining the $p_2 = 1$ and $p_3 = 1$ corners. The probability vectors here can be seen as being supported only on the two coordinates 2 and 3.

Now let us how exactly to simulate from the Dirichlet $(a_1, \dots, a_k)$ distribution (in other words, what is the code underlying numpy’s np.random.dirichlet). One can use Gamma distributed random variables to obtain Dirichlet samples. We discuss this below. Let us first recall basic facts about Gamma distributions. Suppose $X_i \overset{\text{ind}}{\sim} \text{Gamma}(\text{shape} = a_i, \text{rate} = 1)$ for $i = 1, \dots, n$. Recall that the Gamma distribution with shape $a$ and rate $\lambda$ has the density:

We will only be using Gamma distributions with rate parameter set to 1 (i.e., $\lambda = 1$). In this case, the Gamma density becomes:

The simplest $\text{Gamma}(a, \lambda)$ density correspnds to $a = 1$ and it is the standard exponential density $\text{Gamma}(1, \lambda) = \text{Exp}(\lambda)$ with density function $\lambda \exp(-\lambda x) I\{x > 0\}$. In particular, $\text{Gamma}(1, 1)$ is $\text{Exp}(1)$ with density $\exp(-x) I\{x > 0\}$.

Now we can state the connection between Dirichlet $(a_1, \dots, a_k)$ and Gamma random variables:

provided

In other words, to generate a probability vector from the Dirichlet $(a_1, \dots, a_k)$ distribution, we simply generate independent Gamma $(a_i, 1)$ distributed random variables, and then normalize them (by dividing each of them with their sum) to get a probability vector.

This Dirichlet-Gamma connection can be used to understand some of the facts about Dirichlet $(a_1, \dots, a_k)$ stated previously. For example, suppose all the $a_i$’s are equal to a small number i.e., $a_i = \epsilon$ for all $i$. Then it turns out that the if we generate $X_i \sim \text{Gamma}(\epsilon, 1), i = 1, \dots, k$, the largest observation will be much larger than all of the others. This can be easily verified by simulation.

This clearly shows that when we normalize the generated Gamma values by their sum, the largest value will become 1 and all the other values will be near zero. Because any of the $k$ Gamma values can be the largest with equal probability, Dirichlet ($\epsilon, \dots, \epsilon$) with $\epsilon \downarrow 0$ will be uniformly distributed on $e_1, \dots, e_k$ where $e_i$ is the probability vector with 1 in the $i$-th position and 0 otherwise.

The Bayesian Bootstrap¶

Coin Tossing¶

The Bayesian bootstrap can be understood with the following very simple example: Suppose a coin is tossed twice resulting in one head and one tail. What can then be said about the coin’s probability $p$ of landing heads?

Bayesian inference with the prior $p \sim \text{Beta}(a, b)$ gives the answer: $p \mid \text{data} \sim \text{Beta}(a + 1, b + 1)$. This answer depends on $a$ and $b$. If we take $a = 0$ and $b = 0$, then we get the posterior $\text{Beta}(1, 1)$, which is same as the uniform distribution on $(0, 1)$. This means that (after seeing the data of 1 head and 1 tail), our uncertainty is spread equally among all the values of $p$. For example, $p \in (0.9, 1)$ has the same posterior probability as $p \in (0, 0.1)$ as well as $(0.45, 0.55)$. This wide uncertainty can be attributed to the fact that we have only observed the results of 2 coin tosses which resulted in 1 heads and 1 tails. If we get more data, then the posterior will get narrower.

An interesting feature of the $\text{Beta}(0, 0)$ prior is the following. Suppose we only observe 1 coin toss which resulted in a heads. Then the posterior of $p$ is $\text{Beta}(1, 0)$ which is concentrated at 1. This means that we are fully confident that $p = 1$.

This might appear quite peculiar because after 1 toss, we have a fully confident concentrated posterior, but after two tosses (which resulted in two different outcomes), we have a very flat posterior.

Jaynes (Subsection 12.4.3) explains this phenomenon using the following analogy. Imagine ordering coffee at a cafe and looking for sugar. You see a jar containing a white granular substance. You taste a single granule and find it sweet. You then you immediately conclude all the granules in the jar are sweet. This corresponds to the $n = 1$ fully concentrated posterior above. But suppose you sample another granule and it turns out to be salty, then you are fully confused and have no idea how much proportion of granules in the jar are sweet and how much proportion are salty. Then your posterior on the proportion of sugar granules in the jar is the uniform distribution on $(0, 1)$.

It is important to specially note that the posterior after 1 head and 1 tail is not concentrated at the observed proportion $1/2$. Rather, it is fully diffuse and uniform over the entire interval $(0, 1)$.

Rolling dice¶

The same phenomenon applies in the Dirichlet case. Suppose we roll a six-faced die six times and observe each face exactly once. What then can we say about the probabilities $p_1, \dots, p_k$ ($k = 6$) associated with the six die faces? If our prior is $(p_1, \dots p_k) \sim \text{Dirichlet}(0, \dots, 0)$, then the posterior is $\text{Dirichlet}(1, \dots, 1)$. This means that we are fully uncertain about the probabilities and have assigned a uniform posterior for all probability vectors of length 6.

Numerical Measurements¶

The same ideas can be applied more generally. Consider the following problem (simpler case of problem discussed in Lectures 12 and 13): A scientist is repeatedly running an experiment which produced the observations $26.6, 38.5, 34.4$ in the first three runs. Based on these, what can we say about the probabilities with which the experiment produces various observations?

The idea is to treat this problem as the die rolling problem. In the die rolling problem, the total number of outcomes is 6. Now, we consider a much larger set of outcomes, say, $\{0, 0.1, 0.2, \dots, 99.8, 99.9, 100\}$ ($k$ is the size of this set). The experiment has been repeated three times resulting in the three distinct outcomes $26.6, 38.5, 34.4$. The goal is to estimate the probabilities $p_1, p_2, \dots, p_k$ where $p_i$ is the probability associated with the $i$-th outcome in the set $\{0, 0.1, 0.2, \dots, 99.8, 99.9, 100\}$.

The prior on $(p_1, \dots, p_k)$ is $\text{Dirichlet}(0, \dots, 0)$ and the posterior will be $\text{Dirichlet}(a_1, \dots, a_k)$ where $a_i$ equals zero for all $i$ except when the outcome is actually observed (in this case, $a_i$ will be 1). So the posterior is supported on probability vectors $(p_1, \dots, p_k)$ for which $p_i \neq 0$ only when $i$ is observed. So the posterior can be viewed as being supported on probability vectors on the three outcome space $\{26.6, 38.5, 34.4\}$. Further, the posterior is uniform on these three-outcome probability vectors:

This posterior refuses to consider any outcome that has not been actually been observed. However, on the observed outcomes, it is displaying maximal posterior uncertainty by allowing all weights with equal probability.

The code below generates posterior samples.

The posterior samples are uniformly distributed over the entire probability triangle.

Comparison to classical bootstrap samples¶

In the classical bootstrap, one obtains re-samples by sampling with replacement from the observed data.

Note that the order of numbers in each bootstrap sample is irrelevant. Each bootstrap sample contains a subset of the three observed numbers $26.6, 38.5, 34.4$. We can obtain an equivalent representation of the bootstrap samples by simply tracking the proportion of times 26.6 appears in the sample, the proportion of times 38.5 appears and the proportion of times 34.4 appears. For example, in this tracking scheme, the bootstrap sample $34.4, 34.4, 26.6$ becomes $1/3, 0, 2/3$. Note that these proportions (because they are nonnegative and sum to one) can also be viewed as probability vectors, so they can be plotted on the probability triangle. This is done below.

Here we see the discrete nature of the classical bootstrap samples. The proportion vectors take only 10 distinct values (even though we generated $N = 1000$ of these proportion vectors). Below we plot both the Bayesian bootstrap posterior probability vectors and the classical bootstrap proportion vectors in the same plot.

For classical bootstrap, we can generate the proportion vectors directly without actually resampling original data. This is because each proportion vector is itself distributed according to the following distribution:

In practice, the goal in these kinds of problems is not to estimate the whole probability distribution of the experiment generating the observations, but instead, to estimate some parameter such as the mean or median $\theta$ of that probability distribution. In each of these methods (classical bootstrap or Bayesian bootstrap), posterior samples for $\theta$ are obtained by calculating the mean or mode for each bootstrap sample.

It is interesting that the posterior of $\theta$ is not uniform over its range, even though the posterior of $(w_1, w_2, w_3)$ was uniform over all probability vectors of length 3. The posterior of $\theta$ is peaked near the mean of the observed data.

Next we plot the sample of $\theta$ obtained from classical bootstrap.

The classical bootstrap histogram for $\theta$ appears very choppy. This is because of its discrete nature and the fact that we have only 3 observations here. It is also not uniform, and is peaked at the mean of the observed data.

When $n$ is large, both Bayesian bootstrap and the classical bootstrap will give similar results.

Ice Cream Problem¶

Suppose you and a friend are planning to open an ice cream shop. You will sell three kinds of ice cream: Vanilla, Chocolate and Coffee. In the first week of opening, you sell

• 19 Vanilla,

• 49 Chocolate, and

• 29 Coffee ice creams.

Problem 1¶

Using a Bayesian model, which proportion of each of the three flavors should you buy for the next week? You may assume that at the beginning of the week, you were unsure which ice cream would sell best (e.g., use a $\textsf{Uniform}$ prior).

Problem 2¶

Suppose you friend wants to make a bet with you:

• If the next customer buys a chocolate ice cream, you will pay your friend $\$20$.

• If the next customer buys a vanilla ice cream, your friend will pay you $\$10$.

• If the next customer buys a coffee ice cream, nothing happens.

Should you take the bet?

Problem 1 Solution¶

For this problem, we use the $\textsf{Dirichlet}-\textsf{Multinomial}$ model with a $\textsf{Uniform}$ prior:

where $p_1,p_2,p_3$ represent the proportion of Vanilla, Chocolate and Coffee flavored ice cream.

We can use the posterior mean as a reasonable estimate for what proportion of ice cream flavor we should order for next week:

Therefore, we should order $2/10, 5/10, 3/10$ proportion of Vanilla, Chocolate and Coffee flavored ice cream for next week.

Problem 2 Solution¶

For this problem, we will consider the Posterior Predictive Distribution.

Applying the law of total probability, we have

In this case, the posterior predictive distribution (of the next draw) is equal to the posterior mean.

Using this result, the expected payoff of the proposed bet is

where

Plugging these probabilities in, we find that the expected payoff is $-\$8$, so we should $\textbf{not}$ take the bet.