04 - Noisy channels and capacity: from error probability to Shannon's second theorem

Now and then one bit in ten arrives flipped: it happens on an ageing disk, on a noisy radio link, on a Wi-Fi connection far from the router. You want to deliver the data without errors, and the obvious fix is to repeat every bit three times and take a majority vote. It works (the error drops to 2.8%), but the speed falls to a third. The question this article answers is: how fast can you communicate reliably over a channel that errs with probability p? Shannon's answer is a number, the capacity C = 1 - H2(p), which for p = 0.1 is 0.531 bits per channel use. Below that rate there are codes with an error as small as you like, and repetition is very far from the best; above it, no code can have an arbitrarily small error.

The theorem is, however, a statement of existence: it does not hand you the code, and it does not say how long the code must be. Let us look at the precise statement, compute the capacity in three independent ways, and simulate a channel to watch the numbers move.

What a channel is and what capacity measures

A discrete memoryless channel is a table Q(y|x): for each input symbol x, the probability of receiving each y, independently of what happened in earlier uses. The binary symmetric channel (BSC) has inputs and outputs in {0, 1} and flips the bit with probability p, the same in both directions. A block code (N, K) is a list of 2K codewords of length N; its rate is R = K/N bits per channel use. With these definitions, which follow MacKay (§9.6), capacity is defined as the maximum of the mutual information over the input distribution: C = maxP(x) I(X;Y).

The definition looks abstract, but it has a concrete meaning: the mutual information I(X;Y) = H(Y) - H(Y|X) is how much uncertainty about the input disappears when you see the output. For the BSC, if the input is 1 with probability q, the output is 1 with probability q(1-p) + (1-q)p, while H(Y|X) = H2(p) does not depend on q, because the noise has the same strength whichever bit you send. So I = H2(q(1-p) + (1-q)p) - H2(p), and the first term reaches its maximum, 1 bit, when the output is equiprobable, that is for q = 1/2. It follows that C = 1 - H2(p), MacKay's formula (1.35).

Three independent ways to compute the same capacity

A formula I have just derived by hand is not yet a verified formula. The listing below checks it in three ways that do not share the reasoning: the closed form with SymPy, the search for the maximum of I(X;Y) computed from the definition over 2001 values of q, and the Blahut-Arimoto algorithm, which finds the capacity of any channel by iterating on the input distribution without knowing anything about symmetries. The functions mutual_info and blahut_arimoto accept any matrix Q: that is why, in the second half, I also use them on two channels that are not symmetric.

import math
import numpy as np
import sympy as sp
from scipy.optimize import minimize_scalar


def h2(p):
    """Binary entropy in bits, with the convention 0 log 0 = 0."""
    if p in (0.0, 1.0):
        return 0.0
    return -p * math.log2(p) - (1 - p) * math.log2(1 - p)


def mutual_info(px, Q):
    """I(X;Y) in bits from the definition: sum of p(x) Q(y|x) log2( Q(y|x) / p(y) )."""
    py = px @ Q
    total = 0.0
    for x in range(Q.shape[0]):
        for y in range(Q.shape[1]):
            if px[x] > 0 and Q[x, y] > 0:
                total += px[x] * Q[x, y] * math.log2(Q[x, y] / py[y])
    return total


def blahut_arimoto(Q, iterations=3000):
    """Capacity of a discrete memoryless channel, by iteration (Blahut-Arimoto)."""
    n_x = Q.shape[0]
    r = np.full(n_x, 1.0 / n_x)
    for _ in range(iterations):
        q = r @ Q
        c = np.ones(n_x)
        for x in range(n_x):
            for y in range(Q.shape[1]):
                if Q[x, y] > 0:
                    c[x] *= (Q[x, y] / q[y]) ** Q[x, y]
        r = r * c / np.sum(r * c)
    return mutual_info(r, Q), r


def bsc(f):
    return np.array([[1 - f, f], [f, 1 - f]])


def bec(f):  # outputs: 0, ?, 1
    return np.array([[1 - f, f, 0.0], [0.0, f, 1 - f]])


def z_channel(f):  # 0 always arrives as 0; 1 arrives as 0 with probability f
    return np.array([[1.0, 0.0], [f, 1 - f]])


