Learning to Be Thoughtless in Python with Mesa
I recently read Who's Counting? Uniting Numbers and Narratives with Stories from Pop Culture, Puzzles, Politics, and More by John Allen Paulos. In the Chapter 6 section titled "Voting Blocs: Red States, Blue State, and a Model for Thoughtless Voting," Paulos mentions an interesting paper by Joshua Epstein, Learning To Be Thoughtless: Social Norms And Individual Computation. In this paper, Epstein proposes a simple but compelling agent-based model for how social norms evolve in a community. This model is interesting because it captures two fundametal forces that motivate human behavior quite elegantly:
- the social influence the behavior of others exerts on our own, and
- laziness.
Paulos uses Epstein's model as a potential explanation for the increasing polarization of American Politics, which, while plausible, is not the focus of this post. After following the citation to Epstein's paper, I thought it would be fun and instructive to implement this model in Python and reproduce the main results of the paper.
Epstien's paper uses the decision to drive on the left- or right-hand side of the road as the social norm under analysis. For the purpose of this post, we will use a boolean flag norm to indicate whether or not an agent conforms to the social norm (defined, for instance, by driving on the right-hand side of the road in the United States) or does not (driving on the left-hand side of the road in the United States).
Epstein names his agents "lazy statisticans." These lazy statisticians determine whether or not they follow the social norm based on the norm-following behavior of the other agents in their vicinity. The size of the neighborhood that influences a statistician's behavior is their radius. With this terminology, we can begin to implement lazy statisticians in Python. Throughout this post we use Mesa, an agent-based modeling framework that
allows users to quickly create agent-based models using built-in core components (such as spatial grids and agent schedulers) or customized implementations; visualize them using a browser-based interface; and analyze their results using Python’s data analysis tools.
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from fastprogress.fastprogress import progress_bar
from matplotlib import pyplot as plt
import mesa
import seaborn as sns
plt.rc("figure", figsize=(8, 6))
plt.rc("font", size=14)
sns.set(color_codes=True)
CMAP = "viridis"
Each lazy statistician is a Mesa Agent. Following the paper, if no initial value for norm is past, we choose one with a coin flip, and if no initial radius is passed, we choose one uniformly between one and fitfy.
class LazyStatistician(mesa.Agent):
def __init__(self, unique_id, model, norm=None, radius=None, tol=0.05):
super().__init__(unique_id, model)
self.norm = self.random.choice([False, True]) if norm is None else norm
self.radius = self.random.randint(1, 51) if radius is None else radius
self.tol = tol
When a lazy statistician updates their norm-following behavior, they do so based on the proportion of norm-following agents in their vicinity as determined by their radius (optionally offset by Δ).
class LazyStatistician(LazyStatistician):
def get_norm_pct(self, *, Δ=0):
radius = self.radius + Δ
nhbrs = self.model.grid.iter_neighbors(self.pos, moore=True, radius=radius)
return sum(nhbr.norm for nhbr in nhbrs) / (2 * radius + 1)
When a lazy statistician updates their norm-following behavior, they do so in two steps. First they update the radius of neighboring agents that influence them, then they update their norm.
When updating their radius, a lazy statistician first compares the proportion of norm-followers in their current neighborhood, defined by their radius, to the proportion of norm-followers in the neighborhood defined by radius + 1. If these not close (for the purposes of this post, we mostly follow the guideline of the paper and define "close" to be within 5% of each other), the lazy statistician increases their radius by one, since more data is better. If they are close, the lazy statistician compares the proportion of norm-followers in their current neighborhood to the proportion of norm-followers in the neighborhood by radius - 1. If these are close, the lazy statisician decreases their radius by one (due to their laziness).
class LazyStatistician(LazyStatistician):
def update_radius(self):
norm_pct = self.get_norm_pct()
more_pct = self.get_norm_pct(Δ=1)
if abs(norm_pct - more_pct) > self.tol:
self.radius += 1
elif self.radius > 1:
less_pct = self.get_norm_pct(Δ=-1)
if abs(norm_pct - less_pct) <= self.tol:
self.radius -= 1
After updating their radius, a lazy statistician updates their norm to agree with the majority of the agents in their neighborhood (including themselves). Note that including themsevles means that there are an odd number of agents (2 * radius + 1) agents in their neighborhood, so their will never be an ambiguous case where exactly 50% of agents in the neighborhood are norm-followers.
class LazyStatistician(LazyStatistician):
def update_norm(self):
self.norm = self.get_norm_pct() > 0.5
def step(self):
self.update_radius()
self.update_norm()
If there are $N$ agents, at each time step $N$ agents are drawn randomly with replacement and their status is updated as described above. Mesa does not include this sort of bootstrap scheduler out-of-the-box, but it is simple to implement.
class BootstrapActivation(mesa.time.BaseScheduler):
def step(self):
for _ in range(self.get_agent_count()):
self.model.random.choice(self.agents).step()
self.steps += 1
self.time += 1
With LazyStatistician and BoostrapActivation defined, we are ready to implement our model.
class ThoughtlessModel(mesa.Model):
def __init__(self, n_agent, norm=None, seed=None,
agent_cls=LazyStatistician, **agent_kwargs):
self.grid = mesa.space.SingleGrid(n_agent, 1, torus=True)
self.schedule = BootstrapActivation(self)
for unique_id in range(n_agent):
agent = agent_cls(unique_id, self, norm=norm, **agent_kwargs)
self.schedule.add(agent)
self.grid.place_agent(agent, (unique_id, 0))
self.datacollector = mesa.DataCollector(agent_reporters={"norm": "norm", "radius": "radius"})
def step(self):
self.datacollector.collect(self)
self.schedule.step()
The most important feature of of this model is that the agents are placed on a circle (as indicated by the torus=True keyword argument to SingleGrid.) Another point of note is that the constructor take the agent class as a parameter. This flexibility will be useful in reconstructing examples from the paper that add noise and/or temporal shocks to the lazy statistician's update rules.
With ThoughtlessModel in hand, we can start reproducing Epstein's computational results. The paper presents six runs of the model, exploring
stylized facts regarding the evolution of norms: [l]ocal conformity, global diversity, and punctuated equilibria.
Run 1. Monolithic Social Norm, Individual Computing Dies Out¶
Following Epstein, for the first simulation we set each lazy statistician's initial norm to False. In this case, we anticipate that norm should remain False for all agents for all time and each agent's radius should shrink to the minimum value of one.
We choose a seed value for reproducibility.
SEED = 1234567890
Instead of simulating 190 lazy statisticians as in the paper, we simulate 500. We find that this yields interesting outcomes in more complex situations far more often (as the seed varies). As in the paper, we simulate 275 time steps.
N_AGENT = 500
N_STEP = 275
def simulate(model, n_step):
for _ in progress_bar(range(n_step)):
model.step()
return model
model_1a = simulate(ThoughtlessModel(N_AGENT, norm=False, seed=SEED), N_STEP)
Visualizing the evolution of each agent's norm and radius over time confirms our expectations.
def plot_results(model, cmap=CMAP, figsize=(16, 6)):
agent_df = model.datacollector.get_agent_vars_dataframe()
fig, axes = plt.subplots(ncols=2, sharex=True, sharey=True,
figsize=figsize, gridspec_kw={"width_ratios": (0.85, 1)})
norm_ax, radius_ax = axes
sns.heatmap(agent_df["norm"].unstack(),
cmap=cmap, vmin=0, vmax=1, cbar=False,
ax=norm_ax)
norm_ax.set_xticks([])
norm_ax.set_xlabel("Agent")
norm_ax.set_yticks([])
norm_ax.set_ylabel("Step")
norm_ax.set_title("norm")
rhm = sns.heatmap(agent_df["radius"].unstack(),
cmap=cmap, vmin=1, vmax=10,
ax=radius_ax)
radius_ax.set_xticks([])
radius_ax.set_xlabel("Agent")
radius_ax.set_yticks([])
radius_ax.set_ylabel(None)
radius_ax.set_title("radius")
fig.tight_layout()
return fig, axes
plot_results(model_1a);
Though the paper does show this scenario, we get an analagous result when all of the initial norms are True.
model_1b = simulate(ThoughtlessModel(N_AGENT, norm=True, seed=SEED), N_STEP)
plot_results(model_1b);
In each of these cases, "individual computing dies out" is reflected by the fact that every agent's radius eventually decreases to one, the smallest possible.
Run 2. Random Initial Norms, Individual Computing At Norm Boundaries¶
Run one verified our intuition in the basic but unrealistic case of a uniform ("monolithic") initial value of norm. In run two we explore the more interesting and realistic situation where the initial values of norm are random.
model_2a = simulate(ThoughtlessModel(N_AGENT, seed=SEED), N_STEP)
plot_results(model_2a);
We see that there are two disjoint regions where norm is True and two where it is False. (Recall that the agents are on a circle so the regions at the left and right ends of this plot are the same from their perspective.) "Individual computing at norm boundaries" corresponds to the fact that all agents sufficiently far in the interior of the constant-norm regions eventually have radius one, and it is only agents near the borders of these regions that "compute" in the sense of consulting more than their direct neighbors in their norm update decisions.
In order to understand the factors influencing the width of these edges that compute, we reduce the tolerance used when updating an agent's radius from 5% to 2.5%.
model_2b = simulate(ThoughtlessModel(N_AGENT, tol=0.025, seed=SEED), N_STEP)
plot_results(model_2b);
We see that a lower tolerance causes an agent to expand their radius more often (remember the agent checks a larger radius before considering a smaller radius) and therefore we end up with more agents with nontrivial radius. In the opposite direction, doubling the tolerance to 10% causes agents to increase their radius less often and therefore the boundary where computing occurs shrinks.
model_2c = simulate(ThoughtlessModel(N_AGENT, tol=0.1, seed=SEED), N_STEP)
plot_results(model_2c);
Run 4. Modest Noise Level and Endogenous Neighborhood Norms¶
For implementation reasons we will replicate run three in the paper last and skip to run four now. The first two scenarios treated a noiseless case where agents never spontaneously change their norm regardless of the input from their neighbors. This noiselessness is unrealistic, so we introduce a NoisyLazyStatistician subclass of LazyStatistician that randomly draws a new norm occasionally according to its noise parameter.
class NoisyLazyStatistician(LazyStatistician):
def __init__(self, unique_id, model, norm=None, radius=None, noise=0):
super().__init__(unique_id, model, norm=norm, radius=radius)
self.noise = noise
def update_norm(self, noise=None):
if self.random.random() < (self.noise if noise is None else noise):
self.norm = self.random.choice([False, True])
else:
super().update_norm()
We simulate behavior when agents randomize their norm due to noise 15% of the time and visualize the results below.
model_4 = simulate(
ThoughtlessModel(N_AGENT, seed=SEED, noise=0.15,
agent_cls=NoisyLazyStatistician),
N_STEP
)
plot_results(model_4);
These results are visually quite interesting but not hugely different from those in previous runs. For the most part we get two blocks each for each value of norm. Interestingly, in this case small endogenous dissenting islands occasionally apear inside these blocks and manage to persist for a few timesteps before disappearing.
Run 5. Higher Noise and Endogenous Neighborhood Norms¶
The next simulation doubles the noise level from 15% to 30%.
model_5 = simulate(
ThoughtlessModel(N_AGENT, seed=SEED, noise=0.3,
agent_cls=NoisyLazyStatistician),
N_STEP
)
plot_results(model_5);