Capstone D: Random Structures Laboratory

Random models with controlled parameters

The laboratory distinguishes two questions: when does a random graph contain a component occupying a positive fraction of its vertices, and when are all its vertices connected? A giant component can coexist with many isolated vertices. The corresponding thresholds occur at different edge-probability scales, so one plot labeled “phase transition” is not an adequate specification.

All graphs are finite, simple, undirected, and labeled on $\{0,\ldots,n-1\}$, with $n\ge1$. Write $N=\binom n2$. In $G(n,p)$, each of the $N$ possible edges is present independently with probability $p$. In $G(n,m)$, choose a uniformly random $m$-element subset of the possible edges. The first model has random edge count $M\sim\operatorname{Bin}(N,p)$; the second has exactly $m$ edges and dependent edge indicators. Matching their expected counts does not make the distributions equal.

The prerequisites the probabilistic method, random graphs and random combinatorial structures, information-theoretic and entropy methods, and randomized algorithms and approximate counting support the project. The finite calculations and statistical assumptions used here are stated explicitly.

Independent generation and model validation

Under ideal independent uniform draws $U_e$ on $[0,1)$, the rule “include edge $e$ when $U_e\lt p$” produces independent Bernoulli indicators. Uniform sampling without replacement produces $G(n,m)$ because each $m$-subset is equally likely. The implementation uses random.Random and sample; the statistical interpretation assumes that these pseudorandom streams adequately approximate the ideal draws. A seed reproduces a computation, but it is neither a mathematical source of randomness nor evidence that the sampling rule is correct.

Validate the model at several levels. Boundary tests at $p=0,1$ and $m=0,N$ check empty and complete graphs. Exact enumeration on four vertices checks the distributions and graph measurements without relying on any seed. A separate 10,000-sample audit records edge-count histograms, two edge marginals, their joint frequency, and empirical covariance. Agreement is diagnostic evidence, not proof of independence; a few successful moment checks cannot identify a joint distribution.

Every main sample's component-size list is computed by both graph traversal and disjoint-set union. The algorithms must agree on the entire sorted list, not merely connectivity. The code also checks distinct edges, endpoint ordering, and fixed edge counts in $G(n,m)$. Its model audits deliberately report random discrepancies rather than asserting that sampled moments equal their theoretical values.

Save the panel as random_lab.py and run python3 random_lab.py in an empty output directory. It writes random-data.json and the standalone plot random-structures.svg. It needs only the Python standard library; the reference run used Python 3.12.14. The plot has separate panels and pointwise confidence bars, with lines connecting settings rather than fitting a threshold.

from collections import Counter
from itertools import combinations
from math import comb, exp, log, sqrt
from pathlib import Path
import hashlib
import json
import random
import sys

BASE, REPEATS = 20261013, 200


def gnp(n, p, rng):
    assert n >= 1 and 0 <= p <= 1
    return [(u, v) for u, v in combinations(range(n), 2) if rng.random() < p]


def gnm(n, m, rng):
    pairs = list(combinations(range(n), 2))
    assert n >= 1 and 0 <= m <= len(pairs)
    return rng.sample(pairs, m)


def measurements(n, edges):
    adj = [set() for _ in range(n)]
    parent = list(range(n))
    def root(v):
        while parent[v] != v:
            parent[v] = parent[parent[v]]
            v = parent[v]
        return v
    assert len(set(edges)) == len(edges)
    for u, v in edges:
        assert 0 <= u < v < n
        adj[u].add(v)
        adj[v].add(u)
        parent[root(u)] = root(v)
    seen, sizes = set(), []
    for start in range(n):
        if start in seen:
            continue
        stack, size = [start], 0
        seen.add(start)
        while stack:
            v = stack.pop()
            size += 1
            for w in adj[v] - seen:
                seen.add(w)
                stack.append(w)
        sizes.append(size)
    sizes.sort(reverse=True)
    independent = sorted(Counter(root(v) for v in range(n)).values(), reverse=True)
    assert sizes == independent and sum(sizes) == n
    return {'edges': len(edges), 'connected': int(len(sizes) == 1),
            'L1': sizes[0], 'L2': sizes[1] if len(sizes) > 1 else 0,
            'isolated': sum(not row for row in adj)}