# 1) Closed form with exact fractions (SymPy)
p = sp.symbols("p", positive=True)
C = 1 + p * sp.log(p, 2) + (1 - p) * sp.log(1 - p, 2)
closed = sp.Rational(3, 4) * sp.log(3, 2) - 1  # H2(1/4) = 2 - (3/4) log2 3
print("p = 1/4 : C - ((3/4) log2 3 - 1) =", sp.simplify(C.subs(p, sp.Rational(1, 4)) - closed), "; value", sp.N(closed, 12))
print("p = 1/2 :", sp.simplify(C.subs(p, sp.Rational(1, 2))))
print("C(p) - C(1-p) at p = 1/10 :", sp.simplify(C.subs(p, sp.Rational(1, 10)) - C.subs(p, sp.Rational(9, 10))))

# 2) Three independent ways to the same number: formula, maximum on a grid, Blahut-Arimoto
print()
print("   f   formula  grid     Blahut-Arimoto")
for f in (0.01, 0.1, 0.15, 0.2, 0.5):
    formula = 1 - h2(f)
    Q = bsc(f)
    grid = max(mutual_info(np.array([1 - q, q]), Q) for q in np.linspace(0, 1, 2001))
    ba, _ = blahut_arimoto(Q)
    print(f"{f:5.2f}   {formula:.6f}  {grid:.6f}  {ba:.6f}")

# 3) Non-symmetric channels: erasure and Z channel, f = 0.15 (MacKay, example 9.11 and exercise 9.13)
f = 0.15
cap_bec, _ = blahut_arimoto(bec(f))
print()
print(f"BEC f=0.15: Blahut-Arimoto {cap_bec:.6f}   formula 1-f = {1 - f:.6f}")
cap_z, r = blahut_arimoto(z_channel(f))
# Second route: maximum of H2(p1 (1-f)) - p1 H2(f), MacKay eq. 9.14, with a scalar optimizer
res = minimize_scalar(lambda p1: -(h2(p1 * (1 - f)) - p1 * h2(f)), bounds=(0, 1), method="bounded",
                      options={"xatol": 1e-12})
print(f"Z   f=0.15: Blahut-Arimoto {cap_z:.6f} with P(X=1)={r[1]:.4f}   optimizer {-res.fun:.6f} with P(X=1)={res.x:.4f}")
print(f"            with uniform input: {mutual_info(np.array([0.5, 0.5]), z_channel(f)):.6f}")

# Output:
# p = 1/4 : C - ((3/4) log2 3 - 1) = 0 ; value 0.188721875541
# p = 1/2 : 0
# C(p) - C(1-p) at p = 1/10 : 0
#
#    f   formula  grid     Blahut-Arimoto
#  0.01   0.919207  0.919207  0.919207
#  0.10   0.531004  0.531004  0.531004
#  0.15   0.390160  0.390160  0.390160
#  0.20   0.278072  0.278072  0.278072
#  0.50   0.000000  0.000000  0.000000
#
# BEC f=0.15: Blahut-Arimoto 0.850000   formula 1-f = 0.850000
# Z   f=0.15: Blahut-Arimoto 0.685418 with P(X=1)=0.4449   optimizer 0.685418 with P(X=1)=0.4449
#             with uniform input: 0.678788

The first line of the output shows that the capacity at p = 1/4 is exactly (3/4) log2 3 - 1, about 0.1887, and the second that at p = 1/2 it is zero: a channel that flips half of the bits carries nothing. The third says that C(p) = C(1-p): a channel that is almost always wrong is as good as one that is almost always right, because you only have to invert what you receive. The table shows that the three columns agree to the sixth digit; the values 0.919, 0.531 and 0.278 are those for p equal to 0.01, 0.1 and 0.2, and p = 0.15 gives 0.390, MacKay's figure 0.39 (example 9.9). The listing uses only SymPy 1.14, NumPy 2.3 and SciPy 1.16 under Python 3.14, and it is deterministic.

Non-symmetric channels: erasure and the Z channel

In the second half of the output two channels with f = 0.15 appear. The binary erasure channel (BEC) never flips a bit: with probability f it replaces it with a symbol "?" that the receiver recognises. Its capacity is 1 - f = 0.85, and Blahut-Arimoto recovers it to six digits. The comparison with the BSC is instructive: a BSC with p = 0.15 reaches 0.390 bits, the BEC with the same probability 0.85, because in the BEC you know where the information was lost and you do not have to guess which bits are wrong.

The Z channel is not symmetric: a 0 always arrives as 0, a 1 arrives as 0 with probability f. Here the uniform input is not optimal: the mutual information is 0.6788 with P(X=1) = 0.5 and rises to 0.6854 with P(X=1) ≈ 0.445, a value that matches MacKay's (example 9.11, 0.685 with 0.445) and that two different procedures, Blahut-Arimoto and a scalar optimizer on formula (9.14), both find. The lesson is that "the uniform input is optimal" is a property of the BSC and the BEC, not of capacity: the maximization over the input is part of the definition and must not be skipped.

