Information Theory

CS1302 Introduction to Computer Programming

As mentioned in previous lectures, the following two lists coin_flips and dice_rolls simulate the random coin flips and rollings of a dice:

# Do NOT modify any variables defined here because some tests rely on them
import random
import math
import matplotlib.pyplot as plt
import ipywidgets as widgets
%matplotlib widget
%reload_ext divewidgets

random.seed(0)  # for reproducible results.
num_trials = 200
coin_flips = ["H" if random.random() <= 1 / 2 else "T" for i in range(num_trials)]
dice_rolls = [random.randint(1, 6) for i in range(num_trials)]
print("coin flips: ", *coin_flips)
print("dice rolls: ", *dice_rolls)

What is the difference of the two random processes?
Can we say one process has more information content than the other?

In this Lab, we will use dictionaries to store their distributions and then compute their information content using information theory, which was introduced by Claude E. Shannon. It has numerous applications:

compression (to keep files small)
communications (to send data to mobile phones), and
machine learning (to identify relevant features to learn from).

Entropy¶

Mathematically, we denote a distribution as $\mathbf{p}=[p_i]_{i\in \mathcal{S}}$ , where

$\mathcal{S}$ is the set of distinct outcomes, and
$p_i$ denotes the probability (chance) of seeing outcome $i$ .

The following code shown in the lecture uses a dictionary to store the distribution for a sequence efficiently without storing outcomes with zero counts:

# Do NOT modify any variables defined here because some tests rely on them
def distribute(seq):
    """Returns a dictionary where each value in a key-value pair is
    the probability of the associated key occurring in the sequence.
    """
    p = {}
    for i in seq:
        p[i] = p.get(i, 0) + 1 / len(seq)
    return p


# tests
coin_flips_dist = distribute(coin_flips)
dice_rolls_dist = distribute(dice_rolls)
print("Distribution of coin flips:", coin_flips_dist)
print("Distribution of dice rolls:", dice_rolls_dist)

For $\mathbf{p}$ to be a valid distribution, the probabilities $p_i$ 's have to sum to $1$ , i.e.,

\sum_{i\in \mathcal{S}} p_i = 1,

(1)

which can be verified as follows:

assert math.isclose(sum(coin_flips_dist.values()), 1) and math.isclose(
    sum(dice_rolls_dist.values()), 1
)

How to measure the information content?

Entropy

For instance, if $\mathbf{p}=(p_{H},p_{T})=(0.5,0.5)$ , then

\begin{aligned} H(\mathbf{p}) &= 0.5 \log_2 \frac{1}{0.5} + 0.5 \log_2 \frac{1}{0.5} \\ &= 0.5 + 0.5 = 1 \text{ bit,}\end{aligned}

(3)

i.e., an outcome of a fair coin flip has one bit of information content, as expected.

On the other hand, if $\mathbf{p}=(p_{H},p_{T})=(1,0)$ , then

\begin{aligned} H(\mathbf{p}) &= 1 \log_2 \frac{1}{1} + 0 \log_2 \frac{1}{0} \\ &= 0 + 0 = 0 \text{ bits,}\end{aligned}

(4)

i.e., an outcome of a biased coin flip that always comes up head has no information content, again as expected.

def entropy(dist):
    # YOUR CODE HERE
    raise NotImplementedError()

# tests
assert math.isclose(entropy({"H": 0.5, "T": 0.5}), 1)
assert math.isclose(entropy({"H": 1, "T": 0}), 0)
assert math.isclose(entropy(dict.fromkeys(range(1, 7), 1 / 6)), math.log2(6))

# HIDDEN

Uniform distribution maximizes entropy¶

Intuitively,

for large enough numbers of fair coin flips, we should have $\mathcal{S}=\{H,T\}$ and $p_H=p_T=0.5$ , i.e., equal chance for head and tail.
for large enough numbers of fair dice rolls, we should have $p_i=\frac16$ for all $i\in \mathcal{S}=\{1,2,3,4,5,6\}$ .

