Schelling Segregation Model

Thomas Schelling's 1971 model demonstrates how mild individual preferences for same-type neighbors can produce dramatic macro-level residential segregation -- a landmark result in complexity science showing emergent collective behaviour from simple micro-level rules.

1. The Model

Consider a two-dimensional square lattice of size $N \times N$ where each cell is occupied by an agent of type A, an agent of type B, or is empty. Each agent examines its Moore neighbourhood (8 surrounding cells) and counts the fraction of occupied neighbours that are of the same type:

$$f_i = \frac{\text{number of same-type neighbours of } i}{\text{total occupied neighbours of } i}$$

An agent is satisfied if $f_i \geq T$ where $T$ is the tolerance threshold. The dynamics proceed as follows:

Select a random unsatisfied agent.
Move that agent to a random empty cell.
Repeat until all agents are satisfied or a maximum number of iterations is reached.

The key insight is that even a modest threshold -- say $T = 1/3$, meaning an agent only requires one-third of its neighbours to be similar -- leads to near-complete segregation of the grid. Nobody individually demands a segregated neighbourhood, yet the collective outcome is stark separation.

2. Entropy-Based Segregation Index

To quantify how segregated the lattice has become, we use the Shannon entropy of the local neighbourhood compositions. For each occupied cell $i$, define the local entropy:

$$H_i = -\sum_{k \in \{A,B\}} p_k^{(i)} \ln p_k^{(i)}$$

where $p_k^{(i)}$ is the fraction of occupied neighbours of cell $i$ belonging to type $k$. The maximum possible entropy for two types is:

$$H_{\max} = -2 \cdot \tfrac{1}{2} \ln \tfrac{1}{2} = \ln 2$$

We average over all occupied cells to get the global entropy $\bar{H} = \frac{1}{M}\sum_i H_i$ where $M$ is the number of occupied cells. The normalised segregation index is:

$$S = 1 - \frac{\bar{H}}{H_{\max}} = 1 - \frac{\bar{H}}{\ln 2}$$

When $S = 0$, neighbourhoods are perfectly mixed (maximum entropy); when $S = 1$, every neighbourhood is completely homogeneous (zero entropy).

Derivation: Why Mild Preferences Segregate

Consider a well-mixed initial state. Each agent has roughly 50% same-type neighbours, so $f_i \approx 0.5$. If the threshold $T = 1/3$, most agents are already satisfied. But the few who are not (those near the boundary of a cluster of the other type) will move.

When an unsatisfied agent vacates a cell, it reduces the same-type fraction for remaining neighbours of the same type. This can cause a cascade: one departure makes several previously-satisfied agents unsatisfied. Mathematically, if agent $i$ leaves neighbourhood $\mathcal{N}_j$ of same-type agent $j$:

$$f_j \to f_j' = \frac{n_{\text{same}} - 1}{n_{\text{total}} - 1}$$

For a typical well-mixed neighbourhood with 4 same-type out of 8 total ($f_j = 0.5$), losing one same-type neighbour gives $f_j' = 3/7 \approx 0.43$. Still above $T=1/3$. But if the departing agent is replaced by an opposite-type agent, $f_j' = 3/8 = 0.375$. A second such swap gives $f_j' = 2/8 = 0.25 < 1/3$ -- agent $j$ now becomes unsatisfied and must move, triggering further cascades. This positive feedback drives the system toward segregation.

3. Phase Transition

The Schelling model exhibits a phase transition as the threshold $T$ increases. Below a critical threshold $T_c$, the system remains integrated; above it, the system rapidly segregates.

For a square lattice with Moore neighbourhoods and equal proportions of two types with 10% vacancy, numerical experiments show $T_c \approx 0.30$--$0.35$. The transition is characterised by a rapid increase in the segregation index $S$:

$$S(T) \approx \begin{cases} 0 & T < T_c \\ 1 - \exp\!\bigl[-\alpha(T - T_c)^{\beta}\bigr] & T \geq T_c \end{cases}$$

This resembles a continuous (second-order) phase transition. Near the critical point, the correlation length of same-type clusters diverges. We can define the cluster correlation function:

$$C(r) = \langle \delta_{\sigma_i, \sigma_j} \rangle_{|r_i - r_j| = r} - \langle \delta_{\sigma_i, \sigma_j} \rangle_{\text{random}}$$

where $\sigma_i$ is the type of agent at cell $i$ and the average is over all pairs at distance $r$. At $T_c$, this decays as a power law $C(r) \sim r^{-\eta}$ rather than exponentially, a hallmark of criticality.

4. Python Simulation: Schelling on a Grid

The following simulation implements the full Schelling model on a 50x50 grid. It tracks the segregation index over time for a given threshold and visualises the final state of the grid.

