Agent Based Models of Segregation#
Examples and exercises for the tutorial Complexity Science, based on Think Complexity, 2nd edition.
Click here to run this notebook on Colab. For the exercises, use the version with blanks.
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
np.set_printoptions(legacy='1.25')
def decorate(**options):
plt.gca().set(**options)
plt.tight_layout()
Schelling’s model#
In 1969 Thomas Schelling published “Models of Segregation”, which proposed a simple model of racial segregation.
The Schelling model is a grid where each cell represents a house. The houses are occupied by two kinds of agents, labeled red and blue, in roughly equal numbers. About 10% of the houses are empty.
At any point in time, an agent might be happy or unhappy, depending on the other agents in the neighborhood, where the “neighborhood” of each house is the set of eight adjacent cells.
In one version of the model, agents are happy if they have at least three neighbors like themselves, and unhappy if they have fewer than three.
The simulation proceeds by choosing an agent at random and checking to see whether they are happy. If so, nothing happens; if not, the agent chooses one of the unoccupied cells at random and moves.
You will not be surprised to hear that this model leads to some segregation, but you might be surprised by the degree.
Implementing the model#
We will explore Schelling’s model by implementing it in Python. The following function initializes the grid.
def make_grid(n):
"""Make an array with two types of agents.
n: width and height of the array
return: NumPy array
"""
choices = np.array([0, 1, 2], dtype=np.int8)
probs = [0.1, 0.45, 0.45]
return np.random.choice(choices, (n, n), p=probs)
grid = make_grid(n=10)
grid
array([[1, 2, 0, 1, 1, 2, 1, 2, 0, 1],
[2, 1, 1, 2, 2, 2, 1, 2, 1, 0],
[2, 2, 2, 1, 2, 2, 2, 0, 1, 1],
[0, 2, 2, 1, 1, 1, 1, 2, 2, 0],
[2, 2, 2, 1, 1, 1, 1, 2, 2, 2],
[1, 2, 2, 2, 2, 2, 2, 1, 1, 2],
[0, 1, 2, 1, 2, 1, 2, 1, 2, 2],
[0, 0, 1, 1, 1, 2, 1, 1, 0, 1],
[2, 2, 0, 1, 0, 2, 2, 2, 2, 0],
[1, 1, 1, 1, 1, 0, 2, 2, 1, 2]], dtype=int8)
The following function plots a visualization of the grid.
from matplotlib.colors import LinearSegmentedColormap
# make a custom color map
palette = sns.color_palette('muted')
colors = 'white', palette[1], palette[0]
cmap = LinearSegmentedColormap.from_list('cmap', colors)
def draw(grid):
"""Draws the grid.
grid: NumPy array
"""
# Make a copy because some implementations
# of step perform updates in place.
a = grid.copy()
n, m = a.shape
plt.axis([0, m, 0, n])
plt.xticks([])
plt.yticks([])
options = dict(interpolation='none', alpha=0.8)
options['extent'] = [0, m, 0, n]
return plt.imshow(a, cmap, **options)
Here’s what the random initial condition looks like.
im = draw(grid)
The step function#
To simulate a time step in Schelling’s model, we have to
Choose an unhappy agent.
Choose a random empty cell.
Move the agent.
The next few cells develop the functions we’ll need. As an exercise, you will assemble them into a step function.
First I’ll compute empty, which is a boolean array that is True for the empty cells.
empty = (grid==0)
empty.sum()
14
To find the unhappy agents, we have to compute the number of same-color neighbors for each agent. We’ll do that using 2-D correlation.
As an example, we can compute the number of empty neighbors for each cell.
from scipy.signal import correlate2d
kernel = np.array([[1, 1, 1],
[1, 0, 1],
[1, 1, 1]], dtype=np.int8)
options = dict(mode='same', boundary='wrap')
correlate2d(empty, kernel, **options)
array([[1, 1, 0, 1, 1, 1, 1, 1, 1, 2],
[1, 1, 1, 1, 0, 0, 1, 2, 3, 1],
[3, 1, 0, 0, 0, 0, 1, 0, 3, 3],
[1, 1, 0, 0, 0, 0, 1, 1, 2, 1],
[2, 1, 0, 0, 0, 0, 0, 0, 1, 2],
[1, 1, 0, 0, 0, 0, 0, 0, 0, 1],
[2, 3, 1, 0, 0, 0, 0, 1, 1, 3],
[3, 3, 2, 2, 1, 1, 0, 1, 1, 4],
[3, 3, 1, 2, 1, 2, 1, 1, 2, 2],
[1, 2, 2, 3, 2, 1, 1, 1, 2, 2]], dtype=int8)
The following function computes the number of same-color neighbors for each agent.
def compute_num_same(grid):
"""For each cell, the number of same-color neighbors.
grid: NumPy array
return: new NumPy array
"""
red = grid==1
blue = grid==2
# count red neighbors, blue neighbors, and total
num_red = correlate2d(red, kernel, **options)
num_blue = correlate2d(blue, kernel, **options)
num_same = np.where(red, num_red, num_blue)
return num_same
Here’s what that looks like for the initial grid.
# for each cell, compute the fraction of neighbors with the same color
num_same = compute_num_same(grid)
num_same
array([[4, 1, 2, 5, 3, 3, 1, 3, 4, 4],
[3, 2, 3, 3, 5, 5, 1, 2, 3, 2],
[3, 5, 4, 3, 4, 4, 4, 4, 2, 2],
[6, 7, 5, 4, 6, 5, 3, 4, 4, 5],
[5, 6, 6, 3, 5, 5, 4, 4, 5, 4],
[1, 5, 5, 5, 3, 4, 3, 3, 2, 5],
[3, 2, 3, 3, 4, 2, 3, 4, 2, 2],
[3, 3, 4, 4, 4, 4, 3, 2, 4, 0],
[2, 1, 1, 6, 2, 3, 5, 4, 3, 3],
[3, 3, 4, 5, 4, 4, 6, 5, 1, 2]], dtype=int8)
The mean of this array measures the average level of segregation.
num_same.mean()
3.52
locs_where is a wrapper on np.nonzero that returns a list of pairs, where each pair is the coordinate of an agent that meets a condition.
def locs_where(condition):
"""Find cells where a boolean array is True.
condition: NumPy array
return: list of coordinate pairs
"""
ii, jj = np.nonzero(condition)
return list(zip(ii, jj))
An agent is “unhappy” if they have fewer than 3 same-color neighbors.
Here are the locations of the unhappy agents.
threshold = 3
unhappy_locs = locs_where(num_same < threshold)
unhappy_locs
[(0, 1),
(0, 2),
(0, 6),
(1, 1),
(1, 6),
(1, 7),
(1, 9),
(2, 8),
(2, 9),
(5, 0),
(5, 8),
(6, 1),
(6, 5),
(6, 8),
(6, 9),
(7, 7),
(7, 9),
(8, 0),
(8, 1),
(8, 2),
(8, 4),
(9, 8),
(9, 9)]
During each time step, we’ll choose a random unhappy agent.
The following function chooses a random tuple from a list.
def random_loc(locs):
"""Choose a random element from a list of tuples.
locs: list of tuples
return: tuple
"""
index = np.random.choice(len(locs))
return locs[index]
source is the random unhappy agent who decides to move.
source = random_loc(unhappy_locs)
source
(6, 8)
dest is the random empty cell they choose to move to.
empty_locs = locs_where(empty)
dest = random_loc(empty_locs)
dest
(8, 2)
The following function swaps source and dest.
def move(grid, source, dest):
"""Swap the agents at source and dest.
grid: NumPy array
source: location tuple
dest: location tuple
"""
grid[dest], grid[source] = grid[source], grid[dest]
move(grid, source, dest)
Exercise: Pull the code from the previous cells into a function that chooses an unhappy agent and moves them to a random empty cell. It should compute and return the mean of num_same.
Note: If there are no unhappy cells, step should do nothing.
def step(grid, threshold=3):
"""Simulate one time step.
grid: NumPy array
threshold: number of same-color neighbors needed to be happy
return: average number of same-color neighbors
"""
num_same = compute_num_same(grid)
# FILL THIS IN!
return num_same.mean()
# Solution
def step(grid, threshold=3):
"""Simulate one time step.
grid: NumPy array
threshold: number of same-color neighbors needed to be happy
return: average number of same-color neighbors
"""
num_same = compute_num_same(grid)
unhappy_locs = locs_where(num_same < threshold)
if len(unhappy_locs):
source = random_loc(unhappy_locs)
empty_locs = locs_where(grid==0)
dest = random_loc(empty_locs)
move(grid, source, dest)
return num_same.mean()
Test your step function here.
step(grid)
3.51
Then run the following cells to see how segregation changes over time.
def decorate_seg():
decorate(xlabel='Number of steps',
ylabel='Average number of same-color neighbors',
title='Schelling model')
grid = make_grid(n=10)
segs = [step(grid) for i in range(200)]
plt.plot(segs)
decorate_seg()
Exercise: Experiment with different values of threshold and see what effect they have on the result.
Bigger grid#
Let’s see how that looks on a bigger grid.
grid = make_grid(n=50)
im = draw(grid)
segs = [step(grid, threshold=3) for i in range(3000)]
plt.plot(segs)
decorate_seg()
im = draw(grid)
The following figure shows three steps in the evolution of the grid.
np.random.seed(17)
grid = make_grid(n=50)
# draw the initial grid
plt.figure(figsize=(9,3))
plt.subplot(1,3,1)
draw(grid)
# first update
plt.subplot(1,3,2)
for i in range(1000):
step(grid)
draw(grid)
# second update
plt.subplot(1,3,3)
for i in range(1000):
step(grid)
draw(grid)
plt.tight_layout()
Exercise: One more time, try out different values of threshold. What value yields the biggest difference between threshold and the average number of same-color neighbors?
Animation#
The following functions animate the evolution of the grid.
from time import sleep
from IPython.display import clear_output
def animate(grid, frames, interval=None):
"""Animate the automaton.
grid: NumPy array
frames: number of frames to draw
interval: time between frames in seconds
"""
plt.figure()
try:
for i in range(frames-1):
draw(grid)
plt.show()
if interval:
sleep(interval)
step(grid)
clear_output(wait=True)
draw(grid)
plt.show()
except KeyboardInterrupt:
pass
grid = make_grid(n=50)
# animate(grid, frames=1000)
Exercise: Experiment with different starting conditions: for example, more or fewer empty cells, or unequal numbers of red and blue agents.
The Big Sort#
Bill Bishop, author of The Big Sort, argues that American society is increasingly segregated by political opinion, as people choose to live among like-minded neighbors.
The mechanism Bishop hypothesizes is not that people, like the agents in Schelling’s model, are more likely to move if they are isolated, but that when they move for any reason, they are likely to choose a neighborhood with people like themselves.
Let’s write a version of Schelling’s model to simulate this kind of behavior and see if it yields similar degrees of segregation.
We’ll start with a random initial grid.
grid = make_grid(n=10)
im = draw(grid)
Write a few lines of code to choose a random occupied location, and call it source.
empty = (grid==0)
empty_locs = locs_where(empty)
empty_locs
[(0, 7), (2, 5), (2, 8), (4, 3), (5, 2), (8, 3), (8, 4), (8, 8)]
full_locs = locs_where(~empty)
source = random_loc(full_locs)
source
(9, 0)
Here’s a list of random empty locations the source agent could move to.
num_comps = 2
dests = [random_loc(empty_locs) for i in range(num_comps)]
dests
[(8, 3), (8, 3)]
The following cell computes the number of neighbors in each location that are the same color as source.
same_color = (grid==grid[source])
num_same = correlate2d(same_color, kernel, **options)
num_same
array([[5, 6, 5, 2, 3, 2, 3, 3, 1, 4],
[2, 5, 6, 5, 3, 2, 1, 2, 1, 2],
[3, 5, 4, 2, 1, 1, 2, 3, 2, 2],
[4, 3, 4, 2, 1, 1, 3, 4, 4, 4],
[4, 3, 3, 1, 2, 4, 5, 6, 5, 3],
[7, 5, 4, 3, 2, 4, 4, 5, 6, 6],
[4, 3, 2, 3, 3, 6, 6, 5, 4, 3],
[6, 5, 4, 3, 3, 4, 5, 4, 6, 4],
[5, 3, 2, 1, 3, 4, 6, 4, 5, 3],
[4, 5, 6, 3, 4, 3, 5, 4, 3, 4]], dtype=int8)
Now we can make a list with one tuple per destination, where the first element in each tuple is the number of same-color neighbors.
choices = [(num_same[dest], np.random.random(), dest) for dest in dests]
for choice in choices:
print(choice)
(1, 0.9606424068236474, (8, 3))
(1, 0.6324782497028868, (8, 3))
Using the max function, we get the location with the highest number of same-color neighbors, with ties broken at random.
num, rand, dest = max(choices)
num, rand, dest
(1, 0.9606424068236474, (8, 3))
move(grid, source, dest)
Now pull all of that code into a step function.
def step(grid, num_comps=2):
"""Simulate one time step.
grid: NumPy array
num_comps: number of possible destinations the agent looks at
return: average number of same-color neighbors
"""
# FILL THIS IN
return compute_num_same(grid).mean()
# Solution
def step(grid, num_comps=2):
"""Simulate one time step.
grid: NumPy array
num_comps: number of possible destinations the agent looks at
return: average number of same-color neighbors
"""
empty = (grid==0)
empty_locs = locs_where(empty)
full_locs = locs_where(~empty)
source = random_loc(full_locs)
same_color = (grid==grid[source])
num_same = correlate2d(same_color, kernel, **options)
dests = [random_loc(empty_locs) for i in range(num_comps)]
choices = [(num_same[dest], np.random.random(), dest) for dest in dests]
_, _, dest = max(choices)
move(grid, source, dest)
return compute_num_same(grid).mean()
step(grid)
3.45
grid = make_grid(n=10)
segs = [step(grid) for i in range(200)]
plt.plot(segs)
decorate_seg()
Exercise: Experiment with different values of num_comps and see what effect they have on the result.
Bigger grid#
Let’s see how that looks on a bigger grid.
grid = make_grid(n=50)
im = draw(grid)
segs = [step(grid) for i in range(2000)]
plt.plot(segs)
decorate_seg()
im = draw(grid)
The following figure shows three steps in the evolution of the grid.
np.random.seed(17)
grid = make_grid(n=50)
# draw the initial grid
plt.figure(figsize=(9,3))
plt.subplot(1,3,1)
draw(grid)
# first update
plt.subplot(1,3,2)
for i in range(1000):
step(grid)
draw(grid)
# second update
plt.subplot(1,3,3)
for i in range(1000):
step(grid)
draw(grid)
plt.tight_layout()
Exercise: One more time, try out different values of num_comps.
Animation#
Let’s see what the animation looks like.
grid = make_grid(n=50)
# animate(grid, frames=1000)
Copyright 2019–2026 Allen Downey
Code license: MIT License
Text license: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International