def plot_distribution(seq, *args, **kwargs):
    fig = plt.figure(*args, **kwargs)
    plt.xlabel("Outcomes")
    plt.title("Distribution")
    plt.ylim(0, 1)
    dist_plot = ()

    @widgets.interact(
        n=widgets.IntSlider(
            value=1,
            min=1,
            max=len(seq),
            step=1,
            description="n:",
            continuous_update=False,
        )
    )
    def plot(n):
        plt.figure(fig.number)
        nonlocal dist_plot
        dist_plot and dist_plot.remove()
        dist = distribute(seq[:n])
        dist_plot = plt.stem(dist.keys(), dist.values())
        plt.show()

plot_distribution(coin_flips, num=1, clear=True)

plot_distribution(dice_rolls, num=2, clear=True)

def uniform(outcomes):
    """Returns the uniform distribution (dict) over distinct items in outcomes."""
    # YOUR CODE HERE
    raise NotImplementedError()

# tests
assert uniform("HT") == {"H": 0.5, "T": 0.5}
assert uniform("HTH") == {"H": 0.5, "T": 0.5}
fair_dice_dist = uniform(range(1, 7))
assert all(
    math.isclose(fair_dice_dist[k], v)
    for k, v in {
        1: 0.16666666666666666,
        2: 0.16666666666666666,
        3: 0.16666666666666666,
        4: 0.16666666666666666,
        5: 0.16666666666666666,
        6: 0.16666666666666666,
    }.items())

# HIDDEN

What is the entropy for uniform distributions?

By definition,

\begin{aligned} H(\mathbf{p}) &:= \sum_{i\in \mathcal{S}} p_i \log_2 \tfrac{1}{p_i} \\ &= \sum_{i\in \mathcal{S}} \frac{1}{|\mathcal{S}|} \log_2 |\mathcal{S}| = \log_2 |\mathcal{S}| \kern1em \text{(bits)} \end{aligned}

(6)

Entropy reduces to the formula in Lecture 1 and Lab 1 regarding the number of bits required to encode a set of integers or characters. It is the maximum possible entropy for a given finite set of possible outcomes.

Use this result to test whether you have implemented both entropy and uniform correctly:

assert all(
    math.isclose(entropy(uniform(range(n))), math.log2(n)) for n in range(1, 100)
)

Joint distribution and its entropy¶

If we duplicate a sequence of outcomes, the total information content should remain unchanged, NOT doubled, because the duplicate contains the same information as the original. We will verify this fact by creating a joint distribution

\mathbf{p}=[p_{ij}]_{i\in \mathcal{S},j\in \mathcal{T}}

(7)

where $\mathcal{S}$ and $\mathcal{T}$ are sets of outcomes; and
$p_{ij}$ is the chance we see outcomes $i$ and $j$ simultaneously.

The subscript $ij$ in $p_{ij}$ denotes a tuple $(i,j)$ , NOT the multiplication $i\times j$ . We also have

\sum_{i\in \mathcal{S}} \sum_{j\in \mathcal{T}} p_{ij} = 1.

(8)

def jointly_distribute(seq1, seq2):
    """Returns the joint distribution of the tuple (i,j) of outcomes from zip(seq1,seq2)."""
    # YOUR CODE HERE
    raise NotImplementedError()

# tests
assert jointly_distribute("HT", "HT") == {("H", "H"): 0.5, ("T", "T"): 0.5}
assert jointly_distribute("HHTT", "HTHT") == {
    ("H", "H"): 0.25,
    ("H", "T"): 0.25,
    ("T", "H"): 0.25,
    ("T", "T"): 0.25}
coin_flips_duplicate_dist = {
    ("T", "T"): 0.5350000000000004,
    ("H", "H"): 0.4650000000000003}
coin_flips_duplicate_ans = jointly_distribute(coin_flips, coin_flips)
assert all(
    math.isclose(coin_flips_duplicate_ans[i], pi)
    for i, pi in coin_flips_duplicate_dist.items())

# HIDDEN

If you have implemented entropy and jointly_distribute correctly, you can verify that duplicating a sequence will give the same entropy.

assert math.isclose(
    entropy(jointly_distribute(coin_flips, coin_flips)), entropy(distribute(coin_flips))
)
assert math.isclose(
    entropy(jointly_distribute(dice_rolls, dice_rolls)), entropy(distribute(dice_rolls))
)

However, for two sequences generated independently, the joint entropy is roughly the sum of the individual entropies.

