Skip to content

Instantly share code, notes, and snippets.

Last active October 25, 2023 11:07
  • Star 41 You must be signed in to star a gist
  • Fork 6 You must be signed in to fork a gist
Star You must be signed in to star a gist
What would you like to do?
A script to generate contour plots of Dirichlet distributions
'''Functions for drawing contours of Dirichlet distributions.'''
# Author: Thomas Boggs
from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
_corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
_AREA = 0.5 * 1 * 0.75**0.5
_triangle = tri.Triangulation(_corners[:, 0], _corners[:, 1])
# For each corner of the triangle, the pair of other corners
_pairs = [_corners[np.roll(range(3), -i)[1:]] for i in range(3)]
# The area of the triangle formed by point xy and another pair or points
tri_area = lambda xy, pair: 0.5 * np.linalg.norm(np.cross(*(pair - xy)))
def xy2bc(xy, tol=1.e-4):
'''Converts 2D Cartesian coordinates to barycentric.
`xy`: A length-2 sequence containing the x and y value.
coords = np.array([tri_area(xy, p) for p in _pairs]) / _AREA
return np.clip(coords, tol, 1.0 - tol)
class Dirichlet(object):
def __init__(self, alpha):
'''Creates Dirichlet distribution with parameter `alpha`.'''
from math import gamma
from operator import mul
self._alpha = np.array(alpha)
self._coef = gamma(np.sum(self._alpha)) / \
np.multiply.reduce([gamma(a) for a in self._alpha])
def pdf(self, x):
'''Returns pdf value for `x`.'''
from operator import mul
return self._coef * np.multiply.reduce([xx ** (aa - 1)
for (xx, aa)in zip(x, self._alpha)])
def sample(self, N):
'''Generates a random sample of size `N`.'''
return np.random.dirichlet(self._alpha, N)
def draw_pdf_contours(dist, border=False, nlevels=200, subdiv=8, **kwargs):
'''Draws pdf contours over an equilateral triangle (2-simplex).
`dist`: A distribution instance with a `pdf` method.
`border` (bool): If True, the simplex border is drawn.
`nlevels` (int): Number of contours to draw.
`subdiv` (int): Number of recursive mesh subdivisions to create.
kwargs: Keyword args passed on to `plt.triplot`.
from matplotlib import ticker, cm
import math
refiner = tri.UniformTriRefiner(_triangle)
trimesh = refiner.refine_triangulation(subdiv=subdiv)
pvals = [dist.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]
plt.tricontourf(trimesh, pvals, nlevels, cmap='jet', **kwargs)
plt.xlim(0, 1)
plt.ylim(0, 0.75**0.5)
if border is True:
plt.triplot(_triangle, linewidth=1)
def plot_points(X, barycentric=True, border=True, **kwargs):
'''Plots a set of points in the simplex.
`X` (ndarray): A 2xN array (if in Cartesian coords) or 3xN array
(if in barycentric coords) of points to plot.
`barycentric` (bool): Indicates if `X` is in barycentric coords.
`border` (bool): If True, the simplex border is drawn.
kwargs: Keyword args passed on to `plt.plot`.
if barycentric is True:
X =
plt.plot(X[:, 0], X[:, 1], 'k.', ms=1, **kwargs)
plt.xlim(0, 1)
plt.ylim(0, 0.75**0.5)
if border is True:
plt.triplot(_triangle, linewidth=1)
if __name__ == '__main__':
f = plt.figure(figsize=(8, 6))
alphas = [[0.999] * 3,
[5] * 3,
[2, 5, 15]]
for (i, alpha) in enumerate(alphas):
plt.subplot(2, len(alphas), i + 1)
dist = Dirichlet(alpha)
title = r'$\alpha$ = (%.3f, %.3f, %.3f)' % tuple(alpha)
plt.title(title, fontdict={'fontsize': 8})
plt.subplot(2, len(alphas), i + 1 + len(alphas))
print('Wrote plots to "dirichlet_plots.png".')
Copy link

Really useful - thanks!

Copy link

kwahcarioca commented Aug 29, 2023

Your code is very nice showing how to implement Dirichlet straight from the formula. I ve also tried to experiment it calling scipy.stats.dirichlet library instead. It worked well but we needed to change the tolerance of the xy2bc generator from 1e-4 to 1e-9. Otherwise the assertion of the library code wont let us to run.
if (np.abs(np.sum(x, 0) - 1.0) > 10e-10).any():
raise ValueError("The input vector 'x' must lie within the normal "
"simplex. but np.sum(x, 0) = %s." % np.sum(x, 0))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment