How can averages of a skewed source acquire a nearly Gaussian shape?

Predict which control changes the finite-sample mean distribution and which mainly changes how noisily its histogram estimates it. Inspect the first sample blocks before trusting the bell-shaped picture.

Intuition: study many means, not one running path

The previous chapter followed a single growing average. Here we repeatedly draw blocks of the same size $n$ and calculate one mean per block. Their collection has a distribution of its own. The source can remain skewed while the standardized mean distribution becomes more nearly Gaussian as $n$ grows.

A raw mean’s spread shrinks with $n$. To compare shape across block sizes, remove its theoretical center and divide by its theoretical standard deviation. This standardization keeps the reference spread fixed; it does not force observed bars to match a normal curve.

Definition

Let $X_1,X_2,\ldots$ be iid with finite mean $\mu$ and finite positive variance $\sigma^2$. For the mean of a block of size $n$, define

The classical Central Limit Theorem states $Z_n\xrightarrow{d}\mathcal N(0,1)$. Convergence in distribution means that, at every real $z$ here, $\mathbb P(Z_n\leq z)$ tends to the standard-normal CDF $\Phi(z)$. It does not make finite-$n$ observations exactly Gaussian or specify a universal accuracy threshold such as $n=30$.

If $\sigma^2=0$, the source is constant almost surely and this standardization is undefined. A normal overlay must not be presented as a valid standardized limit in that state.

Mechanics: one stream, contiguous sample blocks

Choose Bernoulli, Uniform$[0,1)$, or rate-$1$ exponential observations. All have known finite moments; the exponential source supplies a clearly nonnormal skewed population. Block size ranges from $1$ to $100$, and repetitions from $10$ to $1000$, bounding an iid run at $100000$ draws.

The canonical run uses the rate-$1$ exponential source, seed $42$, block size $n=8$, and $240$ repetitions. A uint32 seed drives the same explicit LCG used in the LLN chapter. Each iid block consumes its own contiguous portion of one stream. Open the sample details to inspect the first three complete actual samples retained along with their sums and means; all repeated means are standardized using the population moments. Reset restores canonical parameters, and replay reproduces the full run.

The dependence preset uses a uniform source and repeats one draw throughout each block. Its displayed iid standardization is explicitly a violated reference. In that example the block mean is just the single draw and the correct mean variance is $\sigma^2$, not $\sigma^2/n$.

Derivation: why the center and scale take this form

Under iid assumptions, linearity and zero cross-covariances give

Subtracting $\mu$ centers the mean, and division by $\sigma/\sqrt n$ gives $\mathbb E[Z_n]=0$ and $\operatorname{Var}(Z_n)=1$. These moment identities explain the scaling but do not prove the Gaussian limit. Many nonnormal distributions also have mean zero and variance one; convergence in distribution is the additional theorem.

An equivalent sum formulation subtracts $n\mu$ and divides by $\sigma\sqrt n$. Confusing sum and mean scales changes the variance of the standardized result and can create a misleadingly narrow or wide histogram.

Worked examples

Four rate-$1$ exponential observations with sum $6$ have mean $3/2$. Their population mean and standard deviation are both $1$, so

This is one standardized mean, not a proof about its distribution. For Bernoulli probability $1/4$ and $n=100$, the repeated means have theoretical expectation $1/4$ and variance $3/1600$.

For a fixed-window histogram with $1200$ repetitions, a width-$1/3$ bin containing $60$ observations has density $60/(1200\times1/3)=3/20$ and area $1/20$. If $1140$ observations lie inside the window, the displayed histogram area is $19/20$; the remaining empirical mass $1/20$ belongs to separately reported tails.

Executable reference code

The standard-library program generates contiguous sample blocks, retains three actual samples, computes sums and means, and handles zero variance explicitly. It verifies literal seeded blocks, replay, dependence, count conservation, and the full bounded workload without requiring favorable random agreement with the normal reference.

The NumPy program exposes a two-dimensional block array and averages along its sample axis. Its histogram first filters to the declared half-open window: NumPy otherwise includes the final right edge in the last bin. It divides counts by all repetitions and bin width, avoiding density=True renormalization of retained observations.

