Statistical analysis - Part 2¶

#Imports
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import scipy.stats as ss

Representative analytic distributions¶

In this part we will look at how to generate and plot analytic distributions.

The Python code below generates and plots the PDF and CDF of a normal (Gaussian) distribution whose parameters are mu and sigma.

x = np.linspace(-5, 5, 5000)
mu = 0
sigma = 1

y_pdf = ss.norm.pdf(x, mu, sigma) # the normal pdf
y_cdf = ss.norm.cdf(x, mu, sigma) # the normal cdf

plt.plot(x, y_pdf, label='pdf')
plt.plot(x, y_cdf, label='cdf')
plt.legend();

Next we write code to:

Plot the PDF of three exponential distributions (with lambda equal to 0.5, 1, and 2) on the same plot.
Plot the CDF of three exponential distributions (with lambda equal to 0.5, 1, and 2) on the same plot.

xs = np.linspace(0, 5, 5000)

lambdas=[0.5, 1, 2]

for lamb in lambdas:
    y_pdf = [lamb*np.e**(-lamb*x) for x in xs]
    plt.plot(x, y_pdf, label= 'lambda= ' + str(lamb))
    plt.legend()

for lamb in lambdas:
    y_cdf = [1-np.e**(-lamb*x) for x in xs]
    plt.plot(x, y_cdf, label= 'lambda= ' + str(lamb))
    plt.legend()

Empirical distributions vs. analytic distributions¶

The question we are trying to answer in this part is: How well can we model empirical distributions with analytic distributions?

Let's start by asking the question Are the Belmont Stakes results normally distributed?

