pareto_world.py

"""This file contains code used in "Think Stats",
by Allen B. Downey, available from greenteapress.com

Copyright 2011 Allen B. Downey
License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html
"""

from __future__ import print_function

import numpy as np
import random
import scipy.stats

import thinkstats2
import thinkplot

from collections import Counter


def ParetoCdf(x, alpha, xmin):
    """Evaluates CDF of the Pareto distribution with parameters alpha, xmin."""
    if x < xmin:
        return 0
    return 1 - pow(x / xmin, -alpha)


def ParetoMedian(xmin, alpha):
    """Computes the median of a Pareto distribution."""
    return xmin * pow(2, 1/alpha)


def MakeParetoCdf():
    """Generates a plot of the CDF of height in Pareto World."""
    n = 50
    max = 1000.0
    xs = [max*i/n for i in range(n)]
    
    xmin = 100
    alpha = 1.7
    ps = [ParetoCdf(x, alpha, xmin) for x in xs]
    print('Median', ParetoMedian(xmin, alpha))
    
    pyplot.clf()
    pyplot.plot(xs, ps, linewidth=2)
    myplot.Save('pareto_world1',
                title = 'Pareto CDF',
                xlabel = 'height (cm)',
                ylabel = 'CDF',
                legend=False)


def MakeFigure(xmin=100, alpha=1.7, mu=150, sigma=25):
    """Makes a figure showing the CDF of height in ParetoWorld.

    Compared to a normal distribution.

    xmin: parameter of the Pareto distribution
    alpha: parameter of the Pareto distribution
    mu: parameter of the Normal distribution
    sigma: parameter of the Normal distribution
    """

    t1 = [xmin * random.paretovariate(alpha) for i in range(10000)]
    cdf1 = Cdf.MakeCdfFromList(t1, name='pareto')

    t2 = [random.normalvariate(mu, sigma) for i in range(10000)]
    cdf2 = Cdf.MakeCdfFromList(t2, name='normal')

    myplot.Clf()
    myplot.Cdfs([cdf1, cdf2])
    myplot.Save(root='pareto_world2',
                title='Pareto World',
                xlabel='height (cm)',
                ylabel='CDF')


def TallestPareto(iters=2, n=10000, xmin=100, alpha=1.7):
    """Find the tallest person in Pareto World.

    iters: how many samples to generate
    n: how many in each sample
    xmin: parameter of the Pareto distribution
    alpha: parameter of the Pareto distribution
    """
    tallest = 0
    for i in range(iters):
        t = [xmin * random.paretovariate(alpha) for i in range(n)]
        tallest = max(max(t), tallest)
    return tallest


def ParetoSample(alpha, xmin, n):
    t = scipy.stats.pareto.rvs(b=alpha, size=n) * xmin
    return t.round()


def main():
    
    counter = Counter()
    for i in range(10000):
        sample = ParetoSample(1.7, 0.001, 10000)
        counter.update(Counter(sample))

    print(len(counter))
    return

    pmf = thinkstats2.Pmf(counter)
    print('mean', pmf.Mean())
    for x, prob in pmf.Largest(10):
        print(x)

    cdf = thinkstats2.Cdf(pmf)
    thinkplot.Cdf(cdf, complement=True)
    thinkplot.Show(xscale='log',
                   yscale='log')
    return

    MakeFigure()
    MakeParetoCdf()
    print(TallestPareto(iters=2))


if __name__ == "__main__":
    main()