Exercise from Think Stats, 2nd Edition (thinkstats2.com)
Allen Downey


In [106]:
from __future__ import print_function, division

import thinkstats2
import thinkplot
from brfss import *
import populations as p
import random
import pandas as pd
import test_models

%matplotlib inline

Exercise 5.1

In the BRFSS (see Section 5.4), the distribution of heights is roughly normal with parameters µ = 178 cm and σ = 7.7 cm for men, and µ = 163 cm and σ = 7.3 cm for women.

In order to join Blue Man Group, you have to be male between 5’10” and 6’1” (see http://bluemancasting.com). What percentage of the U.S. male population is in this range? Hint: use scipy.stats.norm.cdf.

scipy.stats contains objects that represent analytic distributions


In [2]:
import scipy.stats

For example scipy.stats.norm represents a normal distribution.


In [3]:
mu = 178
sigma = 7.7
dist = scipy.stats.norm(loc=mu, scale=sigma)
type(dist)


Out[3]:
scipy.stats._distn_infrastructure.rv_frozen

A "frozen random variable" can compute its mean and standard deviation.


In [3]:
dist.mean(), dist.std()


Out[3]:
(178.0, 7.7000000000000002)

It can also evaluate its CDF. How many people are more than one standard deviation below the mean? About 16%


In [4]:
dist.cdf(mu-sigma)


Out[4]:
0.15865525393145741

How many people are between 5'10" and 6'1"?


In [6]:
dist.cdf(185.42) - dist.cdf(177.8)


Out[6]:
0.34274683763147368

In [32]:
thinkstats2.RandomSeed(17)

nrows = int(1000)    
df = brfss.ReadBrfss(nrows=10000)

In [36]:
MakeNormalPlot(df.age)



In [39]:
p.MakeFigures()


Number of cities/towns 19515
Writing populations_pareto.pdf
Writing populations_pareto.eps
Writing populations_normal.pdf
Writing populations_normal.eps
<matplotlib.figure.Figure at 0x10a23ba50>

Exercise 5.2

To get a feel for the Pareto distribution, let’s see how different the world would be if the distribution of human height were Pareto. With the parameters $x_m = 1$ m and $α = 1.7$, we get a distribution with a reasonable minimum, 1 m, and median, 1.5 m.

Plot this distribution. What is the mean human height in Pareto world? What fraction of the population is shorter than the mean? If there are 7 billion people in Pareto world, how many do we expect to be taller than 1 km? How tall do we expect the tallest person to be?

scipy.stats.pareto represents a pareto distribution. In Pareto world, the distribution of human heights has parameters alpha=1.7 and xmin=1 meter. So the shortest person is 100 cm and the median is 150.


In [40]:
alpha = 1.7
xmin = 1
dist = scipy.stats.pareto(b=alpha, scale=xmin)
dist.median()


Out[40]:
1.5034066538560549

In [41]:
xs, ps = thinkstats2.RenderParetoCdf(xmin, alpha, 0, 10.0, n=100) 
thinkplot.Plot(xs, ps, label=r'$\alpha=%g$' % alpha)
thinkplot.Config(xlabel='height (m)', ylabel='CDF')


What is the mean height in Pareto world?


In [42]:
dist.mean()


Out[42]:
2.4285714285714288

What fraction of people are shorter than the mean?


In [43]:
dist.cdf(dist.mean())


Out[43]:
0.77873969756528805

Out of 7 billion people, how many do we expect to be taller than 1 km? You could use dist.cdf or dist.sf.


In [44]:
(1 - dist.cdf(1000))*7000000000


Out[44]:
55602.976430479954

How tall do we expect the tallest person to be? Hint: find the height that yields about 1 person.


In [56]:
dist.isf(1/7000000000)


Out[56]:
618349.61067595053

Exercise 5.3

The Weibull distribution is a generalization of the exponential distribution that comes up in failure analysis (see http://wikipedia.org/wiki/Weibull_distribution). Its CDF is

$CDF(x) = 1 − \exp(−(x / λ)^k)$

Can you find a transformation that makes a Weibull distribution look like a straight line? What do the slope and intercept of the line indicate?

Use random.weibullvariate to generate a sample from a Weibull distribution and use it to test your transformation.


In [98]:
alpha = 100
lam = 1
sample = [random.weibullvariate(alpha, lam) for i in xrange(1000)]

In [99]:
cdf = thinkstats2.Cdf(sample)
thinkplot.Cdf(np.log(np.log(cdf)), complement=True)
thinkplot.Show()


---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-99-e2895d30577b> in <module>()
      1 cdf = thinkstats2.Cdf(sample)
----> 2 thinkplot.Cdf(np.log(np.log(cdf)), complement=True)
      3 thinkplot.Show()

/Users/nathankiner/Documents/Data Science/ThinkStats2/code/thinkstats2.pyc in __getitem__(self, x)
   1004 
   1005     def __getitem__(self, x):
-> 1006         return self.Prob(x)
   1007 
   1008     def __setitem__(self):

/Users/nathankiner/Documents/Data Science/ThinkStats2/code/thinkstats2.pyc in Prob(self, x)
   1076         if x < self.xs[0]:
   1077             return 0.0
-> 1078         index = bisect.bisect(self.xs, x)
   1079         p = self.ps[index-1]
   1080         return p

KeyboardInterrupt: 

Exercise 5.4

For small values of n, we don’t expect an empirical distribution to fit an analytic distribution exactly. One way to evaluate the quality of fit is to generate a sample from an analytic distribution and see how well it matches the data.

For example, in Section 5.1 we plotted the distribution of time between births and saw that it is approximately exponential. But the distribution is based on only 44 data points. To see whether the data might have come from an exponential distribution, generate 44 values from an exponential distribution with the same mean as the data, about 33 minutes between births.

Plot the distribution of the random values and compare it to the actual distribution. You can use random.expovariate to generate the values.


In [14]:
import analytic

df = analytic.ReadBabyBoom()
diffs = df.minutes.diff()
cdf = thinkstats2.Cdf(diffs, label='actual')
thinkplot.Cdf(cdf, complement=True)
thinkplot.Config(yscale='log')



In [103]:
sample = [random.expovariate(1/33) for i in xrange(44)]

In [105]:
cdf = thinkstats2.Cdf(sample)
thinkplot.Cdf(cdf, complement=True)
thinkplot.Config(yscale='log')



In [117]:
test_models.main("test_models.py", "mystery2.dat")


<matplotlib.figure.Figure at 0x10d892410>

In [ ]: