| Download

📚 The CoCalc Library - books, templates and other resources

Project: 📚 The Library - Shared Public Version

Path: cocalc-examples / think-stats-2ed / homeworks / homework02.ipynb

Views: ⁹⁶¹⁶⁰
License: OTHER

Kernel: Python 3

Homework 2

Visualize, describe, and model distributions

Allen Downey

MIT License

In [1]:

%matplotlib inline

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='white')

from utils import decorate
from thinkstats2 import Pmf, Cdf

import thinkstats2
import thinkplot

Here are some of the functions from Chapter 5.

In [2]:

def MakeNormalModel(values, label=''):
    """Plots a CDF with a Normal model.

    values: sequence
    """
    cdf = thinkstats2.Cdf(values, label=label)

    mean, var = thinkstats2.TrimmedMeanVar(values)
    std = np.sqrt(var)
    print('n, mean, std', len(values), mean, std)

    xmin = mean - 4 * std
    xmax = mean + 4 * std

    xs, ps = thinkstats2.RenderNormalCdf(mean, std, xmin, xmax)
    thinkplot.Plot(xs, ps, label='model', linewidth=4, color='0.8')
    thinkplot.Cdf(cdf)

In [3]:

def MakeNormalPlot(values, label=''):
    """Generates a normal probability plot.

    values: sequence
    """
    mean, var = thinkstats2.TrimmedMeanVar(values, p=0.01)
    std = np.sqrt(var)

    xs = [-5, 5]
    xs, ys = thinkstats2.FitLine(xs, mean, std)
    thinkplot.Plot(xs, ys, color='0.8', label='model')

    xs, ys = thinkstats2.NormalProbability(values)
    thinkplot.Plot(xs, ys, '+', alpha=0.3, label=label)

Read the GSS data again.

In [4]:

%time gss = pd.read_hdf('gss.hdf5', 'gss')
gss.shape

In [5]:

gss.head()

Most variables use special codes to indicate missing data. We have to be careful not to use these codes as numerical data; one way to manage that is to replace them with NaN, which Pandas recognizes as a missing value.

In [6]:

def replace_invalid(df):
    df.realinc.replace([0], np.nan, inplace=True)                  
    df.educ.replace([98,99], np.nan, inplace=True)
    # 89 means 89 or older
    df.age.replace([98, 99], np.nan, inplace=True) 
    df.cohort.replace([9999], np.nan, inplace=True)
    df.adults.replace([9], np.nan, inplace=True)

replace_invalid(gss)

Distribution of age

Here's the CDF of ages.

In [7]:

cdf_age = Cdf(gss.age)
thinkplot.cdf(cdf_age, label='age')

decorate(title='Distribution of age', 
         xlabel='Age (years)', 
         ylabel='CDF')

Exercise: Each of the following cells shows the distribution of ages under various transforms, compared to various models. In each text cell, add a sentence or two that interprets the result. What can we say about the distribution of ages based on each figure?

Here's the CDF of ages compared to a normal distribution with the same mean and standard deviation.

Interpretation:

In [8]:

MakeNormalModel(gss.age.dropna(), label='')

decorate(title='Distribution of age', 
         xlabel='Age (years)', 
         ylabel='CDF')

Here's a normal probability plot for the distribution of ages.

Interpretation:

In [9]:

MakeNormalPlot(gss.age.dropna(), label='')

decorate(title='Normal probability plot', 
         xlabel='Standard normal sample', 
         ylabel='Age (years)')

Here's the complementary CDF on a log-y scale.

Interpretation:

In [10]:

thinkplot.cdf(cdf_age, label='age', complement=True)

decorate(title='Distribution of age', 
         xlabel='Age (years)', 
         ylabel='Complementary CDF, log scale',
         yscale='log')

Here's the CDF of ages on a log-x scale.

Interpretation:

In [11]:

thinkplot.cdf(cdf_age, label='age')

decorate(title='Distribution of age', 
         xlabel='Age (years)', 
         ylabel='CDF',
         xscale='log')

Here's the CDF of the logarithm of ages, compared to a normal model.

Interpretation:

In [12]:

values = np.log10(gss.age.dropna())
MakeNormalModel(values, label='')

decorate(title='Distribution of age', 
         xlabel='Age (log10 years)', 
         ylabel='CDF')

Here's a normal probability plot for the logarithm of ages.

Interpretation:

In [13]:

MakeNormalPlot(values, label='')

decorate(title='Distribution of age', 
         xlabel='Standard normal sample', 
         ylabel='Age (log10 years)')

Here's the complementary CDF on a log-log scale.

Interpretation:

In [14]:

thinkplot.cdf(cdf_age, label='age', complement=True)

decorate(title='Distribution of age', 
         xlabel='Age (years)', 
         ylabel='Complementary CDF, log scale',
         xscale='log',
         yscale='log')

Here's a test to see whether ages are well-modeled by a Weibull distribution.