"""Bounded repeated means, honest fixed-window density, and an iid counterexample."""
from math import erfc
from math import log1p, sqrt


def uniforms(seed):
    if type(seed) is not int or not 0 <= seed <= 0xffffffff:
        raise ValueError("Seed must be a uint32 integer")
    state=seed
    while True:
        state=(1664525*state+1013904223)&0xffffffff
        yield state/4294967296


def population(source, p=50):
    if source not in ('bernoulli','uniform','exponential'):
        raise ValueError("Unknown population")
    if isinstance(p,bool) or not isinstance(p,(int,float)) or not 0 <= p <= 100 or abs(p*100-round(p*100))>1e-8:
        raise ValueError("Probability percent must be in [0,100] with at most two decimals")
    q=p/100
    mean,variance={'bernoulli':(q,q*(1-q)),'uniform':(.5,1/12),'exponential':(1.,1.)}[source]
    draw=lambda u: float(u<q) if source=='bernoulli' else u if source=='uniform' else -log1p(-u)
    return mean,variance,draw


def experiment(source='exponential',p=50,seed=42,n=8,reps=240,dependent=False):
    if type(n) is not int or not 1<=n<=100 or type(reps) is not int or not 10<=reps<=1000:
        raise ValueError("Require n in [1,100] and repetitions in [10,1000]")
    if type(dependent) is not bool:
        raise ValueError("Dependent flag must be boolean")
    mu,variance,draw=population(source,p)
    rng=uniforms(seed)
    means,sums,first=[],[],[]
    for rep in range(reps):
        if dependent:
            sample=[draw(next(rng))]*n
        else:
            sample=[draw(next(rng)) for _ in range(n)]
        if rep<3:
            first.append(sample)
        total=sum(sample)
        sums.append(total)
        means.append(total/n)
    if variance==0:
        return means,sums,first,None,None,None,None,variance
    z=[sqrt(n)*(x-mu)/sqrt(variance) for x in means]
    counts=[0]*24
    under=over=0
    for value in z:
        if value < -4:
            under+=1
        elif value >= 4:
            over+=1
        else:
            counts[min(23,int((value+4)*3))]+=1
    width=1/3
    density=[count/(reps*width) for count in counts]
    assert sum(counts)+under+over==reps
    assert abs(sum(density)*width-sum(counts)/reps)<1e-14
    actual_mean_variance=variance if dependent else variance/n
    return means,sums,first,z,counts,density,(under,over),actual_mean_variance


def main():
    r=experiment(n=1,reps=10,source='bernoulli')
    assert r[2]==[[1.],[1.],[0.]] and r[3][:3]==[1.,1.,-1.]
    assert r==experiment(n=1,reps=10,source='bernoulli')
    paired=experiment(n=2,reps=10,source='bernoulli')
    assert paired[2][0]==[1.,1.] and paired[2][1]==[0.,1.]
    assert experiment(p=0,source='bernoulli')[3] is None
    assert experiment(p=100,source='bernoulli')[4] is None
    dependent=experiment(n=10,reps=10,source='bernoulli',dependent=True)
    assert all(len(set(sample))==1 for sample in dependent[2])
    assert dependent[-1]==.25  # Not .25/10: the entries are perfectly dependent.
    assert 2*(1-.5*erfc(-4/sqrt(2)))>0  # Normal tails are small but nonzero.
    for source in ('bernoulli','uniform','exponential'):
        for n in (1,100):
            r=experiment(source,n=n,reps=1000)
            assert len(r[0])==1000 and len(r[2][0])==n
            assert sum(r[4])+sum(r[6])==1000
    for kwargs in ({'n':0},{'n':101},{'reps':9},{'reps':1001},{'seed':-1}):
        try:
            experiment(**kwargs)
        except ValueError:
            pass
        else:
            raise AssertionError("Invalid experiment accepted")
    print("First actual samples:",paired[2],"means:",paired[0][:3])
    print("Bin density = count/(all repetitions × width); omitted tails stay visible.")
    print("Increase n to study CLT shape; repetitions reduce histogram noise.")
    print("Zero variance: no standardization. Dependent blocks: no iid CLT claim.")


