Distributions and Descriptive Statistics in `Python`

This tutorial introduces distributions and descriptive statistics in Python using pandas and helper functions that mirror R's syntax.

Getting Started

Install and Import Packages

%pip install pandas plotnine scipy

import os, sys
import pandas as p
from plotnine import *
sys.path.append(os.path.abspath('functions'))
from functions_distributions import *

Our Data

sw = p.Series([4.5, 5, 5.5, 5, 5.5, 6.5, 6.5, 6, 5, 4])
sw

Size

Length

len(sw)

Location

Mean and Median

sw.mean()

sw.median()

Mode

sw.mode()

Spread (1)

Percentiles

sw.quantile(q=0)  # min
sw.quantile(q=1)  # max
sw.quantile(q=.75)

Spread (2)

Standard Deviation, Variance, CV, SE

# Manual SD (sample)
x = ((sw - sw.mean())**2).sum()
x = x / (len(sw) - 1)
x**0.5

sw.std()
sw.var()
sw.std()**2
sw.std() / sw.mean()  # CV

se = sw.std() / (len(sw)**0.5)
se

Shape

Skewness and Kurtosis

diff = sw - sw.mean()
n = len(sw) - 1
sigma = sw.std()
sum(diff**3) / (n * sigma**3)
sum(diff**4) / (n * sigma**4)

# Using helper functions mirroring R
skewness(sw)
kurtosis(sw)

Finding Parameters for Your Distributions

sw = p.Series([4.5, 5, 5.5, 5, 5.5, 6.5, 6.5, 6, 5, 4])
mymean = sw.mean()
mysd = sw.std()

Common Distributions

Normal

mynorm = rnorm(n=1000, mean=mymean, sd=mysd)
h1 = hist(mynorm)
h1.save("plotnine_figures/02_hist_normal.png", dpi=100, width=6, height=4)

Poisson

mypois = rpois(n=1000, mu=mymean)
h2 = hist(mypois)
h2.save("plotnine_figures/02_hist_poisson.png", dpi=100, width=6, height=4)

Exponential

myrate_e = 1 / sw.mean()
myexp = rexp(n=1000, rate=myrate_e)
h3 = hist(myexp)
h3.save("plotnine_figures/02_hist_exponential.png", dpi=100, width=6, height=4)

Gamma

myshape = sw.mean()**2 / sw.var()
myrate = 1 / (sw.var() / sw.mean())
mygamma = rgamma(n=1000, shape=myshape, rate=myrate)
h4 = hist(mygamma)
h4.save("plotnine_figures/02_hist_gamma.png", dpi=100, width=6, height=4)

Weibull

from scipy import stats as fitdistr
myshape_w, loc, myscale_w = fitdistr.weibull_min.fit(sw, floc=0)
myweibull = rweibull(n=1000, shape=myshape_w, scale=myscale_w)
h5 = hist(myweibull)
h5.save("plotnine_figures/02_hist_weibull.png", dpi=100, width=6, height=4)

Comparing Distributions

mysim = p.concat([
  p.DataFrame({'x': sw, 'type': "Observed"}),
  p.DataFrame({'x': mynorm, 'type': "Normal"}),
  p.DataFrame({'x': mypois, 'type': "Poisson"}),
  p.DataFrame({'x': mygamma, 'type': "Gamma"}),
  p.DataFrame({'x': myexp, 'type': "Exponential"}),
  p.DataFrame({'x': myweibull, 'type': "Weibull"})
])

g1 = (ggplot(mysim, aes(x='x', fill='type')) +
  geom_density(alpha=0.5) +
  labs(x='Seawall Height (m)', y='Density (Frequency)', subtitle='Which distribution fits best?', fill='Type'))
g1.save("plotnine_figures/02_density_comparison.png", dpi=100, width=8, height=6)

g2 = g1 + xlim(0,10)
g2.save("plotnine_figures/02_density_comparison_xlim.png", dpi=100, width=8, height=6)

Distribution comparison with x-axis limits

LC 01

Simulating a normal distribution from sw

Simulate 1000 draws from a normal distribution using your sw mean and standard deviation. What are the simulated mean and sd? How close are they to sw's?

Use rnorm() fed with sw.mean() and sw.std(), then compare the simulated Series' own .mean() and .std() back to sw's.

mymean = sw.mean(); mysd = sw.std()
m = rnorm(1000, mean=mymean, sd=mysd)
[m.mean(), m.std()]

Conclusion

You computed size, location, spread, and shape statistics and compared common simulated distributions using helper functions that mirror R.

Distributions and Descriptive Statistics in Python