Schelling Segregation: Full Grid Simulation

Python

script.py134 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Parameters ---
N = 50              # grid size
frac_A = 0.45       # fraction type A
frac_B = 0.45       # fraction type B (rest is empty)
threshold = 0.40    # tolerance threshold
max_steps = 5000

# --- Initialise grid ---
# 0 = empty, 1 = type A, 2 = type B
total = N * N
n_A = int(frac_A * total)
n_B = int(frac_B * total)
n_empty = total - n_A - n_B

cells = np.array([1]*n_A + [2]*n_B + [0]*n_empty)
np.random.shuffle(cells)
grid = cells.reshape((N, N))

def get_neighbours(grid, i, j, N):
    """Return array of neighbour values (Moore neighbourhood)."""
    neighbours = []
    for di in [-1, 0, 1]:
        for dj in [-1, 0, 1]:
            if di == 0 and dj == 0:
                continue
            ni, nj = (i + di) % N, (j + dj) % N
            neighbours.append(grid[ni, nj])
    return np.array(neighbours)

def compute_segregation(grid, N):
    """Compute normalised entropy-based segregation index."""
    H_total = 0.0
    count = 0
    for i in range(N):
        for j in range(N):
            if grid[i, j] == 0:
                continue
            nbrs = get_neighbours(grid, i, j, N)
            occupied = nbrs[nbrs > 0]
            if len(occupied) == 0:
                continue
            n_same = np.sum(occupied == grid[i, j])
            n_diff = len(occupied) - n_same
            p_same = n_same / len(occupied)
            p_diff = n_diff / len(occupied)
            H = 0.0
            if p_same > 0:
                H -= p_same * np.log(p_same)
            if p_diff > 0:
                H -= p_diff * np.log(p_diff)
            H_total += H
            count += 1
    H_bar = H_total / count if count > 0 else 0
    S = 1.0 - H_bar / np.log(2)
    return S

def find_unsatisfied(grid, N, threshold):
    """Return list of (i,j) for unsatisfied agents."""
    unsatisfied = []
    for i in range(N):
        for j in range(N):
            if grid[i, j] == 0:
                continue
            nbrs = get_neighbours(grid, i, j, N)
            occupied = nbrs[nbrs > 0]
            if len(occupied) == 0:
                continue
            frac_same = np.sum(occupied == grid[i, j]) / len(occupied)
            if frac_same < threshold:
                unsatisfied.append((i, j))
    return unsatisfied

# --- Run simulation ---
seg_history = [compute_segregation(grid, N)]
steps_history = [0]

for step in range(1, max_steps + 1):
    unsatisfied = find_unsatisfied(grid, N, threshold)
    if len(unsatisfied) == 0:
        break
    # Pick a random unsatisfied agent
    idx = np.random.randint(len(unsatisfied))
    ai, aj = unsatisfied[idx]
    agent_type = grid[ai, aj]
    # Find empty cells
    empties = list(zip(*np.where(grid == 0)))
    if len(empties) == 0:
        break
    # Move to random empty cell
    eidx = np.random.randint(len(empties))
    ei, ej = empties[eidx]
    grid[ei, ej] = agent_type
    grid[ai, aj] = 0
    # Record segregation periodically
    if step % 200 == 0 or step == 1:
        seg_history.append(compute_segregation(grid, N))
        steps_history.append(step)

seg_history.append(compute_segregation(grid, N))
steps_history.append(step)

# --- Plot ---
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Grid visualisation
cmap = plt.cm.colors.ListedColormap(['#1a1a2e', '#10b981', '#f59e0b'])
axes[0].imshow(grid, cmap=cmap, interpolation='nearest')
axes[0].set_title(f'Final Grid State (T={threshold})', fontsize=13, color='white')
axes[0].set_xlabel('Column', color='white')
axes[0].set_ylabel('Row', color='white')
axes[0].tick_params(colors='white')

# Segregation over time
axes[1].plot(steps_history, seg_history, 'o-', color='#10b981', linewidth=2, markersize=4)
axes[1].set_xlabel('Simulation Step', fontsize=12, color='white')
axes[1].set_ylabel('Segregation Index S', fontsize=12, color='white')
axes[1].set_title('Segregation Index Over Time', fontsize=13, color='white')
axes[1].set_ylim(0, 1)
axes[1].tick_params(colors='white')
axes[1].grid(True, alpha=0.3)

fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()
print(f"Final segregation index: S = {seg_history[-1]:.4f}")
print(f"Converged after {step} steps, {len(find_unsatisfied(grid, N, threshold))} unsatisfied agents remain")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Phase Transition: Entropy vs Threshold

