# Statistical functions and Statistical Distributions

## Table of contents

- Statistical functions
- Averages and measures of central location
- Measures of Spread
- Statistical distributions
- Discrete distributions
- Continuous distributions
- Law of Large Numbers
- Central Limit Theorem

## Statistical functions

## Averages and measures of central location

## Mean

The mean is the sum of all values divided by the number of values.

import statistics as st

x = [4,5,1,2,7,2,6,9,3,7,2]

print('Mean:', st.mean(x))Mean: 4.363636363636363

##1) ✍ Implement mean from scratch.

def mean(X):

sum = 0

for i, n in enumerate(X):

sum += X[i]

return sum/(i + 1)

print('Mean:', mean(x))Mean: 4.363636363636363x[4, 5, 1, 2, 7, 2, 6, 9, 3, 7, 2]sum = 0sum(x)48

## Median

The median is the middle value in the list of numbers.

print('Median:', st.median(x))Median: 4

##2) ✍ Implement median from scratch.

def median(X):

N = len(X) #0 1 2 3 4 5 6 7

#X.sort()

X = sorted(X)

if N%2 != 0:

return X[N//2]

else:

return round((X[N//2] + X[N//2 - 1])/2, 1)

print('Median:', median(x))Median: 4

## Mode

The mode is the value that occurs most often. If no number in the list is repeated, then there is no mode for the list.

x = [1,2,3,3,5]

print(x)

print('Mode:', st.mode(x))[1, 2, 3, 3, 5]

Mode: 3

## ✍ 3) Implement mode from scratch.

#mode = max(set(x), key=x.count)

def mode(X):

x = [(i, X.count(i)) for i in X]

y = list(set(x))

y.sort(key=lambda y: y[0], reverse=True)

y.sort(key=lambda y: y[1], reverse=True)

return(y[0][0])

x = [1, 2, 2, 3]

print('Mode:', mode(x))Mode: 2X = [2,2,2,3,3,5,5,5,5,8,8]

X.sort()

count = 1

count_number = X[0]

count_max = 1

for i in range(1,len(X)):

if X[i] == X[i-1]:

count = count + 1

if count > count_max:

count_max = count

count_number = X[i]

else:

count = 1

print(count_number)5from collections import Counter

Counter(X).most_common()[0][0]5l = [2,5,4,8,3,8,6,8,6]

s = list(set(l))

c = []

for i in s:

c.append(l.count(i))

d = {k:v for k, v in zip(s,c)}

for k in d.keys():

if(d[k] == max(c)):

print(k)8

## Measures of Spread

## Range

Range is the differene between the largest data value and the smallest data value in the set.

import numpy as np

x = [1, 4, 8]

print(x)

print('Range:', np.ptp(x))[1, 4, 8]

Range: 7

##4) ✍ Implement range from scratch.

np.max(x) - np.min(x)7def my_range(l):

l_min, l_max = float("inf"), -float("inf")

for elem in l:

l_min = min(l_min, elem)

l_max = max(l_max, elem)

return l_max - l_minmy_range([1,2,4,3])3

## Interquartile Range (IQR)

The interquartile range (IQR), also called midspread, is the length of the middle 50% of data interval. The IQR is calculated as the difference between the third or upper quartile ($Q_3$) and the first or lower quartile ($Q_1$).

import numpy as np

from scipy.stats import iqr

#x = [4,5,1,2,7,2,6,9,3,7,2]

x = [1, 4, 5, 6, 6, 7, 11, 12, 15, 17]

print('IQR of data:', iqr(x, ))IQR of data: 6.5q3, q1 = np.percentile(x, [75 , 25], interpolation = "nearest")

print('IQR of data:', q3 - q1)IQR of data: 7print(q1,q3)5 12

##5) ✍ Implement IQR from scratch.

def median(X, N):

X = sorted(X)

if N%2 != 0:

return X[N//2]

else:

return (X[N//2] + X[N//2 - 1])/2.

def Interquartile_Range(X):

N = len(X)

X.sort()

Q1 = median(X[:N//2], N//2)

print(Q1)

if N%2 == 0:

Q3 = median(X[N//2:], N//2)

else:

Q3 = median(X[N//2+1:], N//2)

print(Q3)

return (Q3 - Q1)

x = [1, 4, 5, 6, 6, 7, 11, 12, 15, 17]

print("%0.1f" % Interquartile_Range(x))5

12

7.0

## Variance

Variance is the expectation of the squared deviation of a random variable from its mean. Informally, it measures how far a set of (random) numbers are spread out from their average value.

print(x)

print('Population variance:', st.pvariance(x))

print('Sample variance:', st.variance(x))[1, 4, 5, 6, 6, 7, 11, 12, 15, 17]

Population variance: 23.64

Sample variance: 26.266666666666666

## Standard Deviation

The standard deviation is the square root of its variance.

?np.stdprint('Population standard deviation:', st.pstdev(x))

print('Sample standard deviation:', st.stdev(x))Population standard deviation: 4.862098312457287

Sample standard deviation: 5.125101625008686print('Population standard deviation:',np.std(x, ddof=0)) # ddof: deduct degrees of freedom

print('Sample standard deviation:', np.std(x, ddof=1))Population standard deviation: 4.862098312457287

Sample standard deviation: 5.125101625008686

## 6) ✍ Implement population standard deviation from scratch.

def mean(X, N):

sum = 0

for i in range(len(X)):

sum += X[i]

return sum/N

def standard_deviation(X):

N = len(X)

mu = mean(X, N)

tmp = 0

for i in range(N):

tmp += (X[i] - mu)**2

return (tmp/N)**0.5

x = [1,2,3,4,5,6]

standard_deviation(x)1.707825127659933

## Skewness

Skewness [Ref] is a measure of symmetry, or more precisely, the lack of symmetry. A distribution, or data set, is symmetric if it looks the same to the left and right of the center point.

- The skewness for a normal distribution is zero, and any symmetric data should have a skewness near zero.
- Negative values for the skewness indicate data that are skewed left and positive values for the skewness indicate data that are skewed right. By skewed left, we mean that the left tail is long relative to the right tail.
- Similarly, skewed right means that the right tail is long relative to the left tail.
- If the data are multi-modal, then this may affect the sign of the skewness.

from scipy.stats import skew

import scipy.stats as st

#data = np.random.normal(0, 1, 1000000)

data = np.random.uniform(0, 1, 100000)

print("skewness:",st.skew(data))skewness: -0.0030195334962366182

## Kurtosis

Kurtosis is a measure of whether the data are heavy-tailed or light-tailed relative to a normal distribution. That is, data sets with high kurtosis tend to have heavy tails, or outliers. Data sets with low kurtosis tend to have light tails, or lack of outliers. A uniform distribution would be the extreme case.

from scipy.stats import kurtosis

print("kurtosis:", kurtosis(data))kurtosis: -1.198085296060503

## Statistical distributions

## Discrete distributions

## Discrete uniform distribution

The discrete uniform distribution is a symmetric probability distribution whereby a finite number of values are equally likely to be observed. Its probability mass function is $$ p(x) = \frac{1}{b — a} $$ for $x = a, …, b — 1$.

Python syntax:

scipy.stats.randint(low, high+1)

np.random.randint(low, high+1)

from scipy.stats import randint

low = 3

high = 8

X = randint(low, high + 1) # Declare X to be a uniform random variable

X.rvs(size= 10) # Get 10 random samples form Xarray([8, 6, 7, 6, 6, 5, 4, 3, 5, 7])?sns.distplotimport seaborn as sns

import matplotlib.pyplot as plt

%matplotlib inline

def plot_dist(X, title):

ax = sns.distplot(X, kde=False, color='crimson', hist_kws={"linewidth": 15,'alpha':1})

ax.set(xlabel= 'x', ylabel='Frequency', title = title);

X = randint.rvs(low, high + 1, size= 10000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Uniform Distribution')skewness: -0.019059873221003085

kurtosis: -1.2545266059524007

/usr/local/lib/python3.6/dist-packages/seaborn/distributions.py:2551: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).

warnings.warn(msg, FutureWarning)

## Bernoulli distribution

A Bernoulli trial is an experiment that has two outcomes that can be encoded as success ($x = 1$) or failure ($x = 0$). The Bernoulli Probability Mass Function (PMF) is $$ p(x) = {p^x}{(1 — p)^{(1-x)}} $$ where $x \in {0, 1}$.

Python syntax:

scipy.stats.bernoulli(theta)

np.random.choice([0, 1], p=[1-theta, theta])

from scipy.stats import bernoulli

p = 0.25 # probability of success

X = bernoulli(p) # Declare X to be a bernoulli random variable

X.rvs(10) # Get 10 random samples form Xarray([0, 0, 0, 0, 1, 1, 0, 0, 0, 0])print('P(X = 0) =', X.pmf(0))

print('P(X = 1) =', X.pmf(1))

print('P(X <= 0.5) =', X.cdf(0.5)) # P(X<= 0.5) = P(X = 0) = .75

print('E[X] =', X.mean())

print('Var(X) =', X.var())

print('Std(X) =', X.std())P(X = 0) = 0.75

P(X = 1) = 0.25

P(X <= 0.5) = 0.75

E[X] = 0.25

Var(X) = 0.1875

Std(X) = 0.4330127018922193import seaborn as sns

X = bernoulli.rvs(p= 0.6, size= 1000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Bernoulli Distribution')skewness: -0.37860678521640406

kurtosis: -1.8566569021880999

/usr/local/lib/python3.6/dist-packages/seaborn/distributions.py:2551: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).

warnings.warn(msg, FutureWarning)

## Binomial Distribution

The binomial distribution is used to obtain the probability of observing $x$ successes in $n$ trials, with the probability of success on a single trial denoted by $p$. The binomial distribution assumes that $p$ is fixed for all trials. The formula for the binomial probability mass function is $$ p(x) = {n\choose x} {p^x}{(1 — p)^{(n-x)}} $$ where $x \in {0, 1,…, n}$, and ${n\choose x}$ is the binomial coefficient that calculates different ways of distributing $x$ successes in a sequence of $n$ trials.

Python syntax:

scipy.stats.binom(N, theta)

np.random.binomial(N, theta)

from scipy.stats import binom

n = 7 # Number of trials (number of coins)

p = 0.5 # Probability of success in each trial

sample_size = 5 # Sample size (number of experiments)

X = binom.rvs(n = n, p = p, size = sample_size)

Xarray([6, 4, 2, 3, 2])x = 2 # Number of successes sought

X = binom(n = n, p = p)

print('P(X = %d) = %f' %(x, X.pmf(x)))

print('P(X <= %d) = %f' % (x, X.cdf(x)))

print('E[X] =', X.mean())

print('Var(X) =', X.var())

print('Std(X) =', X.std())X = binom.rvs(n = 10, p = 0.3, size= int(1e6))

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Binomial Distribution')

## Geometric distribution

We perform a series of Bernoulli trials with probability of success $p$ until we get a success. The probability that the $x$-th trial (out of $x$ trials) is the first success is $$ p(x) = {(1 — p)^{(x-1)}}{p} $$ where $x \ge 1$.

Python syntax:

scipy.stats.geom(theta)

np.random.geometric(theta)

`from scipy.stats import geom`

p = 0.2 # probability of success on each trial

X = geom.rvs(p, size= 1000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Geometric Distribution')

## Poisson distribution

The Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time or space if these events occur with a known constant rate and independently of the time since the last event. The Poisson probability that exactly $x$ successes occur in an interval, with the mean number of successes $\lambda$, is $$ p(x) = exp(-\lambda)\frac{{{\lambda}^x}}{{x!}} $$ where $x \ge 1$.

Python syntax:

scipy.stats.poisson(lam)

np.random.poisson(lam)

`from scipy.stats import poisson`

lam = 3 # average number of events per interval

X = poisson.rvs(lam, size= 10000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Poisson Distribution')

## Continuous distributions

## Uniform distribution

Any outcome in a given range has equal probability.

Python syntax:

scipy.stats.uniform(low, high)

np.random.uniform(low, high)

`from scipy.stats import uniform`

low = 3

high = 8

X = uniform.rvs(low, high - low, size= 10000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Uniform Distribution')

## Normal Distribution

The probability density of the normal (Gaussian) distribution is $$ f(x) = \frac{1}{\sqrt{2\pi{{\sigma}²}}} exp(- \frac{(x — \mu)²}{2{\sigma}²}) $$ where $\mu$ is the mean, and ${\sigma}²$ is the variance.

Python syntax:

scipy.stats.norm(mu, sigma)

np.random.normal(mu, sigma)

`from scipy.stats import norm`

print(1e6)

X = norm.rvs(loc=0, scale=1, size=int(1e6))

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Normal Distribution')

## Student’s t distribution

The Student’s t-distribution (or simply the t-distribution) is given by $$ f(x) = \frac{\Gamma(\frac{v+1}{2})}{\sqrt{\pi v} \Gamma(v)} (1 + \frac{x²}{v})^{-\frac{v+1}{2}} $$ where $x$ is a real number, $Î½ > 0$ is the degrees of freedom parameter, and $\Gamma$ is the gamma function.

Python syntax:

scipy.stats.t(nu, mu, sigma)

mu + sigma * np.random.standard_t(nu)

`from scipy.stats import t`

df = 3

X = t.rvs(df, size=1000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 't-distribution')

## Chi-square distribution

The chi-square distribution with $k$ degrees of freedom is the distribution of a sum of the squares of $k$ independent standard normal random variables. Its PDF is $$ f = \frac{1}{2^{\frac{k}{2}} \Gamma(\frac{k}{2})} x^{\frac{k}{2} — 1} exp(-\frac{x}{2}) $$ where $x>0$ and $k>0$.

Python syntax:

scipy.stats.chi2(nu)

np.random.chisquare(nu)

`from scipy.stats import chi2`

df = 55

X = chi2.rvs(df, size=1000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Chi-square Distribution')

## Exponential distribution

The exponential distribution is the probability distribution of the time between events in a Poisson process. The PDF of an exponential distribution is $$ f = \lambda exp(-\lambda x) $$ where $x \ge 0$ and $\lambda > 0$ is the parameter of the distribution, often called the rate parameter.

Python syntax:

scipy.stats.expon(loc=0, scale=1/lambda)

np.random.exponential(1/lambda)

`from scipy.stats import expon`

X = expon.rvs(loc=0, scale=1, size=1000)

print("skewness:",skew(X))

print("kurtosis:", kurtosis(X))

plot_dist(X, 'Exponential Distribution')

## Law of Large Numbers

The law of large numbers (LLN) states that, the average of the results obtained from a large number of trials should be close to the expected value, and will tend to become closer as more trials are performed.

Let ${X_1, X_2, …, X_n}$ be a random sample of size $n$ drawn iid from a distribution with mean $\mu$. According to LLN, we have

$$\bar{X} = \frac{X_1 + X_2 + …+ X_n}{n} \to \mu$$ for large $n$.

For example, a single roll of a fair, six-sided dice produces one of the numbers 1, 2, 3, 4, 5, or 6, each with equal probability. Therefore, the expected value of the average of the rolls is: $$\frac{1+2+3+4+5+6}{6}=3.5$$ According to the law of large numbers, if a large number of six-sided dice are rolled, the average of their values (the sample mean) is likely to be close to 3.5 (the population mean), with the precision increasing as more dice are rolled [Ref]. In the following, we have an illustration of the law of large numbers using a particular run of rolls of a single dice.

np.random.seed(123)

N = 1000 # sample size

rolls = np.random.randint(1, 7, N)

sample_means = []

for i in range(N):

sample_means.append(rolls[:i + 1].mean())

plt.plot(range(1, N+1), 3.5*np.ones(N), 'r', label = 'Population mean')

plt.plot(range(1, N+1), sample_means, 'b', label = 'Sample mean')

plt.xlabel('Sample size (N)')

plt.ylabel("Mean")

plt.legend()

plt.show()#rolls

## 7) ✍ Illustrate LLN for samples drawn iid from a normal distribution with mean of 50 and a variance of 10.

?np.random.normalfrom scipy.stats import norm

import numpy as np

import matplotlib.pyplot as plt

X = norm(loc= 50, scale= np.sqrt(10))

mu = X.mean()

var = X.var()

print('Population mean', mu)

print('Population variance', var)

N = 1000 # sample size

x = X.rvs(size = N)#, random_state= 123)

sample_means = []

for i in range(N):

sample_means.append(sum(x[:i + 1])/(i+1))

plt.plot(range(1, N+1), mu*np.ones(N), 'r', label = 'Population mean')

plt.plot(range(1, N+1), sample_means, 'b', label = 'Sample mean')

plt.xlabel('Sample size (N)')

plt.ylabel("Mean")

plt.legend(loc = 'lower right')

plt.show()Population mean 50.0

Population variance 10.000000000000002

## Central Limit Theorem

The central limit theorem (CLT) dictates that the distribution of sample means (calculated across different random samples from our overall population) will be normally distributed. Mathematically, let ${X_1, X_2, …, X_n}$ be a random sample of size $n$ drawn iid from a distribution with mean $\mu$ and variance ${\sigma}²$, and define the sample average of the random variables as $$ \bar{X} = \frac{X_1 + X_2 + … + X_n}{n} $$ The CLT states that $$ \bar{X} \sim N(\mu, \frac{{\sigma}²}{n}) $$ for large $n$.

## ✍ 8) Test CLT for iid Gamma random variables using scipy.

from scipy.stats import gamma, norm

import numpy as np

import matplotlib.pyplot as plt

a, loc = 2, 3

X = gamma(a, loc)

mu = X.mean()

var = X.var()

print('Population mean:', mu)

print('Population variance:', var)

N = 50 # sample size - corresponding to n

num_trials = 100000 # Number of trials - how many many we want to generate

samples = X.rvs(size = (num_trials, N))

samples_means = samples.mean(axis= 1)

print('Mean of all the sample means:', np.mean(samples_means))

print('Variance of all the sample means =', np.var(samples_means))

print('var/N:', var/N)

plt.hist(samples_means, bins= 100, density=True)

S_N = norm(loc= mu, scale= np.sqrt(var/N))

x = np.linspace(S_N.ppf(0.001), S_N.ppf(0.999), 100)

plt.plot(x, S_N.pdf(x), 'r')

plt.show()Population mean: 5.0

Population variance: 2.0

Mean of all the sample means: 4.999268307684889

Variance of all the sample means = 0.039875471278788675

var/N: 0.04

## 9) ✍ Test CLT for iid uniform random variables in range of [5, 10] using numpy.

import numpy as np

import matplotlib.pyplot as plt

%matplotlib inline

shape, scale = 2., 2.

populaiton_size = 100000

s = np.random.uniform(low= 5, high= 10, size = populaiton_size)

mu = np.mean(s) # mu = (high + low)/2

var = np.var(s) # var = (high - low)**2/12

print('Population mean:', mu)

print('Population variance:', var)

plt.hist(s)

plt.show()

N = 20

num_trials = 10000 # Number of trials

samples = [np.mean(np.random.choice(s, size = N)) for _ in range(num_trials)]

print('Mean of all the sample means:', np.mean(samples))

print('Variance of all the sample means:', np.var(samples))

print('var/N:', var/N)

plt.hist(samples, density=True)

def normal_pdf(x, mu, sigma):

return (1/(sigma*np.sqrt(2*np.pi)) *np.exp( - (x - mu)**2 /(2*sigma**2)))

x = np.linspace(6, 9, 100)

plt.plot(x, normal_pdf(x, mu = mu, sigma = np.sqrt(var/N)), 'r')

plt.show()Population mean: 7.5082694056875

Population variance: 2.0802442142948188Mean of all the sample means: 7.506555371635119

Variance of all the sample means: 0.1031419760228997

var/N: 0.10401221071474094