What Shannon's second theorem states

Shannon states it in section 13 of his 1948 paper (theorem 11); MacKay (§10.1) writes it in three parts for a discrete memoryless channel:

1. for every ε > 0 and every rate R < C, for N large enough there exists a code of length N and rate at least R, with a decoder for which the maximal block error probability is less than ε;
2. if you accept a bit error probability pb, rates up to R(pb) = C / (1 - H2(pb)) are achievable;
3. for any pb, rates greater than R(pb) are not achievable.

The first part is the guarantee; the third is the limit. The second is the part that gets forgotten and the most useful in practice: accepting a bit error of 1% moves the limit from 0.531 to 0.578 bits per use, and an error of 5% to 0.744. This does not contradict the converse: it is the theory of lossy compression that you meet in article 6 of the series, because allowing errors means throwing information away. From the last line of the listing's output you can also see that, with pb = p = 0.1, the limit is exactly 1: you can use the channel without a code, and you receive exactly the 10% of errors.

Just as important is what the theorem does not say:

• how large N must be for the error to fall below a given threshold, nor how N depends on the distance C - R;
• how to construct the code: Shannon's proof is by random choice, and section 14 of his paper observes that, apart from trivial cases, no explicit description is known of codes that approach the limit;
• that the error becomes zero at finite length: it becomes arbitrarily small, never null, for the BSC or the BEC;
• anything about delay and the complexity of encoding and decoding, which grow with N.

A practical consequence is that capacity serves to measure how far you are from the best possible, not to design. When a real code works at a stated rate and a certain error probability, the comparison with C tells you how much theoretical margin remains, and nothing more: codes built under complexity constraints are judged by their own parameters, which is the subject of the next article.

Why repetition is not enough: the price of reliability

The repetition code RN sends every bit N times and decides by majority; with N = 3 the bit error is 3p² - 2p³ (MacKay, eq. 1.36), that is 7/250 = 0.028 for p = 0.1. The listing checks the formula in three ways: the exact fraction, the binomial sum and a fixed-seed simulation of one million bits, and works out how much N has to grow to bring the error below a threshold. The simulation serves to compare the measured frequency with the theoretical one: the standard error of a frequency over n trials is √(p(1-p)/n), and the deviation must be read in those units.

import math
from fractions import Fraction
import numpy as np

SEED = 20261113
p = Fraction(1, 10)
C = 1 + 0.1 * math.log2(0.1) + 0.9 * math.log2(0.9)  # capacity of the BSC with p = 0.1