For context: Since 1926, the Belmont Stakes is a 1.5 mile-long race of 3-year old thoroughbred horses. Secretariat ran the fastest Belmont Stakes in history in 1973. While that was the fastest year, 1970 was the slowest because of unusually wet and sloppy conditions. These two outliers have been removed from the data set, which has been obtained by scraping the Belmont Wikipedia page. (The file belmont.csv is available on Canvas, if you want to learn more about the race's results.)

The code below:

computes mean and standard deviation of Belmont winners' times with the two outliers removed.
takes 10,000 samples out of a normal distribution with this mean and standard deviation using np.random.normal().
computes the CDF of the theoretical samples and the ECDF of the Belmont winners' data, assigning the results to x_theor, y_theor and x, y, respectively.
plots the CDF of your samples with the ECDF, with labeled axes.

import numpy as np
import matplotlib.pyplot as plt

belmont_no_outliers = np.array([148.51,  146.65,  148.52,  150.7,  150.42,  150.88,  151.57,
                                147.54,  149.65,  148.74,  147.86,  148.75,  147.5,  148.26,
                                149.71,  146.56,  151.19,  147.88,  149.16,  148.82,  148.96,
                                152.02,  146.82,  149.97,  146.13,  148.1,  147.2,  146.,
                                146.4,  148.2,  149.8,  147.,  147.2,  147.8,  148.2,
                                149.,  149.8,  148.6,  146.8,  149.6,  149.,  148.2,
                                149.2,  148.,  150.4,  148.8,  147.2,  148.8,  149.6,
                                148.4,  148.4,  150.2,  148.8,  149.2,  149.2,  148.4,
                                150.2,  146.6,  149.8,  149.,  150.8,  148.6,  150.2,
                                149.,  148.6,  150.2,  148.2,  149.4,  150.8,  150.2,
                                152.2,  148.2,  149.2,  151.,  149.6,  149.6,  149.4,
                                148.6,  150.,  150.6,  149.2,  152.6,  152.8,  149.6,
                                151.6,  152.8,  153.2,  152.4,  152.2])

def ecdf(data):
    """Compute ECDF for a one-dimensional array of measurements."""

    # Number of data points: n
    n = len(data)

    # x-data for the ECDF: x
    x = np.sort(data)

    # y-data for the ECDF: y
    y = np.arange(1, n + 1) / n

    return x, y

# Seed random number generator
np.random.seed(42)

# Compute mean and standard deviation: mu, sigma
mu = np.mean(belmont_no_outliers)
sigma = np.std(belmont_no_outliers)

# Sample out of a normal distribution with this mu and sigma: samples
samples = np.random.normal(mu, sigma, 10000)

# Get the CDF of the samples and of the data
x_theor, y_theor = ecdf(samples)
x, y = ecdf(belmont_no_outliers)

# Plot the CDFs and show the plot
_ = plt.plot(x_theor, y_theor)
_ = plt.plot(x, y, marker='.', linestyle='none')
plt.margins(0.02)
_ = plt.xlabel('Belmont winning time (sec.)')
_ = plt.ylabel('CDF')
plt.show()

The empirical distribution of the data looks normally distributed.
Other plots that could be used to represent this data are a histogram, or a normal probability plot, which provides a linear representation of normally distributed data.

Let's try to answer the question: What are the chances of a horse matching or beating Secretariat's record?

Assuming that the Belmont winners' times are Normally distributed (with the 1970 and 1973 years removed), we will attempt to answer the question: What is the probability that the winner of a given Belmont Stakes will run it as fast or faster than Secretariat?

To this end, we will:

Take 1,000,000 samples from the normal distribution using the np.random.normal() function.
Compute the mean mu and standard deviation sigma from the belmont_no_outliers array.
Compute the fraction of samples that have a time less than or equal to Secretariat's time of 144 seconds.
Print the result.

sec_winners=[]

np.random.seed(50)

# Compute mean and standard deviation: mu, sigma
mu = np.mean(belmont_no_outliers)
sigma = np.std(belmont_no_outliers)

# Sample out of a normal distribution with this mu and sigma: samples
samples = np.random.normal(mu, sigma, 1000000)

for sample in samples:
    if sample<=144:
        sec_winners.append(sample)

sec_winners=np.array(sec_winners)
print("Probability of beating or tying Secretariat's record: ", len(sec_winners)/len(samples))

Probability of beating or tying Secretariat's record:  0.000643

Let's investigate whether the speed of light measurements by Michelson are normally distributed.

We will follow a similar sequence of steps as above, namely:

Compute the mean mu and standard deviation sigma from the michelson_speed_of_light array.
Take 10,000 samples out of a normal distribution with this mean and standard deviation using np.random.normal().
Compute the CDF of the theoretical samples and the ECDF of the Michelson speed of light data, assigning the results to x_theor, y_theor and x, y, respectively.
Plot the CDF of your samples with the ECDF, with labeled axes.

For more on Michelson: https://en.wikipedia.org/wiki/Albert_A._Michelson

michelson=pd.read_csv("michelson_speed_of_light.csv")

michelson_speed_of_light=np.array(michelson['velocity of light in air (km/s)'])

# Seed random number generator
np.random.seed(0)

# Compute mean and standard deviation: mu, sigma
mu = np.mean(michelson_speed_of_light)
sigma = np.std(michelson_speed_of_light)

# Sample out of a normal distribution with this mu and sigma: samples
samples = np.random.normal(mu, sigma, 10000)

# Get the CDF of the samples and of the data
x_theor, y_theor = ecdf(samples)
x, y = ecdf(michelson_speed_of_light)

# Plot the CDFs and show the plot
_ = plt.plot(x_theor, y_theor)
_ = plt.plot(x, y, marker='.', linestyle='none')
plt.margins(0.02)
_ = plt.xlabel('Velocity of light in air (km/s)')
_ = plt.ylabel('CDF')
plt.show()

The empirical distribution of the data looks normally distributed.

Next, let's turn our attention to baby births.

In the real world, exponential distributions come up when we look at a series of events and measure the times between events, called interarrival times. If the events are equally likely to occur at any time, the distribution of interarrival times tends to look like an exponential distribution.

We will use the dataset from babies_brisbane.csv containing information about the time of birth for 44 babies born in a hospital in Brisbane, Australia, on December 18, 1997, as reported in the local paper.

We will write code to:

Read the data and build a Pandas dataframe.
Compute the reciprocal of the mean of the sample exponential distribution (call this lam, since lambda is a reserved word in Python).
Take 10,000 samples out of an exponential distribution with this scale using np.random.exponential().
Compute the CDF of the theoretical samples and the ECDF of the sample data, assigning the results to x_theor, y_theor and x, y, respectively.
Plot the CDF of your samples with the ECDF, with labeled axes.
Compute the Complementary CDF (CCDF) and plot the CCDF for both theoretical and sample values, on a log-y scale.

babies_brisbane=pd.read_csv("babies_brisbane.csv")

diffs=babies_brisbane['minutes'].diff()
diffs_mean=np.mean(diffs)
lam=1/diffs_mean

# Seed random number generator
np.random.seed(25)

# Sample out of a normal distribution with this mu and sigma: samples
samples = np.random.exponential(diffs_mean, 10000)

# Get the CDF of the samples and of the data
x_theor, y_theor = ecdf(samples)
x, y = ecdf(diffs)

# Plot the CDFs and show the plot
_ = plt.plot(x_theor, y_theor)
_ = plt.plot(x, y, marker='.', linestyle='none')
plt.margins(0.02)
_ = plt.xlabel('Time between births (min)')
_ = plt.ylabel('CDF')
plt.show()

y_c=1-y
y_theor_c=1-y_theor

_=plt.yscale('log')
_ = plt.plot(x_theor, y_theor_c)
_ = plt.plot(x, y_c, marker='.', linestyle='none')
_ = plt.xlabel('Time between births (min)')
_ = plt.ylabel('CCDF')
plt.show()

Moments and skewness¶

Let's revisit the dataset of salaries from the data science stack notebook and use it to measure skewness.

salaries = pd.read_csv('salaries.csv')
salaries.describe()

We write Python code to:

Compute the median and mean salary for the entire sample.
Compute the first raw moment and show that it is equivalent to computing the mean value.
Compute the second central moment and show that it is equivalent to computing the variance.
Compute the skewness using scipy.stats.skew

salary_act=np.array(salaries['earn'])
salary_mean=np.mean(salary_act)
salary_median=np.median(salary_act)
print("Mean of salaries: \n", salary_mean)
print("Median of salaries: ", salary_median)

Mean of salaries: 
 23154.773489932886
Median of salaries:  20000.0

moment_1=sum(salary_act)/len(salary_act)

print("First raw moment: \n", moment_1)
print("Salary mean: ", salary_mean)

First raw moment: 
 23154.773489932886
Salary mean:  23154.773489932886

moment_2=ss.moment(salary_act, moment=2)
salary_var=np.std(salary_act)**2

print("Second central moment: \n", moment_2)
print("Salary variance: ", salary_var)

Second central moment: 
 378852251.6248677
Salary variance:  378852251.6248677

salary_skew=ss.skew(salary_act)
print("Salary skewness: ", salary_skew)

Salary skewness:  2.880309741267592

The skewness value above means that there are higher values in the right tail of the distribution of this data.

Next we write code to repeat the steps above, this time for male and female employees separately.

male_salaries=[]
female_salaries=[]

for i in range(0, len(salaries)):
    if salaries['sex'][i]=='male':
        male_salaries.append(salaries['earn'][i])
    else:
        female_salaries.append(salaries['earn'][i])

male_salaries=np.array(male_salaries)
female_salaries=np.array(female_salaries)

print("Average salary (male): \n", np.mean(male_salaries))
print("Median  salary (male): \n", np.median(male_salaries))

print("Average salary (female): \n", np.mean(female_salaries))
print("Median  salary (female):", np.median(female_salaries))

Average salary (male): 
 29786.130693069306
Median  salary (male): 
 25000.0
Average salary (female): 
 18280.195050946142
Median  salary (female): 15000.0

moment_1_m=sum(male_salaries)/len(male_salaries)
moment_1_f=sum(female_salaries)/len(female_salaries)

print("First raw moment (male): \n", moment_1_m)
print("Salary mean (male): \n", np.mean(male_salaries))
print("First raw moment (female): \n", moment_1_f)
print("Salary mean (female): ", np.mean(female_salaries))

First raw moment (male): 
 29786.130693069306
Salary mean (male): 
 29786.130693069306
First raw moment (female): 
 18280.195050946142
Salary mean (female):  18280.195050946142

moment_2_m=ss.moment(male_salaries, moment=2)
moment_2_f=ss.moment(female_salaries, moment=2)

print("Second central moment (male): \n", moment_2_m)
print("Salary variance (male): \n", np.std(male_salaries)**2)
print("Second central moment (female): \n", moment_2_f)
print("Salary variance (female): ", np.std(female_salaries)**2)

Second central moment (male): 
 541713403.7611371
Salary variance (male): 
 541713403.7611371
Second central moment (female): 
 203049680.6286218
Salary variance (female):  203049680.62862176

salary_skew_m=ss.skew(male_salaries)
salary_skew_f=ss.skew(female_salaries)

print("Salary skew (male): \n", salary_skew_m)
print("Salary skew (female): \n", salary_skew_f)

Salary skew (male): 
 2.733752806565786
Salary skew (female): 
 2.234266764847592

	earn	height	ed	age
count	1192.000000	1192.000000	1192.000000	1192.000000
mean	23154.773490	66.915154	13.504195	41.378356
std	19472.296925	3.853968	2.420175	15.867428
min	200.000000	57.503219	3.000000	18.000000
25%	10000.000000	64.009746	12.000000	29.000000
50%	20000.000000	66.451265	13.000000	38.000000
75%	30000.000000	69.848100	16.000000	51.000000
max	200000.000000	77.051282	18.000000	91.000000