def exact_checks():
    pairs = list(combinations(range(4), 2))
    graphs = [[e for i, e in enumerate(pairs) if mask >> i & 1]
              for mask in range(64)]
    assert [sum(len(e) == k for e in graphs) for k in range(7)] == [comb(6, k) for k in range(7)]
    fixed = [e for e in graphs if len(e) == 3]
    assert len(fixed) == 20
    assert sum(pairs[0] in e for e in graphs) == 32
    assert sum(pairs[0] in e and pairs[1] in e for e in graphs) == 16
    assert sum(pairs[0] in e for e in fixed) == 10
    assert sum(pairs[0] in e and pairs[1] in e for e in fixed) == 4
    assert sum(measurements(4, e)['connected'] for e in graphs) == 38
    assert sum(measurements(4, e)['connected'] for e in fixed) == 16
    assert gnp(4, 0, random.Random(BASE)) == []
    assert len(gnp(4, 1, random.Random(BASE))) == 6
    assert gnm(4, 0, random.Random(BASE)) == []
    assert len(gnm(4, 6, random.Random(BASE))) == 6
    return {'G4_half_connected': [38, 64], 'G4_3_connected': [16, 20],
            'edge_count_histogram': [1, 6, 15, 20, 15, 6, 1]}


def model_audit(model):
    rng = random.Random(BASE - (1 if model == 'Gnp' else 2))
    histogram, hits, joint = Counter(), [0, 0], 0
    trials, n, p, m = 10000, 6, 0.4, 6
    for _ in range(trials):
        edges = set(gnp(n, p, rng) if model == 'Gnp' else gnm(n, m, rng))
        histogram[len(edges)] += 1
        a, b = int((0, 1) in edges), int((0, 2) in edges)
        hits[0] += a
        hits[1] += b
        joint += a * b
    mean = sum(k * v for k, v in histogram.items()) / trials
    variance = sum((k - mean) ** 2 * v for k, v in histogram.items()) / trials
    return {'model': model, 'trials': trials, 'mean_edges': mean,
            'variance_edges': variance, 'edge_marginals': [h / trials for h in hits],
            'joint_probability': joint / trials,
            'covariance': joint / trials - hits[0] * hits[1] / trials ** 2,
            'histogram': sorted(histogram.items())}


