Mar 16

Probability Podcast Ep2: Imprecise probabilities with Gert de Cooman


I happened to be travelling through Brussels, so I stopped by Ghent, the world hotspot for research into imprecise probabilities, and setup an interview with Gert de Cooman. Gert has been working in imprecise probabilities for more than twenty years, is a founding member and former President of SIPTA, the Society for Imprecise Probability: Theories and Applications, and has helped organize many of the ISIPTA conferences and SIPTA Schools.

Topics include fair betting rates, Dutch books, Monte Carlo methods, Markov chains, utility, and the foundations of probability theory. We had a rich, wide-ranging discussion. You may need to listen two (or more!) times to process everything.

Episode on SoundCloud

Jun 14

A new way to visualize content

Right now I’m working on a project that involves new ways to view units of content and the relationships between them. I’ve posted the comic I worked on, it has a number of stats references throughout. This is early alpha stages for the software, you may run into issues. To see the relationships, go to the puffball menu and make sure that “Show relationships” is clicked.

Feb 13

Population simulation leads to Valentine’s Day a[R]t

Working on a quick-and-dirty simulation of people wandering around until they find neighbors, then settling down. After playing with the coloring a bit I arrived at the above image, which I quite like. Code below:

# Code by Matt Asher for
# Feel free to modify and redistribute, but please keep this notice 
maxSettlers = 150000
# Size of the area
areaW = 300
areaH = 300
# How many small movements will they make to find a neighbor
maxSteps = 200
# Homesteaders, they don't care about finding a neighbor
numbHomesteaders = 10
areaMatrix = matrix(0, nrow=areaW, ncol=areaH)
# For the walk part
adjacents = array(c(1,0,1,1,0,1,-1,1,-1,0,-1,-1,0,-1,1,-1), dim=c(2,8))
# Is an adjacent cell occupied?
hasNeighbor <- function(m,n,theMatrix) {
	toReturn = FALSE
	for(k in 1:8) {
		yCheck = m + adjacents[,k][1]
		xCheck = n + adjacents[,k][2]
		if( !((xCheck > areaW) | (xCheck < 1) | (yCheck > areaH) | (yCheck < 1)) ) {
			if(theMatrix[yCheck,xCheck]>0) {
				toReturn = TRUE
# Main loop
for(i in 1:maxSettlers) {
	steps = 1
	xPos = sample(1:areaW, 1)
	yPos = sample(1:areaH, 1)
	if(i <= numbHomesteaders) {
		# Seed it with homesteaders
		areaMatrix[xPos,yPos] = 1
	} else {
		if(areaMatrix[xPos,yPos]==0 & hasNeighbor(xPos,yPos,areaMatrix)) {
			areaMatrix[xPos,yPos] = 1
		} else {
			spotFound = FALSE
			outOfBounds = FALSE
			while(!spotFound & !outOfBounds & (steps<maxSteps)) {
				# Look for a new location in one of adjacent 9 cells, while still in area
				steps = steps + 1
				movement = adjacents[,sample(1:8,1)]
				xPos = xPos + movement[1]
				yPos = yPos + movement[2]
				if( (xPos > areaW) | (xPos < 1) | (yPos > areaH) | (yPos < 1)) {
					outOfBounds = TRUE
				} else if(hasNeighbor(xPos,yPos,areaMatrix) ) {
					areaMatrix[xPos,yPos] = steps
					spotFound = TRUE
image(areaMatrix, col=rev(rgb(seq(0.01,1,0.01),seq(0.01,1,0.01),seq(0.01,1,0.01))))
# I think this version looks nicer!
# areaMatrix[areaMatrix !=0] = 1
# image(areaMatrix, col=rev(rgb(.5,0,seq(0.2,1,0.2))))

Dec 12

Let it snow!

A couple days ago I noticed a fun piece of R code by Allan Roberts, which lets you create a digital snowflake by cutting out virtual triangles. Go give it a try. Roberts inspired me to create a whole night sky of snowflakes. I tried to make the snowfall look as organic as possible. There are lots of options to adjust. Here’s the code, have fun and Happy Holidays!

# Code by Matt Asher for
# Feel free to modify and redistribute 
# How many flakes do you want to fall?
flakes = 100
# Width and height of your space
width = 800
height = 600
# Initial wind
wind = 0
# Setup the background of the plot and margins
par(bg = "black")
plot(0, 0, col="black", pch=".", xlim=c(0,width), ylim=c(0,height), axes=F)
for(i in 1:flakes) {
    startY = height
    startX = runif(1,1,width)
    xPos = startX
    yPos = startY
    for(j in 1:height) {
		# Optional drift in wind
		wind = wind + rcauchy(1,0,.05)
		# Update snowflake position
        xPos = xPos + rnorm(1,.1,1.5)
        yPos = yPos - runif(1,4,20)
        # Are we in the space, if so display it
        if(xPos &gt; 0 &amp;&amp; xPos &lt;= width &amp;&amp; yPos &gt; 0 &amp;&amp; yPos &lt;= height) {
            points(round(xPos), round(yPos), col="white", pch=8)
            # System dely, slows down the flakes

Oct 12

Comic with stats discussion

I recently finished work on the first issue of a graphic novel. It’s in the form of a fictional first person narrative. The story isn’t directly about statistics, but there are a few digressions on the subject. Here are some samples, make sure to click on the images for a larger view:

If you’re interested, head over to and pick up a copy. Here’s the order page. The comic comes with a full money-back guarantee, including shipping. You don’t even have to send back your copy to claim the refund.