Scale-Free Networks

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.

# Install empiricaldist if needed (Colab doesn't have it)

try:
    import empiricaldist
except ImportError:
    %pip install -q empiricaldist

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx

np.set_printoptions(legacy='1.25')

from empiricaldist import Pmf, Cdf

def decorate(**options):
    plt.gca().set(**options)
    plt.tight_layout()

Graphs

To represent social networks, we’ll use nx.Graph, the graph representation provided by NetworkX.

Each person is represented by a node. Each friendship is represented by an edge between two nodes.

Here’s a simple example with 4 people:

G = nx.Graph()
G.add_edge(1, 0)
G.add_edge(2, 0)
G.add_edge(3, 0)
nx.draw(G)

The number of friends a person has is the number of edges that connect to their node, which is the “degree” of the node.

for node in G.nodes():
    print(node, G.degree(node))

1 1
0 3
2 1
3 1

We are often interested in the “degree distribution” of a graph, which is the number of people who have 0 friends, the number who have 1 friend, and so on.

The following function extracts a list of degrees, one for each node in a graph.

def degrees(G):
    """List of degrees for nodes in `G`.
    
    G: Graph object
    
    returns: list of int
    """
    return [G.degree(node) for node in G]

Here’s the result for the small example.

degrees(G)

[1, 3, 1, 1]

I’ll use Pmf from empiricaldist to make a probability mass function.

pmf = Pmf.from_seq(degrees(G))
pmf

	probs
1	0.75
3	0.25

And bar to display it as a bar plot.

pmf.bar()
decorate(xlabel='Degree', ylabel='PMF', xlim=[0.4, 3.6])

Exercise: Add another node or nodes to the graph above, and add a few edges. Plot the new degree distribution.

# Solution

G.add_edge(4, 1)
G.add_edge(4, 2)
nx.draw(G)

# Solution

pmf = Pmf.from_seq(degrees(G))
pmf.bar()
decorate(xlabel='Degree', ylabel='PMF', xlim=[0.4, 3.6])

Facebook data

The following function reads a file with one edge per line, specified by two integer node IDs.

def read_graph(filename):
    """Read a graph from a file.
    
    filename: string
    
    return: Graph
    """
    G = nx.Graph()
    array = np.loadtxt(filename, dtype=int)
    G.add_edges_from(array)
    return G

We’ll read the Facebook data downloaded from SNAP.

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
        
download('https://github.com/AllenDowney/ComplexityScience/raw/master/facebook_combined.txt.gz')

# https://snap.stanford.edu/data/facebook_combined.txt.gz

fb = read_graph('facebook_combined.txt.gz')
n = len(fb)
m = len(fb.edges())
n, m

(4039, 88234)

To see how popular “you” are, on average, we’ll draw a random sample of 1000 people.

sample = np.random.choice(fb.nodes(), 1000, replace=True)

For each “you” in the sample, we’ll look up the number of friends.

sample_friends = [fb.degree(node) for node in sample]

To plot the degree distribution, I’ll use Pmf.

Pmf.from_seq(sample_friends).plot()
decorate(xlabel='Number of friends', 
         ylabel='PMF',
         title='Degree PMF, friends')

Cdf.from_seq(sample_friends).plot(color='C0')
decorate(xlabel='Number of friends', 
         ylabel='CDF',
         title='Degree CDF, friends')

Now what if, instead of “you”, we choose one of your friends, and look up the number of friends your friend has.

sample_fof = []
for node in sample:
    friends = list(fb.neighbors(node))
    friend = np.random.choice(friends)
    sample_fof.append(fb.degree(friend))

Here’s the degree distribution for your friend’s friends:

Pmf.from_seq(sample_fof).plot()
decorate(xlabel='Number of friends', 
         ylabel='PMF',
         title='Degree PMF, friends of friends')

Cdf.from_seq(sample_fof).plot(color='C1')
decorate(xlabel='Number of friends', 
         ylabel='CDF',
         title='Degree CDF, friends of friends')

Cdf.from_seq(sample_friends).plot(color='C0', label='friends')
Cdf.from_seq(sample_fof).plot(color='C1', label='friends of friends')
decorate(xlabel='Number of friends', 
         ylabel='CDF',
         title='Degree distributions')

The bulk of the distribution is wider, and the tail is thicker. This difference is reflected in the means:

np.mean(sample_friends), np.mean(sample_fof)

(46.642, 97.832)

And we can estimate the probability that your friend has more friends than you.

np.mean([friend > you for you, friend in zip(sample_friends, sample_fof)])

0.649

Power law distributions

As we’ll see below, the degree distribution in the Facebook data looks, in some ways, like a power law distribution. To see what that means, we’ll look at the Zipf distribution, which has a power law tail.

Here’s a sample from a Zipf distribution.

zipf_sample = np.random.zipf(a=2, size=10000)

Here’s what the PMF looks like.

pmf = Pmf.from_seq(zipf_sample)
pmf.plot()
decorate(xlabel='Zipf sample', ylabel='PMF')

Here it is on a log-x scale.

pmf.plot()
decorate(xlabel='Zipf sample (log)', ylabel='PMF', xscale='log')

And on a log-log scale.

pmf.plot()
decorate(xlabel='Zipf sample (log)', ylabel='PMF (log)', 
                 xscale='log', yscale='log')

On a log-log scale, the PMF of the Zipf distribution looks like a straight line (until you get to the extreme tail, which is discrete and noisy).

For comparison, let’s look at the Poisson distribution, which does not have a power law tail. I’ll choose the Poisson distribution with the same mean as the sample from the Zipf distribution.