coin_flips_entropy = entropy(distribute(coin_flips))
dice_rolls_entropy = entropy(distribute(dice_rolls))
cf_dr_entropy = entropy(jointly_distribute(coin_flips, dice_rolls))
print(
    f"""Entropy of coin flip: {coin_flips_entropy}
Entropy of dice roll: {dice_rolls_entropy}
Sum of the above entropies: {coin_flips_entropy + dice_rolls_entropy}
Joint entropy: {cf_dr_entropy}"""
)

Conditional distribution and entropy¶

Mathematically, we denote a conditional distribution as $\mathbf{q}:=[q_{j|i}]_{i\in \mathcal{S}, j\in \mathcal{T}}$ , where

$\mathcal{S}$ and $\mathcal{T}$ are two sets of distinct outcomes, and
$q_{j|i}$ denotes the probability (chance) of seeing outcome $j$ given the condition that outcome $i$ is observed.

For $\mathbf{q}$ to be a valid distribution, the probabilities $q_{j|i}$ 's have to sum to $1$ for every $i$ , i.e.,

\sum_{j\in \mathcal{T}} q_{j|i} = 1 \kern1em \text{for all }i\in \mathcal{S}

(9)

For example, suppose we want to compute the distribution of coin flips given dice rolls, then the following assign q_H_1 and q_T_1 to the values $q_{H|1}$ and $q_{T|1}$ respectively:

coin_flips_1 = [j for i, j in zip(dice_rolls, coin_flips) if i == 1]
q_H_1 = coin_flips_1.count("H") / len(coin_flips_1)
q_T_1 = coin_flips_1.count("T") / len(coin_flips_1)
print("Coin flips given dice roll is 1:", coin_flips_1)
print(
    'Distribution of coin flip given dice roll is 1: {{ "H": {}, "T": {}}}'.format(
        q_H_1, q_T_1
    )
)
assert math.isclose(q_H_1 + q_T_1, 1)

Note that q_H_1 + q_T_1 is 1 as expected. Similarly, we can assign q_H_2 and q_T_2 to the values $q_{H|2}$ and $q_{T|2}$ respectively.

coin_flips_2 = [j for i, j in zip(dice_rolls, coin_flips) if i == 2]
q_H_2 = coin_flips_2.count("H") / len(coin_flips_2)
q_T_2 = coin_flips_2.count("T") / len(coin_flips_2)
print("Coin flips given dice roll is 2:", coin_flips_2)
print(
    'Distribution of coin flip given dice roll is 2: {{ "H": {}, "T": {}}}'.format(
        q_H_2, q_T_2
    )
)
assert math.isclose(q_H_2 + q_T_2, 1)

Finally, we want to store the conditional distribution as a nested dictionary so that q[i] points to the distribution

[q_{j|i}]_{j\in \mathcal{T}} \kern1em \text{for }i\in \mathcal{S}.

(10)

q = {}
q[1] = dict(zip("HT", (q_H_1, q_T_1)))
q[2] = dict(zip("HT", (q_H_2, q_T_2)))
q

Of course, the above dictionary is missing the entries for other possible outcomes of the dice rolls.

Exercise 4

Define a function conditionally_distribute that

takes two sequences seq1 and seq2 of outcomes of the same length, and
returns the conditional distribution of seq2 given seq1 in the form of a nested dictionary efficiently without storing the unobserved outcomes.

In the above example, seq1 is dice_rolls while seq2 is coin_flips.

Hint

For an efficient implementation without traversing the input sequences too many times, consider using the following solution template and the setdefault method.

def conditionally_distribute(seq1, seq2):
    q, count = {}, {}  # NOT q = count = {}
    for i in seq1:
        count[i] = count.get(i, 0) + 1
    for i, j in zip(seq1, seq2):
        q[i][j] = ...
    return q

def conditionally_distribute(seq1, seq2):
    """Returns the conditional distribution q of seq2 given seq1 such that
    q[i] is a dictionary for observed outcome i in seq1 and
    q[i][j] is the probability of observing j in seq2 given the
    corresponding outcome in seq1 is i."""
    # YOUR CODE HERE
    raise NotImplementedError()

# tests
cf_given_dr_dist = {
    4: {"T": 0.5588235294117647, "H": 0.4411764705882353},
    1: {"T": 0.46511627906976744, "H": 0.5348837209302325},
    3: {"H": 0.5135135135135135, "T": 0.4864864864864865},
    6: {"H": 0.5454545454545454, "T": 0.45454545454545453},
    2: {"T": 0.7586206896551724, "H": 0.2413793103448276},
    5: {"T": 0.5416666666666666, "H": 0.4583333333333333}}
