<jupyter_start><jupyter_text>How to Visualize the Central Limit Theorem Using Python> as described in <NAME>'s "Think Stats"- toc: false- badges: true- comments: true- categories: [statistics, central limit theorem, python]- author: <NAME>- image: images/central_limit_theorem.png- search_exclude: false- permalink: /central-limit-theorem > Note:  The code required to replicate the visualizations below is provided in the ```Show Code``` dropdowns, but the modules ```thinkstats2``` and ```thinkplot``` must first be pulled from [here](https://github.com/AllenDowney/ThinkStats2).<jupyter_code>#collapse-hide
%matplotlib inline

import numpy as np
import pandas as pd

from scipy.stats import norm
from scipy.stats import expon

import random

import thinkstats2
import thinkplot

from matplotlib import MatplotlibDeprecationWarning

from warnings import simplefilter
simplefilter('ignore', MatplotlibDeprecationWarning)<jupyter_output><empty_output><jupyter_text>tl;drIf you add up independent variates from a distribution with finite mean and variance, the sum converges on a normal distribution. The sum of exponential variates converges to normal as sample size increases. The lognormal distribution has higher variance, so it requires a larger sample size before it converges to normal. The Pareto distribution violates the requirements of the CLT and does not generally converge to normal. If the random variates are correlated, that also violates the CLT, so the sums don't generally converge.  Distribution of Independent Variates with Finite Mean and VarianceIf you add up independent variates from a distribution with finite mean and variance, the sum converges on a normal distribution. The following plot shows how the sum of exponential variates converges to normal as sample size increases.<jupyter_code>#collapse-hide
# The following function generates samples with difference
# sizes from an exponential distribution.

def MakeExpoSamples(beta=2.0, iters=1000):
    """Generates samples from an exponential distribution.

    beta: parameter
    iters: number of samples to generate for each size

    returns: list of samples
    """
    samples = []
    for n in [1, 10, 100]:
        sample = [np.sum(np.random.exponential(beta, n))
                  for _ in range(iters)]
        samples.append((n, sample))
    return samples

# This function generates normal probability plots for
# samples with various sizes.

def NormalPlotSamples(samples, plot=1, ylabel=''):
    """Makes normal probability plots for samples.

    samples: list of samples
    label: string
    """
    for n, sample in samples:
        thinkplot.SubPlot(plot)
        thinkstats2.NormalProbabilityPlot(sample)

        thinkplot.Config(title='n=%d' % n,
                         legend=False,
                         xticks=[],
                         yticks=[],
                         xlabel='Random normal variate',
                         ylabel=ylabel)
        plot += 1

# Make plots.        
        
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeExpoSamples()
NormalPlotSamples(samples, plot=1,
                  ylabel='Sum of expo values')<jupyter_output><empty_output><jupyter_text><jupyter_code>#hide_input
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeExpoSamples()
NormalPlotSamples(samples, plot=1,
                  ylabel='Sum of expo values')<jupyter_output><empty_output><jupyter_text>Lognormal DistributionThe lognormal distribution has higher variance, so it requires a larger sample size before it converges to normal.<jupyter_code>#collapse-hide
# The following function generates samples with difference
# sizes from an lognormal distribution.
def MakeLognormalSamples(mu=1.0, sigma=1.0, iters=1000):
    """Generates samples from a lognormal distribution.

    mu: parmeter
    sigma: parameter
    iters: number of samples to generate for each size

    returns: list of samples
    """
    samples = []
    for n in [1, 10, 100]:
        sample = [np.sum(np.random.lognormal(mu, sigma, n))
                  for _ in range(iters)]
        samples.append((n, sample))
    return samples

# Make plots.

thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeLognormalSamples()
NormalPlotSamples(samples, ylabel='sum of lognormal values')<jupyter_output><empty_output><jupyter_text><jupyter_code>#hide_input
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeLognormalSamples()
NormalPlotSamples(samples, ylabel='sum of lognormal values')<jupyter_output><empty_output><jupyter_text>Pareto DistributionThe Pareto distribution has infinite variance, and sometimes infinite mean, depending on the parameters.  It violates the requirements of the CLT and does not generally converge to normal.<jupyter_code>#collapse-hide
# The following function generates samples with difference
# sizes from an Pareto distribution.
def MakeParetoSamples(alpha=1.0, iters=1000):
    """Generates samples from a Pareto distribution.

    alpha: parameter
    iters: number of samples to generate for each size

    returns: list of samples
    """
    samples = []

    for n in [1, 10, 100]:
        sample = [np.sum(np.random.pareto(alpha, n))
                  for _ in range(iters)]
        samples.append((n, sample))
    return samples

# Make plots.

thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeParetoSamples()
NormalPlotSamples(samples, ylabel='sum of Pareto values')<jupyter_output><empty_output><jupyter_text><jupyter_code>#hide_input
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeParetoSamples()
NormalPlotSamples(samples, ylabel='sum of Pareto values')<jupyter_output><empty_output><jupyter_text>Distribution of Correlated Random VariatesIf the random variates are correlated, that also violates the CLT, so the sums don't generally converge. To generate correlated values, we generate correlated normal values and then transform to whatever distribution we want.<jupyter_code>#collapse-hide
def GenerateCorrelated(rho, n):
    """Generates a sequence of correlated values from a standard normal dist.
    
    rho: coefficient of correlation
    n: length of sequence

    returns: iterator
    """
    x = random.gauss(0, 1)
    yield x

    sigma = np.sqrt(1 - rho**2)
    for _ in range(n-1):
        x = random.gauss(x * rho, sigma)
        yield x

def GenerateExpoCorrelated(rho, n):
    """Generates a sequence of correlated values from an exponential dist.

    rho: coefficient of correlation
    n: length of sequence

    returns: NumPy array
    """
    normal = list(GenerateCorrelated(rho, n))
    uniform = norm.cdf(normal)
    expo = expon.ppf(uniform)
    return expo

def MakeCorrelatedSamples(rho=0.9, iters=1000):
    """Generates samples from a correlated exponential distribution.

    rho: correlation
    iters: number of samples to generate for each size

    returns: list of samples
    """    
    samples = []
    for n in [1, 10, 100]:
        sample = [np.sum(GenerateExpoCorrelated(rho, n))
                  for _ in range(iters)]
        samples.append((n, sample))
    return samples

# Make plots.

thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeCorrelatedSamples()
NormalPlotSamples(samples, ylabel='Sum of correlated exponential values')<jupyter_output><empty_output><jupyter_text><jupyter_code>#hide_input
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeCorrelatedSamples()
NormalPlotSamples(samples, ylabel='Sum of correlated exponential values')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08583633693765852, 'fraction_numerical': 0.011613936724068883, 'mean_word_length': 3.9913391072618256, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 41, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '165926', 'license': 'Apache-2.0', 'n_tokens_mistral': 2253, 'n_tokens_neox': 2095, 'n_words': 725, 'path': '_notebooks/2020-04-30-central-limit-theorem.ipynb', 'repo_name': 'nickmccarty/til', 'size': '7510'}