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.

Jan 12

R A Fisher illustration

Ronald Aylmer Fisher, statistics badass. Illustration by Rachelle Scarfó for a project I was working on.

Sep 10

Weekend art in R (Part 4)

Computer creations are perfect by design. We put in numbers, and if all goes well we get out an exact result. If we want a line, we want it perfectly straight. If we want a circle, it should conform to the platonic ideal of a circle. From a mathematical standpoint, these perfect shapes and precisely computed numbers are ideal.

Someday, perhaps, we will have true fuzzy computation built right into our hardware. For now, it takes considerable effort to achieve just the right level of imperfection needed for simulating mistakes, or any organic processes.

I sent each of the circles shown above on a random walk. That part was easy, getting each circle to end up where it started (and close the loop) took a bit more effort. To vary the “wigglyness” of the lines, adjust the “sd” parameter in “rnorm”. To change how quickly randomness tapers off, change the “4″ in “i/4″. Here is my code:

# Circle lengths
j = seq(0.1,1.9,.08)
par(bg = "black")
# How many dots around the circle?
dots = 1000
# Create an offkilter circle
rads = seq(0,2*pi,2*pi/dots)
for(aLength in j) {
	# Pick a random color
	myCol = paste("#",paste(sample(c(1:9,"A","B","C","D","E","F"),6,replace=T),collapse=""),collapse="",sep="")
	# Start at length = 1, then walk.
	myLength = rep(aLength,dots)
	for(i in 2:dots) {
		myLength[i] = myLength[(i-1)] + rnorm(1,0,sd=.005)
		# Closer we are to end, faster we return to where started so circle closes
		dist = aLength - myLength[i]
		myLength[i] = aLength - (dist*((dots-(i/4))/(dots)))
	for(i in 1:dots) {

What do your circles look like?