I came across a little problem that lends itself nicely to visualization. And, it’s the kind of visualization that many of us have seen before: Venn diagrams.
import matplotlib.pyplot as plt
%matplotlib inline
# create some circles
circle1 = plt.Circle((-.5,0), 1, color='r', alpha=.2)
circle2 = plt.Circle(( .5,0), 1, color='b', alpha=.2)
# create a figure
fig = plt.figure(figsize=(8,8))
ax = fig.gca()
# add them to the figure
ax = fig.gca()
ax.add_artist(circle1)
ax.add_artist(circle2)
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_aspect('equal')
# display it
plt.plot();
Almost certainly you’ve seen a Venn diagram used to visualize the partitioning of a universe that has two possible events which may (or may not) occur. With two events, there are four atomic or primitive events \(AB\), \(A{\lnot}B\), \(\lnot{A}{B}\), \(\lnot{A}\lnot{B}\). These possibilities, like components of a Boolean logic expression, can be displayed in a pseudo-truth table:
| Event A | Event B | Primitive Event | Primitive \(P()\) | Example \(P()\) |
|---|---|---|---|---|
| \(T\) | \(T\) | \(AB\) | \(P(AB)\) | .25 |
| \(T\) | \(F\) | \(A{\lnot}B\) | \(P(A{\lnot}B)\) | .25 |
| \(F\) | \(T\) | \(\lnot{A}{B}\) | \(P(\lnot{A}{B})\) | .25 |
| \(F\) | \(F\) | \(\lnot{A}\lnot{B}\) | \(P(\lnot{A}\lnot{B})\) | .25 |
I’m not really harping on this to be daft. I’m also using this notebook to test out my ipython \(\rightarrow\) html conversion system. So far, these cells are testing markdown tables and \(\LaTeX\) symbols. If you chuckle at the meta-structure of that sentence, I’ll forgive you.
Fairly often, Venn Diagrams are introduced in a discussion when the events, \(A\) and \(B\) are independent. Two events are independent, mathematically, if and only if \(P(AB)=P(A)P(B)\). When events are independent, knowning \(P(A)\) and \(P(B)\) is sufficent to determine the probabilities of each primitive event. We just saw how to compute \(P(AB)\). To compute \(P(\lnot{A}\lnot{B})\) we notice that the sum of all the events must be 1.0. Then we give a name to the event where either \(A\), \(B\), or both happens: \(A\cup{B}\). So, we can subtract from one to get \(P(\lnot{A}\lnot{B})=1-P(A\cup{B})\) . For reference, \(P(A\cup{B})=P(A)+P(B)-P(AB)\). Back to the main point: if the events are independent, we can compute all the other primative event probabilities.
Often though, events in the real world are not independent. In that case, we need to fill in the complete table of probabilities. We can do that reasonably in a table for two or three events. But for more than two or three events, it becomes difficult very quickly. That’s where Venn diagrams can come back into the picture. You may have also seen Venn diagrams for three events. I won’t draw one here, but it looks a lot like the one above with a third circle.
Now, I had never actually come to the realization that you can simple keep adding circles for more events and keep getting a complete representation of all the possible primitive events (occurance/non-occurance of the events). Once I knew that was possible, I had a great idea to visually represent the probabilities for all the possible states of the world using these (extended) Venn diagrams.
What we’re going to do is use grayscale values to color the primitive events within the probability space.
Creating Lots of Circles¶
To start out (and to add a bit of color), here’s how you can draw the circles that make up the Venn diagram in pure matplotlib. Turns out we need some more tools to do the coloring, so we won’t follow this up. But, I made it and it looks pretty. It reminds me of Spirograph.
from matplotlib.collections import PatchCollection
# walk around a clock from 3:00 (the positive x-axis) backwards, stop one step shy of full circle
# take n number of steps
# recall that x = r cos(theta) and y = r sin(theta) from trig (polar<->rectangular coordinates)
n = 20
thetas = np.linspace(0, 2*np.pi - (2*np.pi/n) , n)
r = 1.0
# create the circles
patches = [plt.Circle((r*np.cos(th), r*np.sin(th)), 2.0) for th in thetas]
collection = PatchCollection(patches, alpha=.2)
collection.set_array(np.arange(n))
# create the figure and add the circles
fig = plt.figure(figsize=(8,8))
ax = fig.gca()
plt.gca().add_collection(collection)
# scale it appropriately
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.gca().set_aspect('equal')
Required Tools¶
To get started, we need to import a few modules. I had originally wanted to use a pure matplotlib solution. However, it turned out that the fundamental task of slicing up pre-drawn circles and figuring out the constituent patches (the different colored areas in the Venn diagram above) was not immediately possible. So, I asked on stackoverflow and got a pointer to shapely and descartes. Turns out that shapely has an API for exactly this sort of task. In fact, you can explictly take the unions and intersections of graphical objects and compute the patches where those set operations occur.
Huzzah: the right tool for the job!
import itertools as it
import numpy as np
np.set_printoptions(precision=4, suppress=True)
import matplotlib.pyplot as plt
import matplotlib as mpl
import descartes
import shapely.geometry as sg
import shapely.ops as so
%matplotlib inline
I define a few functions to give me a basic vocabulary for the subtasks in this problem: creating a bunch of (overlapping) circles, coloring patches of the output graph from patches and colors, computing a running intersection from a GeometryCollection (a Shapely critter) of circles, and creating patch colors from probabilities.
def makeShapelyCircles(numEvents, size=2.0):
thetas = np.linspace(0, 2*np.pi - (2*np.pi/numEvents), numEvents)
r = np.ones(thetas.shape)
centers = np.column_stack([r*np.cos(thetas), r*np.sin(thetas)])
print len(centers)
return [sg.Point(*c).buffer(float(size)) for c in centers]
# patch = [p1, p2, p3, ...]
# color = [c1, c2, c3, ...]
def addPatchesAndColors(ax, patch, color):
addPatchColorPairs(ax, it.izip(patch, color))
# pAndC = [(patch, color), (patch, color)]
def addPatchColorPairs(ax, pAndC):
for patch, color in pAndC:
ax.add_patch(descartes.PolygonPatch(patch,
fc=color,
ec='w'))
# for bonus points, you could try to replace this with
# a clever reduce/functools.reduce call
def chainIntersection(comps):
comps = iter(comps)
intersection = comps.next()
for comp in comps:
intersection = comp.intersection(intersection)
return intersection
A quick note here, if we have \(N\) events there are \(2^N\) primitive events (this is the size of the powerset created from our \(N\) events). We express each possible primitive event as a (reverse) binary representation of a number from \(0\) to \(2^N-1\). This may seem a bit odd (not even? ha!). But, it allows us to quickly add up primitive events as they would occur in the pseudo-truth table above using a numpy array. [Here the array creation is "hard", but the use of it is "easy". The creation process is really just making a series of indicators for which events occur and don’t occur on a given line. The real key is numpy’s indices routine. I recommend playing around with it at the interactive prompt. The [2] is used since our events either occur or don’t occur. Play around with the other values and see what happens! That’s part of the beauty of python.
eventCt = 2
pwrSetCt = 2**eventCt
atomicEvents = (np.indices([2] * eventCt, dtype=np.uint8).reshape(eventCt, pwrSetCt)[::-1].T)
atomicProbs = np.array([.25, .25, .25, .25])
# this is actually the first time (i can recall) doing a nested
# multiple target assignment in python!
for idx, (event, probs) in enumerate(zip(atomicEvents, atomicProbs)):
print idx, event, prob
# since these are disjoin events, to get P(A) we just add up everywhere A happens
primitiveEvents[:,0] # 0-th column is "A"
primitiveEvents[:,0] * primitiveProbs # contribution to total P(A)
(primitiveEvents[:,0] * primitiveProbs).sum()
# more examples of the indexing scheme:
# index ---> event indicator ---> which events occurred
# 0 --> 0,0,0,0 []
# 1 --> 1,0,0,0 [0]
# 2 --> 0,1,0,0 [1]
# 3 --> 1,1,0,0 [0,1]
# 4 --> 0,0,1,0 [2]
# 5 --> 1,0,1,0 [0,2]
# 15 --> 1,1,1,1 [0,1,2,3]
def makePatchColors(events, probs):
eventCt = len(events)
pwrSetCt = 2**eventCt
primEvents = (np.indices([2] * eventCt, dtype=np.uint8)
.reshape(eventCt, pwrSetCt)[::-1].T)
GC = sg.collection.GeometryCollection
res = []
for anEvent, prob in zip(primEvents, probs):
print anEvent, prob
# UGLY!
indices = [idx for idx, val in enumerate(anEvent) if val]
if not indices: continue
components = GC([events[idx] for idx in indices])
patch = chainIntersection(components)
res.append((patch, str(1.0-prob))) # color/prob inverted
return res
Finally, a few utilities to (1) actually put the Venn diagram (and a color bar) on the matplotlib figure and (2) create random probabilities for some testing.
def makeRandomProbs(numEvents, bigNull = False):
probs = np.random.uniform(0.0, 10.0, 2**numEvents)
if bigNull: # hack null event to big prob.
probs[0] = .5 * probs[1:].sum()
probs = probs / probs.sum()
print "sorted probs: %s sum: %d" % (np.array(sorted(probs)), probs.sum())
return probs
def addColorbar(ax, cmapName):
cbax, cbkw = mpl.colorbar.make_axes(ax)
mpl.colorbar.ColorbarBase(cbax,
cmap=plt.cm.get_cmap(cmapName),
norm=mpl.colors.Normalize(clip=True,
vmin=0.0,
vmax=1.0),
**cbkw)
def addVennDiagram(ax, probs):
circs = makeShapelyCircles(np.log2(len(probs)))
pacs = makePatchColors(circs, probs)
addPatchColorPairs(ax, pacs)
ax.set_axis_bgcolor(str(1.0-probs[0]))
# fig.colorbar(mplPatches)
addColorbar(ax, "Greys")
Some Simple Testing¶
One last helper to make testing easier.
def createWith(probs):
#
# control display
#
fig = plt.figure(figsize=(8,8))
ax = plt.gca()
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.set_aspect('equal')
addVennDiagram(ax, probs)
ax.plot()
For my first example, we’ll do the uniform probabilites as above (for two and three events).
numEvents = 2
createWith(np.ones(2**numEvents)/2**numEvents)
numEvents = 3
createWith(np.ones(2**numEvents)/2**numEvents)
Let’s test the colors to make sure the grayscale color bar matches with the graph.
createWith(np.array([0.0, 0.25, .5, .75])) # doesn't have to actually be legit values
The colors seem reasonable, but we might actually do better with an RGB colors to allow for better discrimination from (.75 to 1.0) and (.25, .40). In these ranges, it is a bit hard to see differences.
And a Less Simple Example¶
numEvents=2
createWith(makeRandomProbs(numEvents, bigNull=True))
Wrap Up¶
I never really mentioned why I needed to deal with tables of probabilities. Suffice it to say, I’ll post on that sometime soon. But, here’s the flavor. Using a probabilistic model, I compute the probabilities of many primitive events and I wanted to compare them. Since there are typically a half-dozen events, evaluating \(2^{6}\) table entries seemed up appealing. Instead, I can look at something like this:
numEvents=5
createWith(makeRandomProbs(numEvents))
Additional Resources
You can grab a copy of this notebook.
Even better, you can view it using nbviewer.
License
Unless otherwise noted, the contents of this notebook are under the following license. The code in the notebook should be considered part of the text (i.e., licensed and treated as as follows).
DrsFenner.org Blog And Notebooks by Mark and Barbara Fenner is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Permissions beyond the scope of this license may be available at drsfenner.org/blog/about-and-contacts.