def pb_repetition(N, p):
    """Exact bit error probability of the repetition code R_N (N odd), majority decoding."""
    return sum(math.comb(N, k) * p**k * (1 - p) ** (N - k) for k in range((N + 1) // 2, N + 1))


# 1) The simulated channel: frequency of flipped bits against p
rng = np.random.default_rng(SEED)
n = 1_000_000
flip = rng.random(n) < 0.1
freq = flip.mean()
sigma = math.sqrt(0.1 * 0.9 / n)
print(f"flipped bits: {freq:.5f} out of {n}; theoretical 0.1; standard error {sigma:.5f}; deviation {abs(freq - 0.1) / sigma:.2f} sigma")

# 2) Repetition code R3: closed form 3p^2 - 2p^3 (MacKay, eq. 1.36), exact and simulated
print("R3 formula 3p^2-2p^3 =", 3 * p**2 - 2 * p**3, "=", float(3 * p**2 - 2 * p**3))
print("R3 binomial sum      =", pb_repetition(3, p))
bit = rng.integers(0, 2, n)
rx = np.stack([bit ^ (rng.random(n) < 0.1) for _ in range(3)])
decod = (rx.sum(axis=0) >= 2).astype(int)
sim = (decod != bit).mean()
sigma3 = math.sqrt(0.028 * 0.972 / n)
print(f"R3 simulated         = {sim:.5f}; standard error {sigma3:.5f}; deviation {abs(sim - 0.028) / sigma3:.2f} sigma")

# 3) What repetition costs: the smallest odd N with pb below the threshold, and its rate
for threshold in (Fraction(1, 10**3), Fraction(1, 10**6), Fraction(1, 10**9)):
    N = 1
    while pb_repetition(N, p) >= threshold:
        N += 2
    print(f"pb < {float(threshold):.0e}: N = {N:3d}, rate 1/N = {1 / N:.4f}, capacity {C:.4f}")

# Output:
# flipped bits: 0.09984 out of 1000000; theoretical 0.1; standard error 0.00030; deviation 0.52 sigma
# R3 formula 3p^2-2p^3 = 7/250 = 0.028
# R3 binomial sum      = 7/250
# R3 simulated         = 0.02800; standard error 0.00016; deviation 0.00 sigma
# pb < 1e-03: N =   9, rate 1/N = 0.1111, capacity 0.5310
# pb < 1e-06: N =  23, rate 1/N = 0.0435, capacity 0.5310
# pb < 1e-09: N =  35, rate 1/N = 0.0286, capacity 0.5310

The simulated channel flips 9.984% of the bits against the theoretical 10%, 0.52 standard errors away, hence compatible. The measured error of R3, 0.02800, agrees with 7/250 within the standard error (0.00016). The important part is the last one: to go below 10-9 repetition requires N = 35 and therefore a rate of 0.0286 bits per use, about nineteen times less than the capacity. With repetition, more reliable means slower; the second theorem says that this trade is not compulsory: the rate can stay close to 0.531 while the error falls, provided you accept much longer blocks and much less naive codes.

Why the number is 2^nC: counting the noises

I do not prove the theorem, which requires jointly typical sequences (MacKay §10.2-10.3), but there is a count that makes the formula natural. Over n transmitted bits, the noise flips about np positions, and the ways to choose them are binom(n, np), about 2nH2(p). For two different codewords not to be confused, the "clouds" of 2nH2(p) possible outputs around each codeword must not overlap in a systematic way, and the total outputs are 2n: the number of codewords cannot exceed about 2n / 2nH2(p) = 2n(1 - H2(p)). This is a plausibility argument, not a proof, and the listing checks its numerical part.

import math
import sympy as sp


def h2(p):
    return -p * math.log2(p) - (1 - p) * math.log2(1 - p)


p = 0.1
C = 1 - h2(p)

# 1) How many noise patterns there are: (1/n) log2 binom(n, pn) -> H2(p), and the sum up to pn has the same limit
print("   n    (1/n)log2 C(n,pn)   (1/n)log2 V(n,pn)   H2(0.1)")
for n in (10, 100, 1000, 10000):
    t = round(p * n)
    comb = math.comb(n, t)
    vol = sum(math.comb(n, j) for j in range(t + 1))
    print(f"{n:6d}   {math.log2(comb) / n:.6f}          {math.log2(vol) / n:.6f}           {h2(p):.6f}")

# 2) The Hamming bound with t = pn tends to 1 - H2(p): same value as the capacity
print()
print("   n    (1/n) log2( 2^n / V(n,pn) )   capacity 1 - H2(0.1)")
for n in (10, 100, 1000, 10000):
    t = round(p * n)
    vol = sum(math.comb(n, j) for j in range(t + 1))
    print(f"{n:6d}   {(n - math.log2(vol)) / n:.6f}                      {C:.6f}")

# 3) The weak converse: Fano + data processing give nR <= nC + 1 + Pe nR. Solved for Pe with SymPy.
n_, R_, C_, Pe = sp.symbols("n R C Pe", positive=True)
solution = sp.solve(sp.Eq(n_ * R_, n_ * C_ + 1 + Pe * n_ * R_), Pe)[0]
print()
print("Pe >=", sp.simplify(solution))
R = 0.8
for n in (100, 1000, 10000):
    print(f"R = 0.8, C = {C:.4f}, n = {n:5d}:  Pe >= {1 - C / R - 1 / (n * R):.4f}")

# 4) With an accepted bit error, the capacity can be exceeded (MacKay, eq. 10.2): R(pb) = C / (1 - H2(pb))
print()
for pb in (0.0, 0.001, 0.01, 0.05, 0.1):
    hb = h2(pb) if pb > 0 else 0.0
    print(f"pb = {pb:5.3f}:  R(pb) = {C / (1 - hb):.4f}")

# Output:
#    n    (1/n)log2 C(n,pn)   (1/n)log2 V(n,pn)   H2(0.1)
#     10   0.332193          0.345943           0.468996
#    100   0.439767          0.441423           0.468996
#   1000   0.464423          0.464592           0.468996
#  10000   0.468372          0.468389           0.468996
#
#    n    (1/n) log2( 2^n / V(n,pn) )   capacity 1 - H2(0.1)
#     10   0.654057                      0.531004
#    100   0.558577                      0.531004
#   1000   0.535408                      0.531004
#  10000   0.531611                      0.531004
#
# Pe >= -C/R + 1 - 1/(R*n)
# R = 0.8, C = 0.5310, n =   100:  Pe >= 0.3237
# R = 0.8, C = 0.5310, n =  1000:  Pe >= 0.3350
# R = 0.8, C = 0.5310, n = 10000:  Pe >= 0.3361
#
# pb = 0.000:  R(pb) = 0.5310
# pb = 0.001:  R(pb) = 0.5371
# pb = 0.010:  R(pb) = 0.5777
# pb = 0.050:  R(pb) = 0.7441
# pb = 0.100:  R(pb) = 1.0000

The first table shows that (1/n) log2 binom(n, pn) approaches H2(0.1) = 0.469 from below, 0.4684 at n = 10000, and that summing all the ways up to pn errors changes very little, because for p < 1/2 the last term dominates the sum. The second table compares the Hamming bound, which you meet in article 5, with the capacity: with t = pn the maximum rate of the Hamming sphere tends, from above, to 0.531, the same number. This numerical coincidence does not make the two statements equivalent: the Hamming bound concerns codes that correct every pattern of t errors, capacity concerns the typical noise. It remains a clue consistent with the count.

The converse: why above capacity the error does not fall

The negative side of the theorem has a short proof, which I derive here (it is the standard proof, not a quotation). Choose a message W uniform among 2nR; Fano's inequality gives H(W|Ŵ) ≤ 1 + Pe nR, and the fact that the channel is memoryless gives I(W;Ŵ) ≤ nC. Then nR = H(W) = H(W|Ŵ) + I(W;Ŵ) ≤ 1 + Pe nR + nC, whence Pe ≥ 1 - C/R - 1/(nR). The third block of the output solves it with SymPy and evaluates it for R = 0.8: the error stays above about 0.336 as n grows, never near zero. This is the correct form of the weak converse; the sentence "P_e ≥ (R - C) / log₂|A|" found in some lecture notes is not Fano's inequality.

For many channels the strong converse also holds: above capacity the error tends to 1, not just to a positive constant. Polyanskiy and Wu state it in theorem 22.1 (§22.1) for stationary memoryless channels with finite input or output alphabet, hence for the BSC and the BEC, and they attribute the weak converse to Fano's inequality (theorem 17.3). Section 22.2 instead constructs a channel with a countable input where the strong converse does not hold. I read the statements and the commentary sentences, not the proofs, which I do not reproduce.

The reliability exponent

MacKay (§10.7) gives a version with explicit dependence on N: there exist codes with pB ≤ exp(-N Er(R)), where Er(R) is the random-coding exponent of the channel, a convex, decreasing and positive function for 0 ≤ R < C. It vanishes as R → C: the closer you get to capacity, the longer the block must be to obtain the same error. MacKay also points out that, even for the BSC, there is no simple expression for Er(R); the base of the logarithm is a convention, which the text fixes at e. An analogous lower bound, with the sphere-packing exponent Esp(R), shows that the error cannot decrease faster than a certain exponential.

Common mistakes

• Writing C = 1 - p. For p = 0.01 it would give 0.99, while the capacity is 0.919; the difference is the entropy of the error.
• Reading "below capacity" as "without errors". The theorem guarantees an arbitrarily small error with long blocks, not a null one.
• Forgetting the maximization over the input. On the Z channel the uniform input gives 0.6788 instead of 0.6854.
• Adding capacities without hypotheses. For k memoryless channels used in parallel the capacity is the sum (proposition 19.10 of Polyanskiy and Wu); if the channels are not independent or there is memory, it must be rederived.
• Confusing the theorem with an algorithm. The statement is one of existence: for finite N it does not say which code to use nor what value N should take.

The common thread is that every mistake comes from turning an asymptotic statement, with a maximum to compute, into a rule to apply to a finite block. Remembering which hypotheses you used (memoryless channel, optimized input, block length tending to infinity) avoids all five.

What you can do now

You can compute the capacity of a discrete memoryless channel with three methods that check one another, and you know that for the BSC it is 1 - H2(p) while for the BEC it is 1 - f and for the Z channel a non-uniform input must be sought. You can state the second theorem with its three parts and say what it does not promise. You know that repetition pays for reliability with rate (35 repetitions to go below 10-9) and that there is a distance between this and the limit.

The next step is to see what constructed codes really do, starting from the Hamming distance, in article 5; the previous chapter, on lossless compression, is in article 3, and for the prerequisites the mutual information is in article 2 and the entropy in article 1; the series index gathers them all. The lossy compression limit that this theorem anticipates is in article 6.

Frequently asked questions