if __name__=='__main__':
    main()

"""NumPy blocks and histogram; the upper endpoint 4 belongs to overflow."""
import numpy as np
from math import log1p, sqrt


def uniforms(seed):
    if type(seed) is not int or not 0 <= seed <= 0xffffffff:
        raise ValueError("Seed must be a uint32 integer")
    state=seed
    while True:
        state=(1664525*state+1013904223)&0xffffffff
        yield state/4294967296


def population(source, p=50):
    if source not in ('bernoulli','uniform','exponential'):
        raise ValueError("Unknown population")
    if isinstance(p,bool) or not isinstance(p,(int,float)) or not 0 <= p <= 100 or abs(p*100-round(p*100))>1e-8:
        raise ValueError("Probability percent must be in [0,100] with at most two decimals")
    q=p/100
    mean,variance={'bernoulli':(q,q*(1-q)),'uniform':(.5,1/12),'exponential':(1.,1.)}[source]
    draw=lambda u: float(u<q) if source=='bernoulli' else u if source=='uniform' else -log1p(-u)
    return mean,variance,draw


def experiment(source='exponential',p=50,seed=42,n=8,reps=240,dependent=False):
    if type(n) is not int or not 1<=n<=100 or type(reps) is not int or not 10<=reps<=1000 or type(dependent) is not bool:
        raise ValueError("Require bounded integer n/reps and a boolean dependence flag")
    mu,variance,_=population(source,p)
    draws=reps if dependent else n*reps
    u=np.fromiter(uniforms(seed),dtype=float,count=draws)
    x=(u<p/100).astype(float) if source=='bernoulli' else u if source=='uniform' else -np.log1p(-u)
    blocks=np.repeat(x[:,None],n,axis=1) if dependent else x.reshape(reps,n)
    means=blocks.mean(axis=1)
    sums=blocks.sum(axis=1)
    if variance==0:
        return means,sums,blocks[:3],None,None,None,None,variance
    z=np.sqrt(n)*(means-mu)/np.sqrt(variance)
    inside=(z>=-4)&(z<4)
    edges=np.linspace(-4,4,25)
    counts,_=np.histogram(z[inside],bins=edges)
    density=counts/(reps*(1/3))  # density=True would incorrectly renormalize retained data.
    tails=(int(np.sum(z < -4)),int(np.sum(z >= 4)))
    assert int(counts.sum())+sum(tails)==reps
    np.testing.assert_allclose(np.sum(density)/3,np.mean(inside),rtol=0,atol=1e-14)
    return means,sums,blocks[:3],z,counts,density,tails,variance if dependent else variance/n


def main():
    r=experiment(source='bernoulli',n=1,reps=10)
    np.testing.assert_array_equal(r[2],[[1],[1],[0]])
    np.testing.assert_array_equal(r[3][:3],[1,1,-1])
    np.testing.assert_array_equal(r[0],experiment(source='bernoulli',n=1,reps=10)[0])
    pair=experiment(source='bernoulli',n=2,reps=10)
    np.testing.assert_array_equal(pair[2][:2],[[1,1],[0,1]])
    assert experiment(source='bernoulli',p=0)[3] is None
    dependent=experiment(source='bernoulli',n=10,reps=10,dependent=True)
    assert np.all(dependent[2]==dependent[2][:,:1]) and dependent[-1]==.25
    for source in ('bernoulli','uniform','exponential'):
        for n in (1,100):
            r=experiment(source,n=n,reps=1000)
            assert len(r[0])==1000 and int(r[4].sum())+sum(r[6])==1000
            np.testing.assert_allclose(r[1]/n,r[0],rtol=0,atol=1e-14)
    # np.histogram includes the final right edge; filter it out explicitly first.
    boundary=np.array([-5.,-4.,0.,4.,5.])
    kept=boundary[(boundary>=-4)&(boundary<4)]
    counts,_=np.histogram(kept,bins=np.linspace(-4,4,25))
    assert counts.sum()==2 and np.sum(boundary>=4)==2 and np.sum(boundary < -4)==1
    for kwargs in ({'n':0},{'reps':9},{'seed':2**32},{'p':float('inf')}):
        try:
            experiment(**kwargs)
        except ValueError:
            pass
        else:
            raise AssertionError("Invalid experiment accepted")
    print("Actual repeated means ≈",pair[0][:6],"first samples:",pair[2])
    print("Finite-n Bernoulli means remain on a lattice, however many repetitions are drawn.")
    print("Classical iid CLT requires finite positive variance; Cauchy does not qualify.")