Interpretation:

In [15]:

thinkplot.cdf(cdf_age, label='age', transform='Weibull')

decorate(title='Distribution of age', 
         xlabel='Age (years)', 
         ylabel='log Complementary CDF, log scale',
         xscale='log',
         yscale='log')

Distribution of income

Here's the CDF of realinc.

In [16]:

cdf_realinc = Cdf(gss.realinc)
thinkplot.cdf(cdf_realinc, label='income')

decorate(title='Distribution of income', 
         xlabel='Income (1986 $)', 
         ylabel='CDF')

Exercise: Use visualizations like the ones in the previous exercise to see whether there is an analytic model that describes the distribution of gss.realinc well.

In [17]:

# Solution goes here

Here's a normal probability plot for the values.

In [18]:

# Solution goes here

Here's the complementary CDF on a log-y scale.

In [19]:

# Solution goes here

Here's the CDF on a log-x scale.

In [20]:

# Solution goes here

Here's the CDF of the logarithm of the values, compared to a normal model.

In [21]:

# Solution goes here

Here's a normal probability plot for the logarithm of the values.

In [22]:

# Solution goes here

Here's the complementary CDF on a log-log scale.

In [23]:

# Solution goes here

Here's a test to see whether the values are well-modeled by a Weibull distribution.

Interpretation:

In [24]:

# Solution goes here

BRFSS

In [25]:

%time brfss = pd.read_hdf('brfss.hdf5', 'brfss')
brfss.head()

Let's look at the distribution of height in the BRFSS dataset. Here's the CDF.

In [26]:

heights = brfss.HTM4

cdf_heights = Cdf(heights)
thinkplot.Cdf(cdf_heights)

decorate(xlabel='Height (cm)', ylabel='CDF')

To see whether a normal model describes this data well, we can use KDE to estimate the PDF.

In [27]:

from scipy.stats import gaussian_kde

Here's an example using the default bandwidth method.

In [28]:

kde = gaussian_kde(heights.dropna())

xs = np.linspace(heights.min(), heights.max())
ds = kde.evaluate(xs)
ds /= ds.sum()

plt.plot(xs, ds, label='KDE heights')

decorate(xlabel='Height (cm)', ylabel='PDF')

It doesn't work very well; we can improve it by overriding the bandwidth with a constant.

In [29]:

kde = gaussian_kde(heights.dropna(), bw_method=0.3)

ds = kde.evaluate(xs)
ds /= ds.sum()

plt.plot(xs, ds, label='KDE heights')

decorate(xlabel='Height (cm)', ylabel='PDF')

Now we can generate a normal model with the same mean and standard deviation.

In [30]:

mean = heights.mean()
std = heights.std()

mean, std

Here's the model compared to the estimated PDF.

In [31]:

normal_pdf = thinkstats2.NormalPdf(mean, std)

ps = normal_pdf.Density(xs)
ps /= ps.sum()

plt.plot(xs, ps, color='gray', label='Normal model')
plt.plot(xs, ds, label='KDE heights')

decorate(xlabel='Height (cm)', ylabel='PDF')

The data don't fit the model particularly well, possibly because the distribution of heights is a mixture of two distributions, for men and women.

Exercise: Generate a similar figure for just women's heights and see if the normal model does any better.

In [32]:

# Solution goes here

In [33]:

# Solution goes here

Exercise: Generate a similar figure for men's weights, brfss.WTKG3. How well does the normal model fit?

In [34]:

# Solution goes here

In [35]:

# Solution goes here

Exercise: Try it one more time with the log of men's weights. How well does the normal model fit? What does that imply about the distribution of weight?

In [36]:

# Solution goes here

In [37]:

# Solution goes here

Skewness

Let's look at the skewness of the distribution of weights for men and women.

In [38]:

male = (brfss.SEX == 1)
male_weights = brfss.loc[male, 'WTKG3']

In [39]:

female = (brfss.SEX == 2)
female_weights = brfss.loc[female, 'WTKG3']

As we've seen, these distributions are skewed to the right, so we expect the mean to be higher than the median.

In [40]:

male_weights.mean(), male_weights.median()

We can compute the moment-based sample skewness using Pandas or thinkstats2. The results are almost the same.

In [41]:

male_weights.skew(), thinkstats2.Skewness(male_weights.dropna())

But moment-based sample skewness is a terrible statistic! A more robust alternative is Pearson's median skewness:

In [42]:

thinkstats2.PearsonMedianSkewness(male_weights.dropna())

Exercise: Compute the same statistics for women. Which distribution is more skewed?

In [43]:

# Solution goes here

In [44]:

# Solution goes here

In [45]:

# Solution goes here

Exercise: Explore the GSS or BRFSS dataset and find something interesting!

In [ ]:

Homework 2

Distribution of age

Distribution of income

BRFSS

Skewness

Product

Resources

Company