def finite_bounds(n, p):
    mu = n * (1 - p) ** (n - 1)
    second = mu + n * (n - 1) * (1 - p) ** (2 * n - 3)
    cut = sum(comb(n, k) * (1 - p) ** (k * (n - k)) for k in range(1, n // 2 + 1))
    return {'expected_isolated': mu,
            'connectivity_lower': max(0.0, 1 - cut),
            'connectivity_upper': 1 - mu * mu / second if second else 1.0}


def rho(c):
    if c <= 1:
        return 0.0
    value = 1.0
    for _ in range(500):
        value = 1 - exp(-c * value)
    return value


def interval(mean, repeats, delta=0.05):
    radius = sqrt(log(2 / delta) / (2 * repeats))
    return [max(0.0, mean - radius), min(1.0, mean + radius)]


def plot(points):
    colors = {(60, 'Gnp'): '#d44927', (120, 'Gnp'): '#236ca3',
              (60, 'Gnm'): '#38804d', (120, 'Gnm'): '#855ca5'}
    out = ['<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 1200 510">',
           '<rect width="1200" height="510" fill="white"/>',
           '<g font-family="Arial" font-size="14" fill="#303841">']
    for panel, (regime, title, field) in enumerate([
            ('giant', 'Largest component: sample mean L1/n', 'mean_L1_fraction'),
            ('connectivity', 'Connectivity: sample success fraction', 'connected_fraction')]):
        left, top, width, height = 85 + panel * 585, 80, 430, 290
        lo, hi = (0.5, 2.0) if regime == 'giant' else (-2.0, 2.0)
        x = lambda v: left + width * (v - lo) / (hi - lo)
        y = lambda v: top + height * (1 - v)
        out.append(f'<text x="{left}" y="38" font-size="19">{title}</text>')
        for tick in (0, 0.25, 0.5, 0.75, 1):
            out.append(f'<path d="M {left} {y(tick)} h {width}" stroke="#dde2e5"/>')
            out.append(f'<text x="{left-12}" y="{y(tick)+5}" text-anchor="end">{tick}</text>')
        for parameter in sorted({p['parameter'] for p in points if p['regime'] == regime}):
            out.append(f'<text x="{x(parameter)}" y="395" text-anchor="middle">{parameter:g}</text>')
        for key, color in colors.items():
            group = sorted([p for p in points if (p['n'], p['model']) == key and p['regime'] == regime], key=lambda p: p['parameter'])
            poly = ' '.join(f"{x(p['parameter']):.3f},{y(p[field]):.3f}" for p in group)
            out.append(f'<polyline points="{poly}" fill="none" stroke="{color}" stroke-width="2"/>')
            for p in group:
                cx, cy = x(p['parameter']), y(p[field])
                ci = p['L1_mean_interval'] if regime == 'giant' else p['connectivity_interval']
                out.append(f'<path d="M {cx} {y(ci[0])} V {y(ci[1])}" stroke="{color}" opacity="0.4"/>')
                out.append(f'<circle cx="{cx}" cy="{cy}" r="4" fill="{color}"/>')
        out.append(f'<text x="{left+width/2}" y="425" text-anchor="middle">{"c in p=c/n" if regime == "giant" else "t in p=(log n+t)/n"}; Gnm rounds expected edges</text>')
    for i, (key, color) in enumerate(colors.items()):
        xpos = 90 + i * 280
        out.append(f'<text x="{xpos}" y="463" fill="{color}">n={key[0]}, {key[1]}, 200 trials</text>')
    out.append('<text x="600" y="492" text-anchor="middle">Bars: pointwise 95% Hoeffding intervals; lines connect sampled settings.</text></g></svg>')
    Path('random-structures.svg').write_text('\n'.join(out), encoding='utf-8')


def main():
    checks = exact_checks()
    audits = [model_audit(model) for model in ('Gnp', 'Gnm')]
    points = []
    for n in (60, 120):
        for regime, parameters in [('giant', (0.5, 1.0, 1.5, 2.0)),
                                   ('connectivity', (-2.0, -1.0, 0.0, 1.0, 2.0))]:
            for parameter in parameters:
                p = parameter / n if regime == 'giant' else (log(n) + parameter) / n
                m = round(comb(n, 2) * p)
                for model in ('Gnp', 'Gnm'):
                    seed = BASE + 1009 * (len(points) + 1)
                    rng = random.Random(seed)
                    samples = [measurements(n, gnp(n, p, rng) if model == 'Gnp' else gnm(n, m, rng)) for _ in range(REPEATS)]
                    connected = sum(s['connected'] for s in samples) / REPEATS
                    largest = sum(s['L1'] / n for s in samples) / REPEATS
                    point = {'n': n, 'regime': regime, 'parameter': parameter,
                             'model': model, 'p': p, 'm': m, 'seed': seed,
                             'repeats': REPEATS, 'samples': samples,
                             'connected_fraction': connected, 'mean_L1_fraction': largest,
                             'connectivity_interval': interval(connected, REPEATS),
                             'L1_mean_interval': interval(largest, REPEATS),
                             'finite_Gnp_bounds': finite_bounds(n, p) if model == 'Gnp' else None,
                             'limit_reference': (rho(parameter) if regime == 'giant' else exp(-exp(-parameter))) if model == 'Gnp' else None}
                    if model == 'Gnm':
                        assert all(s['edges'] == m for s in samples)
                    points.append(point)
    payload = {'base_seed': BASE, 'python': sys.version.split()[0],
               'exact_checks': checks, 'model_audits': audits, 'points': points}
    data = json.dumps(payload, sort_keys=True, separators=(',', ':')) + '\n'
    Path('random-data.json').write_text(data, encoding='utf-8')
    plot(points)
    print('Python:', sys.version.split()[0], 'settings:', len(points), 'graphs:', len(points) * REPEATS)
    print('SHA256:', hashlib.sha256(data.encode()).hexdigest())
    print('audits:', json.dumps(audits, sort_keys=True))
    for p in points:
        if p['n'] == 120 and p['model'] == 'Gnp' and ((p['regime'] == 'giant' and p['parameter'] == 2) or (p['regime'] == 'connectivity' and p['parameter'] in (0, 2))):
            print('selected:', json.dumps({k: p[k] for k in ('regime', 'parameter', 'connected_fraction', 'mean_L1_fraction', 'connectivity_interval', 'finite_Gnp_bounds')}, sort_keys=True))


if __name__ == '__main__':
    main()

The exact four-vertex checks have no Monte Carlo error. In $G(4,1/2)$, all 64 graphs are equiprobable, the edge-count histogram is $1,6,15,20,15,6,1$, and connectivity probability is $38/64=19/32$. In $G(4,3)$, there are 20 equiprobable graphs and connectivity probability is $16/20=4/5$. A specified edge has marginal probability $1/2$ in both models, but two specified distinct edges have joint probability $1/4$ and $1/5$, respectively.

For the audit at $n=6$, compare $G(6,0.4)$ with $G(6,6)$. The theoretical mean edge count is six in each. The first model's variance is $15(0.4)(0.6)=3.6$; the second model's variance is zero. Distinct edge indicators in $G(n,m)$ have joint probability $m(m-1)/(N(N-1))$ and covariance $-(m/N)(1-m/N)/(N-1)$. Their independence in $G(n,p)$ does not imply that vertex degrees or component indicators are independent.

The reference audit observed mean 5.9859 and empirical variance 3.64170119 in $G(6,0.4)$; in $G(6,6)$ they were exactly six and zero. The joint edge frequencies were 0.1623 and 0.1455, compared with theoretical values 0.16 and $1/7$. The empirical covariances were approximately 0.004314 and $-0.016391$. The variance here uses divisor 10,000 and the empirical mean; it is the empirical central moment, not an unbiased variance estimator.

Exercise: why does the sample generator choose from unordered pairs rather than independently choosing each directed adjacency entry? Solution: the model is undirected. Independent directed choices would require symmetrization, generally changing edge probabilities and possibly introducing unintended dependencies.

Connectivity and giant-component experiments

For fixed $0\le c\lt1$ in $G(n,c/n)$, the largest component has order $O(\log n)$ with probability tending to one. For fixed $c\gt1$, the normalized largest size converges in probability to the unique positive solution $\rho(c)$ of

$$\rho=1-e^{-c\rho},$$

and the second-largest component remains logarithmic with high probability. The critical setting $c=1$ is excluded from those two size-order assertions; its finite samples are recorded, but their behavior should not be described by the subcritical logarithmic estimate. The reference fraction is zero for $c\le1$. The global giant-component proof is an imported theorem; local branching intuition alone does not establish concentration or uniqueness. Primary background includes Erdős and Rényi's original study and Bollobás and Riordan's exploration treatment.

Connectivity instead occurs near $p=(\log n)/n$. For fixed $t\in\mathbb R$, the classical window theorem gives

$$\Pr\bigl(G(n,(\log n+t)/n)\text{ is connected}\bigr)\longrightarrow e^{-e^{-t}}.$$

Its proof requires both the limiting isolated-vertex law and control of other disconnected components; matching the expected number of isolated vertices is insufficient. This exact limiting law is imported, while finite bounds and the coarser threshold comparison are proved below. The references in the dataset apply these formulas only to $G(n,p)$; transferring them to the fixed-edge model is not proved by rounding $Np$.

Reference results at $n=120$: at $c=2$, the average $L_1/n$ was 0.789167, close to $\rho(2)\approx0.796812$, yet connectivity occurred in zero of 200 samples. At connectivity parameter $t=0$, 88 of 200 samples were connected, giving 0.44, while average $L_1/n$ was 0.992417. At $t=2$, 180 of 200 were connected, giving 0.90, while average $L_1/n$ was 0.999083. These are sample statistics at finite $n$, not exact model probabilities.

The giant panel plots mean $L_1/n$ against $c$; the connectivity panel plots the connected fraction against $t$. The axes are intentionally different. The JSON retains $L_2$ and isolation counts for further examination. A nearly spanning largest component can still leave one vertex isolated, so mean $L_1/n$ near one does not imply connectivity probability near one.

Exercise: at $p=2/n$, what does the expected isolated count suggest as $n$ grows? Solution: it is asymptotic to $ne^{-2}$, which grows linearly. The second-moment calculation below shows that isolates occur with probability tending to one in this regime, even though a giant component exists.

Probability bounds beside empirical estimates

Let $I$ count isolated vertices in $G(n,p)$, with $n\ge2$. A vertex is isolated when its $n-1$ incident edges are absent, so the first moment is

$$\mu:=\mathbb E I=n(1-p)^{n-1}.$$

Two distinct vertices are both isolated when all $2n-3$ edges incident with at least one of them are absent. Consequently

$$\mathbb E I^2=\mu+n(n-1)(1-p)^{2n-3}.$$

Markov's inequality gives $\Pr(I>0)\le\min\{1,\mu\}$. For the complementary direction, Cauchy–Schwarz applied to $I\mathbf1_{\{I>0\}}$ gives

$$\Pr(I>0)\ge\frac{\mu^2}{\mathbb E I^2}.$$

Since connectivity implies no isolated vertices, this is an upper bound $\Pr(\text{connected})\le1-\mu^2/\mathbb E I^2$. It is not a lower bound on connectivity. Handle $p=1$ separately when the moment ratio is $0/0$.

For a lower connectivity bound, every disconnected graph has a component with some vertex set $S$ of size $1\le k\le n/2$. No edge crosses from $S$ to its complement. Union-bounding over such subsets yields

$$\Pr(\text{disconnected})\le\min\left\{1,\sum_{k=1}^{\lfloor n/2\rfloor}\binom nk(1-p)^{k(n-k)}\right\}.$$

The subsets need not be internally connected in this bound, so overcounting can make it weak. The resulting connectivity lower bound is zero when the sum exceeds one. A correct but uninformative bound should remain labeled uninformative rather than replaced by an empirical fit.

Solved finite check: at $n=4,p=1/2$, $\mu=1/2$ and $\mathbb E I^2=7/8$, giving connectivity at most $5/7$. The cut sum is $4/8+6/16=7/8$, giving connectivity at least $1/8$. The independently enumerated probability $19/32$ lies between these bounds. Having no isolated vertex is not sufficient: the three disjoint-edge matchings are disconnected without isolates.

The moment ratio also proves a regime statement. At $p=c/n$ with fixed $c>0$, $\mu\sim ne^{-c}$ and $\mathbb E I^2/\mu^2\to1$. Thus $\Pr(I>0)\to1$ and connectivity probability tends to zero. At $p=(1-\varepsilon)\log n/n$, fixed $0<\varepsilon<1$, $\mu\sim n^\varepsilon$ and the same ratio tends to one, proving disconnection with high probability below the connectivity scale.

For completeness, the cut bound proves connectivity above that scale at $p=(1+\varepsilon)\log n/n$. Put $a=\varepsilon/(2(1+\varepsilon))$. For $k\le an$, each summand is at most $(e n^{-\varepsilon/2}/k)^k$ using $\binom nk\le(en/k)^k$; their sum tends to zero. For $an\lt k\le n/2$, each cut has at least $an^2/2$ potential edges, and summing binomial coefficients gives at most $2^n\exp(-pa n^2/2)$, also tending to zero. This proves the coarse threshold distinction without claiming the fixed-$t$ window law from the same argument.

These formulas use independent absent-edge probabilities in $G(n,p)$. They are not the formulas for $G(n,m)$, where avoidance of a specified cut of $b$ possible edges has probability $\binom{N-b}{m}/\binom Nm$, interpreted as zero if $m>N-b$.

Exercise: why does $\mathbb E I\to\infty$ alone not prove an isolate appears with high probability? Solution: a large mean can be driven by rare outcomes with very large counts. The second moment rules out that explanation in the regimes just analyzed.

Uncertainty and finite-size effects

There are three separate discrepancies: Monte Carlo fluctuation around a finite-model expectation, finite-size disagreement between that expectation and an asymptotic limit, and possible implementation or model errors. Increasing the trial count reduces the first, not the other two. Increasing graph size studies the second, but does not repair a wrong generator or prove a theorem from a plot.

The report uses a finite-sample concentration interval for both observables. If $Y_1,\ldots,Y_R$ are independent observations in $[0,1]$ with common mean $\theta$, Hoeffding's inequality gives

$$\Pr\bigl(\lvert\bar Y-\theta\rvert>r\bigr)\le2e^{-2Rr^2}.$$

Choose $r=\sqrt{\log(2/\delta)/(2R)}$ and clip $[\bar Y-r,\bar Y+r]$ to $[0,1]$. This has coverage at least $1-\delta$ under the independent-trial model. It applies to the Bernoulli connectivity indicator and the bounded fraction $L_1/n$ in either graph model, despite within-graph edge dependence in $G(n,m)$. It concerns a population mean, not an interval containing 95% of individual component sizes.

With $R=200,\delta=0.05$, the radius is approximately 0.096032. The reference connectivity intervals for the three selected settings are approximately $[0,0.096032]$, $[0.343968,0.536032]$, and $[0.803968,0.996032]$. The largest-component mean at $c=2$ has interval approximately $[0.693134,0.885199]$. These conservative intervals remain meaningful at zero observed successes. Under an ideal Bernoulli model, inversion of the zero-success probability gives the sharper one-sided upper value $1-0.05^{1/200}\approx0.014867$; zero successes still does not establish probability zero.

At $n=120,t=2$, the finite theoretical connectivity bounds evaluate to approximately $[0.874416,0.895415]$, while the observed fraction is 0.90. This is no contradiction: the bound constrains the true probability, and a sample fraction may lie above it. The concentration interval contains the theoretical range. Formula-derived bounds are mathematically exact statements, but their displayed floating-point evaluations and rounded endpoints are not certified numerical enclosures.

The plot's confidence bars are pointwise, not a 95% guarantee that every bar simultaneously covers. The dataset contains 72 mean estimates across 36 settings; the two-panel plot displays 36 of them. A union bound can provide simultaneous coverage by using $\delta=0.05/72$ per estimate; that produces wider intervals. No independence between the two observables on one graph is needed for this union bound. Repetition counts and parameter grids were fixed before sampling; do not apply fixed-sample intervals unchanged after an outcome-dependent stopping rule.

The limiting connectivity reference at $t=0$ is $e^{-1}\approx0.367879$. A finite estimate near it supports an illustrative comparison, not a proof of convergence; a discrepancy may combine finite-size bias and sampling error. With only two graph sizes, do not estimate an asymptotic rate or label the finite curve a discontinuous transition. The original data and parameterization are more informative than a smoothed picture without uncertainty bars.

Exercise: what happens to the Hoeffding radius when the trial count is quadrupled? Solution: it is halved. The finite-size bias of the model expectation relative to the limiting theorem is unaffected.

Random structures report and assessment

A complete submission includes random_lab.py, the exact seeded observation file random-data.json, the two-panel plot random-structures.svg, environment and invocation, exact four-vertex checks, generation audits, theoretical derivations, and an uncertainty discussion. The JSON stores every structural sample's five observables, each setting's model parameters and seed, finite $G(n,p)$ bounds, limiting references, and the empirical intervals. It does not store every sampled edge set; replaying the declared stream regenerates those sets.

The reference run used Python 3.12.14 and produced 7,200 structural samples plus 20,000 generation-audit samples. Its data SHA-256 is edec466337ec26cae27e088d85474bec21bd235882f66a7570d1ce77375e65e6. Sorted JSON keys, compact separators, UTF-8, and a final newline define these bytes. The interpreter version is itself stored, and random sampling algorithms or floating-point outputs can change across environments, so record any differing environment rather than treating a seed as a universal cross-version contract.

Interpret the exact small-instance distributions as model checks, the audits as statistical diagnostics, the structural samples as estimates, the moment and cut arguments as proofs, and the imported asymptotic laws as theorems with their stated regimes. The global giant proof and the connectivity-window proof are not reproduced here. Numerical solution of the giant equation supplies a reference curve, not those missing proofs. Any newly suggested regularity outside the declared range must remain exploratory until separately justified.

A cross-disciplinary application is network reliability. Connectivity measures whether every site can communicate, while a large component measures service to most sites. The contrast matters for systems where a small isolated population is unacceptable. Independent edge failures may be an unsuitable physical model when failures share causes; applying this laboratory to infrastructure requires checking those assumptions rather than borrowing its threshold curves unchanged.