if __name__=='__main__':
    main()

Explore

Start with a skewed exponential source and $n=1$. Predict its standardized shape, then increase $n$ while checking the actual first sample blocks. Look for reduced skew without assuming the noisy histogram must improve monotonically at every chosen size.

Compare the short preset at $n=1$ with the large preset at $n=30$, then hold $n$ fixed and increase repetitions. The theoretical mean distribution stays the same; its empirical histogram becomes a better estimate, subject to sampling variation. Switch to Bernoulli: finite-$n$ means remain on a lattice. More repetitions reveal that lattice more clearly rather than turning it into a continuous distribution.

Inspect the dependence example and compare its actual mean variance with the iid reference. At Bernoulli probability $0$ or $1$, identify the zero denominator and explain why there are deterministic means but no defined standardized histogram.

Common mistakes and counterexamples

The LLN concerns a mean approaching its target; the CLT concerns the distribution of its centered, rescaled fluctuations. Neither says the original observations become Gaussian. The population plot remains a separate object throughout the experiment.

Perfect within-block dependence breaks the iid variance calculation. If every entry equals $Y$, then $\bar X_n=Y$ and the iid reference $\sqrt n(Y-\mu)/\sigma$ has variance $n$, not $1$. An iid Cauchy population violates the classical theorem’s finite-variance and finite-mean assumptions; averaging it does not justify this Gaussian standardization.

There is no universal finite sample size at which every source becomes “normal enough.” Tail behavior, skew, event of interest, and desired accuracy matter. One seeded histogram supplies an illustration rather than evidence that all theorem assumptions hold or that an asymptotic approximation is accurate for a particular application.

Numerical considerations

The stage uses $24$ equal-width bins on $[-4,4)$, each of width $1/3$. Values below $-4$ and values at or above $4$ are counted separately. Bin density is

where $R$ is all repetitions and $\Delta$ is bin width. The histogram area is the retained empirical fraction, not automatically one. The analytic standard-normal reference is normalized over the entire real line and has positive tail mass outside every finite window; its omitted mass beyond $\pm4$ is approximately $0.000063342484$.

The normal overlay uses the Gaussian density and CDF policy from the preceding distribution chapter; CDF approximations have a stated absolute error budget of $10^{-9}$, so a difference has budget $2\times10^{-9}$. This numerical error is separate from Monte Carlo noise, finite-$n$ approximation error, binning, finite drawing bounds, and pseudo-random stream limitations. None is removed by aligning plotted heights manually.

Engineering and ML connections

The CLT motivates approximate Gaussian descriptions of accumulated independent finite-variance noise and some Monte Carlo estimation errors. The raw average’s standard-deviation scale is $\sigma/\sqrt n$; using it requires checking the dependence and moment assumptions, and the approximation’s usefulness depends on the event being estimated.

Batch averages in learning systems often contain correlated observations, and heavy-tailed quantities can violate classical assumptions. The normal reference is a model to assess, not a consequence of drawing a smooth curve over data. Statistical confidence procedures build additional decisions on this probability foundation in the next part of the book.

Exercises

Calculate standardization and histogram area exactly where possible. Explain what each sampling control changes and which assumption a counterexample violates.

Summary and checkpoints

The classical iid CLT standardizes repeated means using finite positive population variance and describes their limiting distribution. Check that you can recover a standardized value from an actual sum, separate n from repetitions, preserve histogram tail mass, handle zero variance and dependence, and keep one finite simulation distinct from a theorem.