[This article was first published on R Works, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't. In 1839, the gifted mathematician Peter Gustav Lejeune Dirichlet was attached to the Philosophy department at the University of Berlin working for less than full pay even though he had become a member of the Prussian Academy of Sciences in 1832. At that time, to become a “full professor” at the university it was required that a candidate deliver a Habilitationsschrift lecture in Latin. Apparently, Dirichlet’s facility with Latin wasn’t up to the task, so like many proficient “adjunct professors” today, Dirichlet took a side gig to support his family. He taught math at a military school. Anyway, I digress. It was about that time that Dirichlet began to work on a problem in celestial mechanics which involved this integral:
Here which is attracted to a point where is the force of attraction and is the Euclidean norm.
After a supernaturally insightful series manipulations detailed by Gupta and Richards, Dirichlet arrived at the following integral
which you will recognize as the Beta function, the normalizing constant for the Dirichlet distribution:
with mean and variance
where
and
The Dirichlet distribution is a multivariate generalization of the Beta distribution that is often used in Bayesian statistics as a prior distribution for categorical and multinomial distributions. I illustrated this use of the Dirichlet in a previous post while constructing a Bayesian model for a three-state Markov chain. The Dirichlet distribution is remarkable in that it brings together 18th and 19th century work in analysis as exemplified by the Gamma, Beta and digamma functions with early 20th ideas from geometry and topology (the simplex) and modern Bayesian statistics.
The (2)-SimplexA simplex is a generalization of the notion of a triangle to multiple dimensions. Informally in K dimensions, a simplex is the simplest polygon that is the convex hull of its K vertices. The vectors that comprise the simplex must be non-negative and sum to 1. So, a simplex is a natural way to represent probabilities that sum to 1 in multidimensional spaces.
The support for the three dimensional Dirichlet distribution, the points on which the distribution is defined, is a (2)-simplex the triangular subset of a 2-dimensional plane intersecting the Euclidean axes at the points (1,0,0), (0,1,0), and (0,0,1). (Orient the triangle in the interactive plot below so that the reference plane is on top and the tip is pointing downward and you will see how the axes line up.)
R packages used in this post
library(ggplot2)library(gganimate)library(dplyr)library(magick)library(MCMCpack) # for rdirichletlibrary(gtools) # for ddirichlet#library(patchwork) # for combining plotslibrary(threejs)library(extraDistr)
Show the code
set.seed(42)# Sample from Dirichlet distribution over 3 categoriesn\_samples <- 2000alpha <- c(1, 1, 1) # uniform prior over the simplexsamples <- rdirichlet(n\_samples, alpha)# 3D coordinates: each row is (x, y, z)x <- samples[,1]y <- samples[,2]z <- samples[,3]# Visualize using threejs scatterplotscatterplot3js(x = x, y = y, z = z, color = "steelblue", size = 0.2, bg = "black", main = "2-simplex", axisLabels = c( "(1,0,0)", "(0,1,0)", "(0,0,1)" ))
When , the Dirichlet density is symmetric about the middle of the simplex, . In the special case when , the density is uniform over the simplex. When all the the density is concentrated at the vertices of the simplex, and when , the density is concentrated in the center of the simplex with most of the mass concentrated on a few values.
How the symetric Dirichlet distribution changes as changesThe following animation, which projects the above plot onto two dimensions, shows how the Dirichlet distribution changes as the common value of , called the concentration parameter, moves systematically from (1,1,1), the uniform distribution, to (0.1,0.1,0.1).
Code for helper functions ```
``` These next two plots, the first and last frames of the animation, clearly show how the density moves from being uniformly distributed over the simplex to being concentrated at the vertices of the simplex. When modeling the development of a multi-state Markov chain, as I was doing in the post I alluded to above, it is common practice to select a uniform Dirichlet prior with . However, if you believe that the process is likely to start off uniformly distributed among the states, then a prior with might be appropriate. If you had reason to believe that the process would start off concentrated on particular states, then you might explore using an asymmetric distribution by setting different values for the . The code driving these animations might be helpful.
Show the code
ggplot(subset(animation\_data, frame == 1), aes(x = x, y = y)) + geom\_point(alpha = 0.3, color = "darkblue") + ggtitle("Initial Frame: Uniform Prior: alpha = (1,1,1)")
Show the code
ggplot(subset(animation\_data, frame == max(animation\_data$frame)), aes(x = x, y = y)) + geom\_point(alpha = 0.3, color = "red") + ggtitle("Final Frame: Concetrated Prior: alpha = (.1,.1,.1)")
This concentration of density as increases is very apparent in this next simulation as moves from (0.1, 0.1, 0.1) to (10.0, 10.0, 10.0). Here we see the distribution concentrating on the mean, = (1/3, 1/3, 1/3).
Show the code ```
``` Note that the animation passes through (0.5, 0.5, 0.5) which is the Jeffreys prior for the Dirichlet distribution.
Variance and Differential EntropyThe Wikipedia article for the Dirichlet distribution prominently displays the distribution’s differential entropy:
where and are defined above and is the digamma function. (This equation triggered my mention of the digamma function above.) But please be advised that differential entropy defined as: for continuous distributions is not the same as the Shannon entropy for discrete distributions and does not have a similar interpretation. Among other things, differential entropy can be negative, in not invariant under a change of variables, and probably doesn’t conform to any intuition you may have developed about maximum entropy. The plot on the left below shows the behavior of the differential entropy for the symmetric Dirichlet distribution we have been considering as moves from (0.1, 0.1, 0.1) to (5.0, 5.0, 5.0). Note that the entropy keeps increasing beyond the point which corresponds to the uniform distribution over the simplex.
Show the code ```
Show the code
k <- seq(.1, 10, by = 0.1)mean_dir <- numeric(length(k))var_dir <- numeric(length(k))for (i in seq_along(k)) { alpha <- k[i] sum_alpha <- alpha * 3 mean_dir[i] <- alpha / sum_alpha var_dir[i] <- (alpha * (sum_alpha - alpha)) / (sum_alpha^2 * (sum_alpha + 1))}df <- data.frame(k = k, mean = mean_dir, variance = var_dir)ggplot(df, aes(x = k, y = variance)) + geom_line(color = "blue", size = 1) + labs( title = "Variance of Dirichlet Distribution vs Alpha Parameter", x = "Alpha Parameter (k)", y = "Variance" )
```
The plot on the right shows the variance of the Dirichlet distribution as a function increasing . We see that the variance decreases towards zero as increases. This is reflected in the second animation above which shows the distribution concentrating on the mean. Uncertainty is going to zero but differential entropy is shooting off towards to infinity. If you are nevertheless intrigued by differential entropy, you may want to have a look at the references I have included below.
References Bela A. Frigyik, Amol Kapila, and Maya R. Gupta, Introduction to the Dirichlet Distribution and Related Processes * Thomas M. Cover and Joy A. Thomas, Elements of Information Theory*, Wiley-Interscience Edition: 2nd Edition, (2006) * Rameshwar Gupta and Donald St. P. Richards, The History of the Dirchet and Liouville Distributions international Statistical Review (2001) * Jiayu Lin, On the Dirichlet Distribution
To leave a comment for the author, please follow the link and comment on their blog: R Works.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.Continue reading: A Note on the Dirichlet Distribution