mu, sigma = zipf_sample.mean(), zipf_sample.std()
mu, sigma

(6.0111, 76.52610519809564)

poisson_sample = np.random.poisson(lam=mu, size=10000)
poisson_sample.mean(), poisson_sample.std()

(6.0052, 2.4442121348197254)

Here’s the PMF on a log-log scale. It is definitely not a straight line.

poisson_pmf = Pmf.from_seq(poisson_sample)
poisson_pmf.plot()
decorate(xlabel='Poisson sample (log)', ylabel='PMF (log)', 
                 xscale='log', yscale='log')

So this gives us a simple way to test for power laws. If you plot the PMF on a log-log scale, and the result is a straight line, that is evidence of power law behavior.

This test is not entirely reliable; there are better options. But it’s good enough for an initial exploration.

Barabási and Albert

Let’s see what the degree distribution for the Facebook data looks like on a log-log scale.

pmf_fb = Pmf.from_seq(degrees(fb))

pmf_fb.plot(label='Facebook')
decorate(xscale='log', 
         yscale='log',
         xlabel='Degree (log)', 
         ylabel='PMF (log)')

For degrees greater than 10, it resembles the Zipf sample (and doesn’t look much like the Poisson sample).

We can estimate the parameter of the Zipf distribution by eyeballing the slope of the tail.

plt.plot([10, 1000], [5e-2, 2e-4], color='gray', linestyle='dashed')

pmf_fb.plot(label='Facebook')
decorate(xscale='log', 
         yscale='log',
         xlabel='Degree (log)', 
         ylabel='PMF (log)')

Here’s a simplified version of the NetworkX function that generates BA graphs.

# modified version of the NetworkX implementation from
# https://github.com/networkx/networkx/blob/main/networkx/generators/random_graphs.py

import random

def barabasi_albert_graph(n, k, seed=None):
    """Constructs a BA graph.
    
    n: number of nodes
    k: number of edges for each new node
    seed: random seed
    """
    if seed is not None:
        random.seed(seed)
    
    G = nx.empty_graph(k)
    targets = set(range(k))
    repeated_nodes = []

    for source in range(k, n):

        G.add_edges_from(zip([source]*k, targets))

        repeated_nodes.extend(targets)
        repeated_nodes.extend([source] * k)

        targets = _random_subset(repeated_nodes, k)

    return G

And here’s the function that generates a random subset without repetition.

def _random_subset(repeated_nodes, k):
    """Select a random subset of nodes without repeating.
    
    repeated_nodes: list of nodes
    k: size of set
    
    returns: set of nodes
    """
    targets = set()
    while len(targets) < k:
        x = random.choice(repeated_nodes)
        targets.add(x)
    return targets

I’ll generate a BA graph with the same number of nodes and edges as the Facebook data:

n = len(fb)
m = len(fb.edges())
k = int(round(m/n))
n, m, k

(4039, 88234, 22)

Providing a random seed means we’ll get the same graph every time.

ba = barabasi_albert_graph(n, k, seed=15)

The number of edges is pretty close to what we asked for.

len(ba), len(ba.edges()), len(ba.edges())/len(ba)

(4039, 88374, 21.88016835850458)

So the mean degree is about right.

np.mean(degrees(fb)), np.mean(degrees(ba))

(43.69101262688784, 43.76033671700916)

The standard deviation of degree is pretty close; maybe a little low.

np.std(degrees(fb)), np.std(degrees(ba))

(52.41411556737521, 41.03760075705614)

Let’s take a look at the degree distribution.

pmf_ba = Pmf.from_seq(degrees(ba))

Looking at the PMFs on a linear scale, we see one difference, which is that the BA model has no nodes with degree less than k, which is 22.

plt.figure(figsize=(12,6))
plt.subplot(1, 2, 1)

pmf_fb.plot(label='Facebook')
decorate(xlabel='Degree', ylabel='PMF')

plt.subplot(1, 2, 2)

pmf_ba.plot(label='BA model')
decorate(xlabel='Degree', ylabel='PMF')

If we look at the PMF on a log-log scale, the BA model looks pretty good for values bigger than about 20. And it seems to follow a power law.

plt.figure(figsize=(12,6))
plt.subplot(1, 2, 1)

pmf_fb.plot(label='Facebook')
decorate(xlabel='Degree', ylabel='PMF',
         xscale='log', yscale='log')

plt.subplot(1, 2, 2)

pmf_ba.plot(label='BA model')
decorate(xlabel='Degree', ylabel='PMF',
         xlim=[1, 1e4],
         xscale='log', yscale='log')

Cumulative distributions

Here are the degree CDFs for the Facebook data and the BA model.

cdf_fb = Cdf.from_seq(degrees(fb))

cdf_ba = Cdf.from_seq(degrees(ba))

If we plot them on a log-x scale, we get a sense of how well the model fits the central part of the distribution.

The BA model is ok for values above the median, but not very good for smaller values.

cdf_fb.plot(label='Facebook')
cdf_ba.plot(color='gray', label='BA model')
decorate(xlabel='Degree', xscale='log',
                 ylabel='CDF')

If we plot the complementary CDF on a log-log scale, we see that the BA model fits the tail of the distribution reasonably well.

complementary_cdf_fb = 1-cdf_fb
complementary_cdf_ba = 1-cdf_ba

complementary_cdf_fb.plot(label='Facebook')
complementary_cdf_ba.plot(color='gray', label='BA model')
decorate(xlabel='Degree', xscale='log',
         ylabel='Complementary CDF', yscale='log')

But there is certainly room for a model that does a better job of fitting the whole distribution.

Copyright 2016–2026 Allen Downey

Code license: MIT License

Text license: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International