We sweep the threshold $T$ from 0 to 0.75 and measure the final segregation index for each value. The result reveals a sharp transition around $T_c \approx 0.30$.

Phase Transition: Segregation Index vs Threshold

Python

script.py92 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

N = 20
frac_A, frac_B = 0.45, 0.45
max_steps = 1000
thresholds = np.arange(0.0, 0.76, 0.10)

def run_schelling(N, frac_A, frac_B, threshold, max_steps):
    total = N * N
    n_A = int(frac_A * total)
    n_B = int(frac_B * total)
    cells = np.array([1]*n_A + [2]*n_B + [0]*(total - n_A - n_B))
    np.random.shuffle(cells)
    grid = cells.reshape((N, N))

for step in range(max_steps):
        unsatisfied = []
        for i in range(N):
            for j in range(N):
                if grid[i, j] == 0:
                    continue
                nbrs = []
                for di in [-1,0,1]:
                    for dj in [-1,0,1]:
                        if di==0 and dj==0: continue
                        nbrs.append(grid[(i+di)%N, (j+dj)%N])
                occ = [v for v in nbrs if v > 0]
                if len(occ) == 0: continue
                if sum(1 for v in occ if v == grid[i,j]) / len(occ) < threshold:
                    unsatisfied.append((i, j))
        if not unsatisfied:
            break
        ai, aj = unsatisfied[np.random.randint(len(unsatisfied))]
        empties = list(zip(*np.where(grid == 0)))
        if not empties: break
        ei, ej = empties[np.random.randint(len(empties))]
        grid[ei, ej] = grid[ai, aj]
        grid[ai, aj] = 0

# Compute segregation index
    H_total, count = 0.0, 0
    for i in range(N):
        for j in range(N):
            if grid[i,j] == 0: continue
            nbrs = []
            for di in [-1,0,1]:
                for dj in [-1,0,1]:
                    if di==0 and dj==0: continue
                    nbrs.append(grid[(i+di)%N, (j+dj)%N])
            occ = [v for v in nbrs if v > 0]
            if len(occ) == 0: continue
            p = sum(1 for v in occ if v == grid[i,j]) / len(occ)
            q = 1 - p
            H = 0
            if p > 0: H -= p * np.log(p)
            if q > 0: H -= q * np.log(q)
            H_total += H
            count += 1
    return 1.0 - (H_total / count / np.log(2)) if count > 0 else 0

# Sweep
results = []
n_trials = 3
for T in thresholds:
    vals = [run_schelling(N, frac_A, frac_B, T, max_steps) for _ in range(n_trials)]
    results.append((T, np.mean(vals), np.std(vals)))

Ts = [r[0] for r in results]
Ss = [r[1] for r in results]
errs = [r[2] for r in results]

fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(Ts, Ss, yerr=errs, fmt='o-', color='#10b981', linewidth=2.5,
            capsize=4, markersize=7, markerfacecolor='#34d399')
ax.axvline(x=0.30, color='#f59e0b', linestyle='--', linewidth=1.5, label=r'$T_c \approx 0.30$')
ax.axvspan(0.28, 0.35, alpha=0.15, color='#f59e0b', label='Critical region')
ax.set_xlabel('Threshold T', fontsize=13, color='white')
ax.set_ylabel('Segregation Index S', fontsize=13, color='white')
ax.set_title('Phase Transition in Schelling Model', fontsize=14, color='white')
ax.legend(fontsize=11, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.set_ylim(-0.05, 1.05)
ax.tick_params(colors='white')
ax.grid(True, alpha=0.3)
fig.patch.set_facecolor('#0f172a')
ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()
print("Phase transition sweep complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Urban Implications and Extensions

The Schelling model has profound implications for urban policy:

No discriminatory intent needed: Even agents who are comfortable being a minority (threshold well below 0.5) collectively produce segregated patterns.
Hysteresis: Once segregation forms, it is very difficult to reverse. Lowering the threshold below $T_c$ does not spontaneously de-segregate the grid because satisfied agents have no reason to move.
Tipping points: The introduction of a small number of dissimilar agents into a homogeneous neighbourhood can trigger cascading departures -- the "white flight" phenomenon studied extensively in American urban sociology.
Multi-group extensions: With $K > 2$ groups, the maximum entropy becomes $H_{\max} = \ln K$ and the critical threshold generally decreases as $T_c \sim 1/K$.

Connection to Statistical Mechanics

The Schelling model can be mapped to a zero-temperature Ising model with external field. The agents' utility function plays the role of a local Hamiltonian, and the threshold acts as a coupling constant. The segregation transition is analogous to the ferromagnetic phase transition in the Ising model -- a connection we will explore further in the mean-field theory chapters.

← Part IV Overview Bid-Rent Gradient →