This article is based on Statistics tutorial to learn essential concepts of Statistics, that we need in Data Science. I will try to present the concepts in a fun and interactive way and I encourage you to play with the code to get a better grasp of the concepts.
Download the data set
Lets start with importing the libraries and reading the data set
# Dependencies
# Standard Dependencies
import os
import numpy as np
import pandas as pd
from math import sqrt
# Visualization
from pylab import *
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import seaborn as sns
# Statistics
from statistics import median
from scipy import signal
from scipy.special import factorial
import scipy.stats as stats
from scipy.stats import sem, binom, lognorm, poisson, bernoulli, spearmanr
from scipy.fftpack import fft, fftshift
# Scikit-learn for Machine Learning models
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
# Seed for reproducability
seed = 12345
np.random.seed(seed)
from google.colab import files
uploaded = files.upload()
df = pd.read_csv("toy_dataset.csv")Discrete and Continuous Variables
- A discrete variable is a variable that can only take on a countable number of values. If you can count a set of items, then it’s a discrete variable. An example of a discrete variable is the outcome of a dice.
It can only have 1 of 6 different possible outcomes and is therefore discrete. A discrete random variable can have an infinite number of values. For example, the whole set of natural numbers (1,2,3,etc.) is countable and therefore discrete. - A continuous variable takes on an uncountable number of values. An example of a continuous variable is length. Length can be measured to an arbitrary degree and is therefore continuous.
- In statistics we represent a distribution of discrete variables with PMF’s (Probability Mass Functions) and CDF’s (Cumulative Distribution Functions). We represent distributions of continuous variables with PDF’s (Probability Density Functions) and CDF’s.
- The PMF defines the probability of all possible values x of the random variable. A PDF is the same but for continuous values. The CDF represents the probability that the random variable X will have an outcome less or equal to the value x. The name CDF is used for both discrete and continuous distributions.
- The functions that describe PMF’s, PDF’s and CDF’s can be quite daunting at first, but their visual counterpart often looks quite intuitive.
PMF (Probability Mass Function)
Here we visualize a PMF of a binomial distribution. You can see that the possible values are all integers. For example, no values are between 50 and 51.
The PMF of a binomial distribution in function form:
# PMF Visualization
n = 100
p = 0.5
fig, ax = plt.subplots(1, 1, figsize=(17,5))
x = np.arange(binom.ppf(0.01, n, p), binom.ppf(0.99, n, p))
ax.plot(x, binom.pmf(x, n, p), 'bo', ms=8, label='Binomial PMF')
ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=5, alpha=0.5)
rv = binom(n, p)
#ax.vlines(x, 0, rv.pmf(x), colors='k', linestyles='-', lw=1, label='frozen PMF')
ax.legend(loc='best', frameon=False, fontsize='xx-large')
plt.title('PMF of a binomial distribution (n=100, p=0.5)', fontsize='xx-large')
plt.show()
PDF (Probability Density Functions)
The PDF is the same as a PMF, but continuous. It can be said that the distribution has an infinite number of possible values. Here we visualize a simple normal distribution with a mean of 0 and standard deviation of 1.
PDF of a normal distribution in formula form:
# Plot normal distribution
mu = 0
variance = 1
sigma = sqrt(variance)
x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)
plt.figure(figsize=(16,5))
plt.plot(x, stats.norm.pdf(x, mu, sigma), label='Normal Distribution')
plt.title('Normal Distribution with mean = 0 and std = 1')
plt.legend(fontsize='xx-large')
plt.show()
CDF (Cumulative Distribution Function)
The CDF maps the probability that a random variable X will take a value of less than or equal to a value x (P(X ≤ x)). CDF’s can be discrete or continuous. In this section we visualize the continuous case. You can see in the plot that the CDF accumulates all probabilities and is therefore bounded between 0 ≤ x ≤ 1.
The CDF of a normal distribution as a formula:
# Data
X = np.arange(-2, 2, 0.01)
Y = exp(-X ** 2)
# Normalize data
Y = Y / (0.01 * Y).sum()
# Plot the PDF and CDF
plt.figure(figsize=(15,5))
plt.title('Continuous Normal Distributions', fontsize='xx-large')
plot(X, Y, label='Probability Density Function (PDF)')
plot(X, np.cumsum(Y * 0.01), 'r', label='Cumulative Distribution Function (CDF)')
plt.legend(fontsize='xx-large')
plt.show()
Distributions
- A Probability distribution tells us something about the likelihood of each value of the random variable.
- A random variable X is a function that maps events to real numbers.
- The visualizations in this section are of discrete distributions. Many of these distributions can however also be continuous.
Uniform Distribution
A Uniform distribution is pretty straightforward. Every value has an equal change of occuring. Therefore, the distribution consists of random values with no patterns in them. In this example we generate random floating numbers between 0 and 1.
The PDF of a Uniform Distribution:
CDF:
# Uniform distribution (between 0 and 1)
uniform_dist = np.random.random(1000)
uniform_df = pd.DataFrame({'value' : uniform_dist})
uniform_dist = pd.Series(uniform_dist)
plt.figure(figsize=(18,5))
sns.scatterplot(data=uniform_df)
plt.legend(fontsize='xx-large')
plt.title('Scatterplot of a Random/Uniform Distribution', fontsize='xx-large')
plt.figure(figsize=(18,5))
sns.distplot(uniform_df)
plt.title('Random/Uniform distribution', fontsize='xx-large')
Normal Distribution
A normal distribution (also called Gaussian or Bell Curve) is very common and convenient. This is mainly because of the Central Limit Theorem (CLT), which states that as the amount independent random samples (like multiple coin flips) goes to infinity the distribution of the sample mean tends towards a normal distribution.
PDF of a normal distribution:
CDF:
# Generate Normal Distribution
normal_dist = np.random.randn(10000)
normal_df = pd.DataFrame({'value' : normal_dist})
# Create a Pandas Series for easy sample function
normal_dist = pd.Series(normal_dist)
normal_dist2 = np.random.randn(10000)
normal_df2 = pd.DataFrame({'value' : normal_dist2})
# Create a Pandas Series for easy sample function
normal_dist2 = pd.Series(normal_dist)
normal_df_total = pd.DataFrame({'value1' : normal_dist,
'value2' : normal_dist2})# Scatterplot
plt.figure(figsize=(18,5))
sns.scatterplot(data=normal_df)
plt.legend(fontsize='xx-large')
plt.title('Scatterplot of a Normal Distribution', fontsize='xx-large')
# Normal Distribution as a Bell Curve
plt.figure(figsize=(18,5))
sns.distplot(normal_df)
plt.title('Normal distribution (n=1000)', fontsize='xx-large')
Binomial Distribution
A Binomial Distribution has a countable number of outcomes and is therefore discrete.
Binomial distributions must meet the following three criteria:
- The number of observations or trials is fixed. In other words, you can only figure out the probability of something happening if you do it a certain number of times.
- Each observation or trial is independent. In other words, none of your trials have an effect on the probability of the next trial.
- The probability of success is exactly the same from one trial to another.
An intuitive explanation of a binomial distribution is flipping a coin 10 times. If we have a fair coin our chance of getting heads (p) is 0.50. Now we throw the coin 10 times and count how many times it comes up heads. In most situations we will get heads 5 times, but there is also a change that we get heads 9 times. The PMF of a binomial distribution will give these probabilities if we say N = 10 and p = 0.5. We say that the x for heads is 1 and 0 for tails.
PMF:
CDF:
A Bernoulli Distribution is a special case of a Binomial Distribution.
All values in a Bernoulli Distribution are either 0 or 1.
For example, if we take an unfair coin which falls on heads 60 % of the time, we can describe the Bernoulli distribution as follows:
p (change of heads) = 0.6
1 – p (change of tails) = 0.4
heads = 1
tails = 0
Formally, we can describe a Bernoulli distribution with the following PMF (Probability Mass Function):
# Change of heads (outcome 1)
p = 0.6
# Create Bernoulli samples
bern_dist = bernoulli.rvs(p, size=1000)
bern_df = pd.DataFrame({'value' : bern_dist})
bern_values = bern_df['value'].value_counts()
# Plot Distribution
plt.figure(figsize=(18,4))
bern_values.plot(kind='bar', rot=0)
plt.annotate(xy=(0.85,300),
s='Samples that came up Tails\nn = {}'.format(bern_values[0]),
fontsize='large',
color='white')
plt.annotate(xy=(-0.2,300),
s='Samples that came up Heads\nn = {}'.format(bern_values[1]),
fontsize='large',
color='white')
plt.title('Bernoulli Distribution: p = 0.6, n = 1000')
Poisson Distribution
The Poisson distribution is a discrete distribution and is popular for modelling the number of times an event occurs in an interval of time or space.
It takes a value lambda, which is equal to the mean of the distribution.
PMF:
CDF:Â 
x = np.arange(0, 20, 0.1)
y = np.exp(-5)*np.power(5, x)/factorial(x)
plt.figure(figsize=(15,8))
plt.title('Poisson distribution with lambda=5', fontsize='xx-large')
plt.plot(x, y, 'bs')
plt.show()
Log-Normal Distribution
A log-normal distribution is continuous. The main characteristic of a log-normal distribution is that it’s logarithm is normally distributed. It is also referred to as Galton’s distribution.
PDF:
CDF:
Where Phi is the CDF of the standard normal distribution.
# Specify standard deviation and mean
std = 1
mean = 5
# Create log-normal distribution
dist=lognorm(std,loc=mean)
x=np.linspace(0,15,200)
# Visualize log-normal distribution
plt.figure(figsize=(15,6))
plt.xlim(5, 10)
plt.plot(x,dist.pdf(x), label='log-normal PDF')
plt.plot(x,dist.cdf(x), label='log-normal CDF')
plt.legend(fontsize='xx-large')
plt.title('Visualization of log-normal PDF and CDF', fontsize='xx-large')
plt.show()
Summary Statistics and Moments
Mean, Median and Mode
Note: The mean is also called the first moment.
Moments
A moment is a quantitative measure that says something about the shape of a distribution. There are central moments and non-central moments. This section is focused on the central moments.
The 0th central moment is the total probability and is always equal to 1.
The 1st moment is the mean (expected value).
The 2nd central moment is the variance.
Variance = The average of the squared distance of the mean. Variance is interesting in a mathematical sense, but the standard deviation is often a much better measure of how spread out the distribution is.
Standard Deviation = The square root of the variance
The 3rd central moment is the skewness.
Skewness = A measure that describes the contrast of one tail versus the other tail. For example, if there are more high values in your distribution than low values then your distribution is ‘skewed’ towards the high values.
The 4th central moment is the kurtosis.
Kurtosis = A measure of how ‘fat’ the tails in the distribution are.
The higher the moment, the harder it is to estimate with samples. Larger samples are required in order to obtain good estimates.
# Summary
print('Summary Statistics for a normal distribution: ')
# Median
medi = median(normal_dist)
print('Median: ', medi)
display(normal_df.describe())
# Standard Deviation
std = sqrt(np.var(normal_dist))
print('The first four calculated moments of a normal distribution: ')
# Mean
mean = normal_dist.mean()
print('Mean: ', mean)
# Variance
var = np.var(normal_dist)
print('Variance: ', var)
# Return unbiased skew normalized by N-1
skew = normal_df['value'].skew()
print('Skewness: ', skew)
# Return unbiased kurtosis over requested axis using Fisher's definition of kurtosis
# (kurtosis of normal == 0.0) normalized by N-1
kurt = normal_df['value'].kurtosis()
print('Kurtosis: ', kurt)
Bias, MSE and SE
Bias is a measure of how far the sample mean deviates from the population mean. The sample mean is also called Expected value.
Formula for Bias:
The formula for expected value (EV) makes it apparent that the bias can also be formulated as the expected value minus the population mean:
# Take sample
normal_df_sample = normal_df.sample(100)
# Calculate Expected Value (EV), population mean and bias
ev = normal_df_sample.mean()[0]
pop_mean = normal_df.mean()[0]
bias = ev - pop_mean
print('Sample mean (Expected Value): ', ev)
print('Population mean: ', pop_mean)
print('Bias: ', bias)#Output Sample mean (Expected Value): -0.09989709931399528 Population mean: -0.012836443467010622 Bias: -0.08706065584698465
MSE (Mean Squared Error) is a formula to measure how much estimators deviate from the true distribution. This can be very useful with for example, evaluating regression models.
RMSE (Root Mean Squared Error) is just the root of the MSE.
from math import sqrt
Y = 100 # Actual Value
YH = 94 # Predicted Value
# MSE Formula
def MSE(Y, YH):
return np.square(YH - Y).mean()
# RMSE formula
def RMSE(Y, YH):
return sqrt(np.square(YH - Y).mean())
print('MSE: ', MSE(Y, YH))
print('RMSE: ', RMSE(Y, YH))#Output MSE: 36.0 RMSE: 6.0
The Standard Error (SE) measures how spread the distribution is from the sample mean.
The formula can also be defined as the standard deviation divided by the square root of the number of samples.
# Standard Error (SE)
uni_sample = uniform_dist.sample(100)
norm_sample = normal_dist.sample(100)
print('Standard Error of uniform sample: ', sem(uni_sample))
print('Standard Error of normal sample: ', sem(norm_sample))
# The random samples from the normal distribution should have a higher standard error#Output Standard Error of uniform sample: 0.028951451192463368 Standard Error of normal sample: 0.10887812225934831
Sampling methods
Non-Representative Sampling:
- Convenience Sampling = Pick samples that are most convenient, like the top of a shelf or people that can be easily approached.
- Haphazard Sampling = Pick samples without thinking about it. This often gives the illusion take you are picking out samples at random.
- Purposive Sampling = Pick samples for a specific purpose. An example is to focus on extreme cases. This can be useful but is limited because it doesn’t allow you to make statements about the whole population.
Representative Sampling:
- Simple Random Sampling = Pick samples (psuedo)randomly.
- Systematic Sampling = Pick samples with a fixed interval. For example every 10th sample (0, 10, 20, etc.).
- Stratified Sampling = Pick the same amount of samples from different groups (strata) in the population.
- Cluster Sampling = Divide the population into groups (clusters) and pick samples from those groups.
# Note that we take very small samples just to illustrate the different sampling methods
print('---Non-Representative samples:---\n')
# Convenience samples
con_samples = normal_dist[0:5]
print('Convenience samples:\n\n{}\n'.format(con_samples))
# Haphazard samples (Picking out some numbers)
hap_samples = [normal_dist[12], normal_dist[55], normal_dist[582], normal_dist[821], normal_dist[999]]
print('Haphazard samples:\n\n{}\n'.format(hap_samples))
# Purposive samples (Pick samples for a specific purpose)
# In this example we pick the 5 highest values in our distribution
purp_samples = normal_dist.nlargest(n=5)
print('Purposive samples:\n\n{}\n'.format(purp_samples))
print('---Representative samples:---\n')
# Simple (pseudo)random sample
rand_samples = normal_dist.sample(5)
print('Random samples:\n\n{}\n'.format(rand_samples))
# Systematic sample (Every 2000th value)
sys_samples = normal_dist[normal_dist.index % 2000 == 0]
print('Systematic samples:\n\n{}\n'.format(sys_samples))
# Stratified Sampling
# We will get 1 person from every city in the dataset
# We have 8 cities so that makes a total of 8 samples
df = pd.read_csv('toy_dataset.csv')
strat_samples = []
for city in df['City'].unique():
samp = df[df['City'] == city].sample(1)
strat_samples.append(samp['Income'].item())
print('Stratified samples:\n\n{}\n'.format(strat_samples))
# Cluster Sampling
# Make random clusters of ten people (Here with replacement)
c1 = normal_dist.sample(10)
c2 = normal_dist.sample(10)
c3 = normal_dist.sample(10)
c4 = normal_dist.sample(10)
c5 = normal_dist.sample(10)
# Take sample from every cluster (with replacement)
clusters = [c1,c2,c3,c4,c5]
cluster_samples = []
for c in clusters:
clus_samp = c.sample(1)
cluster_samples.extend(clus_samp)
print('Cluster samples:\n\n{}'.format(cluster_samples)) #Output ---Non-Representative samples:--- Convenience samples: 0 0.652346 1 0.871600 2 0.268216 3 0.947681 4 0.147268 dtype: float64 Haphazard samples: [-0.5010027615175823, -0.5066784173896016, 0.6814697928519314, 0.8359965545029714, -0.5374734651177141] Purposive samples: 3378 3.525865 2792 3.366626 490 3.260383 7169 3.189940 7057 3.082067 dtype: float64 ---Representative samples:--- Random samples: 3442 0.236063 9169 -0.042718 8511 -1.529128 6500 -0.787297 2141 1.134073 dtype: float64 Systematic samples: 0 0.652346 2000 -0.914345 4000 -0.198577 6000 -0.030696 8000 -1.084444 dtype: float64 Stratified samples: [29892.0, 110450.0, 85635.0, 122364.0, 83214.0, 74469.0, 113187.0, 92854.0] Cluster samples: [-1.1500083079696, 1.2888783369868333, 1.2540511669530716, -0.6914086007291129, 2.259205774351393]
Covariance
- Covariance is a measure of how much two random variables vary together.
- Variance is similar to covariance in that variance shows you how much one variable varies.
- If two variables are independent, their covariance is 0. However, a covariance of 0 does not imply that the variables are independent.
# Covariance between Age and Income
print('Covariance between Age and Income: ')
df[['Age', 'Income']].cov()#Output
Covariance between Age and Income:
Age Income
Age 133.922426 -3.811863e+02
Income -381.186341 6.244752e+08Correlation
Correlation is a standardized version of covariance. Here it becomes more clear that Age and Income do not have a strong correlation in our dataset.
Formula for Pearson’s correlation coefficient:
# Correlation between two normal distributions
# Using Pearson's correlation
print('Pearson: ')
df[['Age', 'Income']].corr(method='pearson')#Output
Pearson:
Age Income
Age 1.000000 -0.001318
Income -0.001318 1.000000- Another method for calculating a correlation coefficient is ‘Spearman’s Rho’.
- The formula looks different but it will give similar results as Pearson’s method.
- In this example we see almost no difference, but this is partly because it is obvious that the Age and Income columns in our dataset have no correlation.
- Formula for Spearmans Rho:
# Using Spearman's rho correlation
print('Spearman: ')
df[['Age', 'Income']].corr(method='spearman')#Output
Spearman:
Age Income
Age 1.000000 -0.001452
Income -0.001452 1.000000Linear regression
Linear Regression can be performed through Ordinary Least Squares (OLS) or Maximum Likelihood Estimation (MLE).
Most Python libraries use OLS to fit linear models.
# Generate data
x = np.random.uniform(low=20, high=260, size=100)
y = 50000 + 2000*x - 4.5 * x**2 + np.random.normal(size=100, loc=0, scale=10000)
# Plot data with Linear Regression
plt.figure(figsize=(16,5))
plt.title('Well fitted but not well fitting: Linear regression plot on quadratic data', fontsize='xx-large')
sns.regplot(x, y)
# Linear regression from scratch
import random
# Create data from regression
xs = np.array(range(1,20))
ys = [0,8,10,8,15,20,26,29,38,35,40,60,50,61,70,75,80,88,96]
# Put data in dictionary
data = dict()
for i in list(xs):
data.update({xs[i-1] : ys[i-1]})
# Slope
m = 0
# y intercept
b = 0
# Learning rate
lr = 0.0001
# Number of epochs
epochs = 100000
# Formula for linear line
def lin(x):
return m * x + b
# Linear regression algorithm
for i in range(epochs):
# Pick a random point and calculate vertical distance and horizontal distance
rand_point = random.choice(list(data.items()))
vert_dist = abs((m * rand_point[0] + b) - rand_point[1])
hor_dist = rand_point[0]
if (m * rand_point[0] + b) - rand_point[1] < 0:
# Adjust line upwards
m += lr * vert_dist * hor_dist
b += lr * vert_dist
else:
# Adjust line downwards
m -= lr * vert_dist * hor_dist
b -= lr * vert_dist
# Plot data points and regression line
plt.figure(figsize=(15,6))
plt.scatter(data.keys(), data.values())
plt.plot(xs, lin(xs))
plt.title('Linear Regression result')
print('Slope: {}\nIntercept: {}'.format(m, b))
[…] Data Science is all about applying maths to data. If you don’t know about maths, you will face difficulties with data science. You can start learning Mathematics for Data Science using Numerical Python and Statistics. […]