cf_given_dr_ans = conditionally_distribute(dice_rolls, coin_flips)
assert all(
    math.isclose(cf_given_dr_ans[i][j], v)
    for i, d in cf_given_dr_dist.items()
    for j, v in d.items())

# HIDDEN

# YOUR CODE HERE
raise NotImplementedError()

# tests
cf_given_dr_dist = {
    4: {"T": 0.5588235294117647, "H": 0.4411764705882353},
    1: {"T": 0.46511627906976744, "H": 0.5348837209302325},
    3: {"H": 0.5135135135135135, "T": 0.4864864864864865},
    6: {"H": 0.5454545454545454, "T": 0.45454545454545453},
    2: {"T": 0.7586206896551724, "H": 0.2413793103448276},
    5: {"T": 0.5416666666666666, "H": 0.4583333333333333}}
assert (
    conditional_entropy(
        {"H": 0.5, "T": 0.5}, {"H": {"H": 0.5, "T": 0.5}, "T": {"H": 0.5, "T": 0.5}})
    == 1)
assert (
    conditional_entropy(
        {"H": 0, "T": 1}, {"H": {"H": 0.5, "T": 0.5}, "T": {"H": 0.5, "T": 0.5}})
    == 1)
assert (
    conditional_entropy(
        {"H": 0.5, "T": 0.5}, {"H": {"H": 1, "T": 0}, "T": {"H": 0, "T": 1}})
    == 0)
assert (
    conditional_entropy(
        {"H": 0.5, "T": 0.5}, {"H": {"H": 1, "T": 0}, "T": {"H": 0.5, "T": 0.5}})
    == 0.5)
assert math.isclose(
    conditional_entropy(dice_rolls_dist, cf_given_dr_dist), 0.9664712793722372)

# HIDDEN

The joint probability $p_{ij}$ over $i\in \mathcal{S}$ and $j\in \mathcal{T}$ can be calculated as follows

p_{ij} = p_{i} q_{j|i}

(12)

where $p_i$ is the probability of $i$ and $q_{j|i}$ is the conditional probability of $j$ given $i$ .

def joint_distribution(p, q):
    # YOUR CODE HERE
    raise NotImplementedError()

# tests
assert joint_distribution(
    {"H": 0.5, "T": 0.5}, {"H": {"H": 0.5, "T": 0.5}, "T": {"H": 0.5, "T": 0.5}}) == {("H", "H"): 0.25, ("H", "T"): 0.25, ("T", "H"): 0.25, ("T", "T"): 0.25}
assert joint_distribution(
    {"H": 0, "T": 1}, {"H": {"H": 0.5, "T": 0.5}, "T": {"H": 0.5, "T": 0.5}}) == {("H", "H"): 0.0, ("H", "T"): 0.0, ("T", "H"): 0.5, ("T", "T"): 0.5}
assert joint_distribution(
    {"H": 0.5, "T": 0.5}, {"H": {"H": 1, "T": 0}, "T": {"H": 0, "T": 1}}) == {("H", "H"): 0.5, ("H", "T"): 0.0, ("T", "H"): 0.0, ("T", "T"): 0.5}, {
    "H": 0.5,
    "T": 0.5}

# HIDDEN

Finally, a fundamental information identity relating the joint and conditional entropies is the chain rule:

If you have implemented jointly_distribute, conditionally_distribute, entropy, and conditional_entropy correctly, we can verify the identity as follows.

def validate_chain_rule(seq1, seq2):
    p = distribute(seq1)
    q = conditionally_distribute(seq1, seq2)
    pq = jointly_distribute(seq1, seq2)
    H_pq = entropy(pq)
    H_p = entropy(p)
    H_q_p = conditional_entropy(p, q)
    print(
        f"""Entropy of seq1: {H_p}
Conditional entropy of seq2 given seq1: {H_q_p}
Sum of the above entropies: {H_p + H_q_p}
Joint entropy: {H_pq}"""
    )
    assert math.isclose(H_pq, H_p + H_q_p)

validate_chain_rule(coin_flips, coin_flips)

validate_chain_rule(dice_rolls, coin_flips)