Exploratory Data Analysis: The Math and Python for Every Step

Exploratory Data Analysis (EDA) is what you do before modeling: understand the shape of the data, find its problems, and form the hypotheses your model will later test. It is not a single step but a sequence of nine, each answering one question:

1. Data Overview — what am I working with?
2. Missing Values — what is absent?
3. Descriptive Statistics — where does the data sit?
4. Univariate Analysis — how is each variable distributed on its own?
5. Bivariate Analysis — how do two variables move together?
6. Correlation Analysis — how strongly are they related?
7. Outlier Detection — what is abnormal?
8. Feature Distribution — is the shape normal, skewed, spread out?
9. Insights & Observations — what patterns, trends, and assumptions emerged?

This article works through all nine. For every step you get three things: what it tells you, the math (with the formula), and the exact Python function to run it. Skim it once end-to-end, then keep it open as a checklist next time you face a new dataset.

All examples assume:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

df = pd.read_csv("your_data.csv")   # one dataframe, n rows, p columns

1. Data Overview

What it tells you: the dimensions of the dataset, what each column is, and what type of data it holds. Before any analysis you need to know you are looking at \(n\) rows and \(p\) columns, and whether each column is numeric (int, float), categorical (object/string), or boolean.

The math: none yet — this step is pure inventory. The only quantities worth computing are the number of rows \(n\), the number of columns \(p\), and the cardinality of each categorical column (number of distinct values), which tells you whether a column is a true category or an ID.

The Python:

df.shape              # (n rows, p columns)
df.columns            # column names
df.dtypes             # data type of each column
df.info()             # shape + dtypes + non-null counts in one call
df.head(10)           # first rows, to see actual values
df.nunique()          # distinct values per column (cardinality)

Reading the output: a column with nunique() equal to n is almost certainly an ID (drop it from modeling); a column with 2 distinct values is binary; mixed-type columns (numbers stored as strings) show up here as object and need casting:

df["sales"] = pd.to_numeric(df["sales"], errors="coerce")
df["date"]  = pd.to_datetime(df["date"], errors="coerce")

2. Missing Values Analysis

What it tells you: which columns have gaps, how many, and at what rate. Missingness drives decisions: drop the column, drop the rows, or impute — and the choice depends on how much is missing.

The math: for column \(j\) with \(n\) rows:

The Python:

df.isna().sum()                    # null count per column
df.isna().mean().round(4) * 100    # missing rate % per column
df.isnull().sum().sort_values(ascending=False).head(20)

# visual: heatmap of missingness (rows x columns)
sns.heatmap(df.isna(), cbar=False, cmap="viridis")
plt.show()

# visual: missingno matrix (install: pip install missingno)
import missingno as msno
msno.matrix(df)

Reading the output: a rule of thumb — under 5% missing is noise you can impute or drop; 5–20% warrants careful imputation; above 20–30% the column is usually more trouble than it is worth, unless it is genuinely informative. Also check where the gaps cluster: missing values that form blocks in the same rows often share a root cause.

3. Descriptive Statistics

What it tells you: the location and spread of each numeric column — where the centre is, how far values typically fall from it.

The math: for a column of values \(x_1, x_2, \dots, x_n\):

• Mean (arithmetic average): $$ \bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i $$
• Median: the middle value after sorting. For \(n\) odd it is \(x_{((n+1)/2)}\); for \(n\) even it is the average of the two middle values. The median is robust to outliers; the mean is not.
• Sample standard deviation (the default std() in pandas): $$ s = \sqrt{\frac{\sum_{i=1}^{n}(x_i - \bar{x})^2}{n - 1}} $$ (the \(n-1\), Bessel’s correction, makes \(s^2\) an unbiased estimate of the population variance; pandas uses this by default, ddof=1)
• Variance: \(s^2\), the average squared deviation from the mean.
• Percentiles: the value below which \(p\%\) of observations fall — e.g. the 25th percentile \(Q_1\) and 75th percentile \(Q_3\).

The Python:

df.describe()                                  # count, mean, std, min, 25%, 50%, 75%, max
df.describe(percentiles=[.01, .05, .95, .99])  # add tails — useful for sales/spend data

df["col"].mean()     df["col"].median()   df["col"].std()
df["col"].var()      df["col"].min()       df["col"].max()
df["col"].quantile([0.25, 0.5, 0.75])

# grouped stats — the most useful descriptive call in business analysis
df.groupby("category")["sales"].agg(["count", "mean", "median", "std", "min", "max"])

Reading the output: if mean ≫ median, the distribution is right-skewed (a few large values pull the mean up) — typical of demand, revenue and lead-time data. If std is several times the mean, the column has very high dispersion relative to its level.

4. Univariate Analysis

What it tells you: the shape of each variable on its own — via histograms (binned frequencies) and box plots (quartile summary).

The math:

• Histogram bins. The number of bins \(k\) controls how much detail you see. Two standard rules:
  • Sturges’ rule (good for roughly normal data): $$ k = \lceil \log_2(n) + 1 \rceil $$
  • Freedman–Diaconis rule (better for skewed data), with bin width: $$ h = \frac{2 \cdot \text{IQR}}{n^{1/3}}, \qquad k = \left\lceil \frac{\max(x) - \min(x)}{h} \right\rceil $$
• Box plot geometry: the box spans \(Q_1\) to \(Q_3\); the line inside is the median \(Q_2\); the whiskers extend to the most extreme points within \(1.5 \times \text{IQR}\) of the box: $$ \text{IQR} = Q_3 - Q_1, \qquad \text{lower fence} = Q_1 - 1.5\,\text{IQR}, \qquad \text{upper fence} = Q_3 + 1.5\,\text{IQR} $$ Points beyond the fences are flagged as potential outliers (step 7).

The Python:

sns.histplot(df["col"], bins="auto", kde=True)   # histogram + density curve
plt.show()

df["col"].plot(kind="hist", bins=30)
plt.show()

sns.boxplot(x=df["col"])                          # horizontal box plot
plt.show()

# one figure per numeric column
df.select_dtypes("number").hist(bins=30, figsize=(12, 10))
plt.tight_layout()
plt.show()

Reading the output: look for the number of peaks (one = unimodal, two = bimodal, often a sign of two underlying groups), the tails (long right tail = skewed), and the gap between box-plot fences and the extreme points.

5. Bivariate Analysis

What it tells you: how two variables behave together — direction of relationship, linearity, and clusters. Scatter plots for numeric pairs, pair plots for the whole matrix.

The math: the strength of the linear relationship is captured by the sample covariance and the Pearson correlation coefficient \(r\):

Covariance is in the original units (hard to interpret); \(r\) is unitless and bounded in \([-1, 1]\). The sign tells direction, the magnitude tells strength — but scatter plots are still essential, because \(r\) says nothing about shape (two clusters, a curve, or one extreme point can all produce the same \(r\)).

The Python:

sns.scatterplot(data=df, x="feature_a", y="target")
plt.show()

sns.pairplot(df, diag_kind="kde")          # scatter matrix for all numeric columns
plt.show()

pd.plotting.scatter_matrix(df, figsize=(12, 12), diagonal="hist")
plt.show()

df[["feature_a", "feature_b"]].cov()       # covariance matrix
df[["feature_a", "feature_b"]].corr()      # Pearson correlation matrix

Reading the output: an upward-sloping cloud is positive correlation, downward is negative. Curved patterns mean a non-linear relationship — Pearson \(r\) will understate it; consider a transformation (log, square root) or rank-based measures.

6. Correlation Analysis

What it tells you: the full pairwise relationship matrix between all numeric columns, usually rendered as a heatmap. This is the step that decides which features to trust and which are redundant.

The math: the correlation matrix \(R\) is the \(p \times p\) matrix whose entry \(r_{jk}\) is the Pearson correlation between column \(j\) and column \(k\) (as in step 5). The diagonal is always 1. When the relationship is monotonic but not linear, the Spearman rank correlation is preferred — Pearson computed on the ranks of the data:

The Python:

corr = df.corr()                                  # Pearson
corr = df.corr(method="spearman")                 # Spearman (rank-based, robust to outliers)

# heatmap with values
sns.heatmap(corr, annot=True, fmt=".2f", cmap="RdBu_r",
            vmin=-1, vmax=1, square=True)
plt.show()

# mask the redundant upper triangle
mask = np.triu(np.ones_like(corr, dtype=bool), k=1)
sns.heatmap(corr, mask=mask, annot=True, fmt=".2f", cmap="RdBu_r",
            vmin=-1, vmax=1, square=True)
plt.show()

# flag the strongest pairs
corr_unstacked = corr.where(~np.eye(corr.shape[0], dtype=bool)).stack().reset_index()
corr_unstacked.columns = ["col_a", "col_b", "r"]
corr_unstacked["abs_r"] = corr_unstacked["r"].abs()
corr_unstacked.sort_values("abs_r", ascending=False).head(10)

Reading the output: |r| < 0.3 is weak, 0.3–0.7 moderate, > 0.7 strong. Two features with |r| > 0.8–0.9 are near-duplicates — keep one and drop the other before feeding a model (multicollinearity inflates regression coefficients and confuses tree-based feature importance). A column with near-zero correlation to everything may be useless noise — or may just need a transformation.

7. Outlier Detection

What it tells you: which values are so far from the bulk of the data that they deserve investigation — data-entry errors, genuine extremes, or spikes that matter (a stockout, a promotion, a one-off order).

The math: two standard rules.

• IQR (Tukey) rule — robust, distribution-free. Flag \(x\) as an outlier if:

  $$ x < Q_1 - 1.5 \cdot \text{IQR} \quad \text{or} \quad x > Q_3 + 1.5 \cdot \text{IQR} $$

  For extreme outliers, use \(3 \cdot \text{IQR}\) instead of \(1.5\).

• z-score rule — assumes approximate normality. Flag \(x\) if:

  $$ |z| = \left|\frac{x - \bar{x}}{s}\right| > 3 $$

  (under the normal distribution, only ~0.3% of values fall beyond \(3\sigma\)).

• Modified z-score (robust) — uses the median and MAD instead of mean and std, so outliers don’t distort the very statistic that detects them:

  $$ M_i = \frac{0.6745 \cdot (x_i - \text{median}(x))}{\text{MAD}}, \qquad \text{flag if } |M_i| > 3.5 $$

  where MAD is the median absolute deviation: \(\text{MAD} = \text{median}(|x_i - \text{median}(x)|)\).

The Python:

# IQR method
Q1, Q3 = df["col"].quantile([0.25, 0.75])
iqr = Q3 - Q1
lower, upper = Q1 - 1.5 * iqr, Q3 + 1.5 * iqr
outliers = df[(df["col"] < lower) | (df["col"] > upper)]

# z-score method
z = np.abs(stats.zscore(df["col"]))            # scipy: sample std, ddof=0
outliers = df[z > 3]

# modified z-score (robust)
mad = stats.median_abs_deviation(df["col"])
mod_z = 0.6745 * (df["col"] - df["col"].median()) / mad
outliers = df[mod_z.abs() > 3.5]

# visual
sns.boxplot(x=df["col"])                       # dots beyond whiskers = flagged
plt.show()

Reading the output: never delete outliers silently. Classify them first — data error (fix or drop), genuine extreme (keep, or cap/winsorize), or signal (a spike your model should learn). For tree models outliers barely matter; for linear models and any mean-based statistic they can dominate.

8. Feature Distribution

What it tells you: the shape of each variable — symmetric or skewed, flat or peaked, normal or heavy-tailed. Distribution shape decides which statistics are trustworthy and which transformations to apply.

The math:

• Skewness (sample, adjusted for bias):

  $$ g_1 = \frac{n}{(n-1)(n-2)} \sum_{i=1}^{n} \left( \frac{x_i - \bar{x}}{s} \right)^3 $$

  \(g_1 > 0\): right-skewed (long tail to the right — the typical sales/demand shape); \(g_1 < 0\): left-skewed; \(|g_1| > 1\) is highly skewed.

• Excess kurtosis (sample):

  $$ g_2 = \frac{n(n+1)}{(n-1)(n-2)(n-3)} \sum_{i=1}^{n} \left( \frac{x_i - \bar{x}}{s} \right)^4 - \frac{3(n-1)^2}{(n-2)(n-3)} $$

  A normal distribution has excess kurtosis 0. \(> 0\): heavy tails (more extremes than normal); \(< 0\): light tails (flatter, wider spread).

• Empirical rule (if approximately normal): ~68% of values within \(\bar{x} \pm 1s\), ~95% within \(\pm 2s\), ~99.7% within \(\pm 3s\).

The Python:

df.skew(numeric_only=True).sort_values(ascending=False)   # skewness per column
df.kurt(numeric_only=True).sort_values(ascending=False)   # excess kurtosis per column

stats.skew(df["col"], bias=False)          # scipy, bias-corrected
stats.kurtosis(df["col"], bias=False)      # scipy, excess kurtosis

sns.histplot(df["col"], kde=True)          # density overlay shows shape
plt.show()

# normality check (interpret: p > 0.05 is consistent with normality)
stats.shapiro(df["col"].dropna().sample(min(5000, len(df)), random_state=42))
stats.normaltest(df["col"].dropna())

Reading the output: strongly skewed columns are the norm in supply-chain data (sales, lead times, costs). Common fixes before modeling: log transform \(\log(x + 1)\), square root, or Box–Cox (\(x^{(\lambda)} = (x^\lambda - 1)/\lambda\) for \(\lambda \neq 0\), \(\log x\) for \(\lambda = 0\)). Box–Cox needs strictly positive values.

from scipy.stats import boxcox
df["col_boxcox"], fitted_lambda = boxcox(df["col"] + 1e-9)   # works on positive data

9. Insights & Observations

What it tells you: the payoff step — patterns, trends, and the assumptions you will carry into modeling. EDA ends not with plots but with a written list of findings: which groups differ, which relationships are strong, which columns are problematic.

The math: mostly none — this step is about comparison. The one calculation that consistently produces insight is a grouped comparison: mean, median, or total of a target across categories, plus a quick significance test to check whether the difference is real or chance:

• Two-group test (t-test): $$ t = \frac{\bar{x}_A - \bar{x}_B}{\sqrt{\frac{s_A^2}{n_A} + \frac{s_B^2}{n_B}}} $$
• Categorical association (chi-squared): \(\chi^2 = \sum \frac{(O - E)^2}{E}\) over the contingency table cells.

The Python:

# patterns: groupwise summary
df.groupby("category")["target"].agg(["count", "mean", "median", "sum"]).sort_values("sum", ascending=False)

# trends: target over time
df.groupby("month")["target"].sum().plot(kind="line")
plt.show()

# two-group difference + test
group_a = df.loc[df["flag"] == "A", "target"]
group_b = df.loc[df["flag"] == "B", "target"]
stats.ttest_ind(group_a, group_b, equal_var=False)     # Welch's t-test

# categorical association
table = pd.crosstab(df["category"], df["segment"])
stats.chi2_contingency(table)

# pivot table — the workhorse of business EDA
pd.pivot_table(df, values="target", index="category", columns="segment", aggfunc="median")

Reading the output: write down 5–10 concrete findings before modeling: “sales are right-skewed; category X dominates 40% of volume; lead time correlates 0.6 with order size; 8% of rows are missing warehouse data and they cluster in region Y.” These observations become the feature-engineering and assumption section of your modeling notes — and the sanity checks you use to judge whether a model’s output makes sense.

The complete workflow, end to end

A compact template that runs every step in order:

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
from scipy import stats

df = pd.read_csv("your_data.csv")

# 1. Overview
print(df.shape, df.dtypes.value_counts(), df.nunique(), sep="\n")

# 2. Missing
missing = df.isna().mean().round(4) * 100
print(missing[missing > 0].sort_values(ascending=False))
sns.heatmap(df.isna(), cbar=False)

# 3. Descriptive
print(df.describe(percentiles=[.01, .05, .95, .99]).T)

# 4. Univariate
df.select_dtypes("number").hist(bins=30, figsize=(12, 10))

# 5. Bivariate
sns.pairplot(df.select_dtypes("number").sample(min(2000, len(df)), random_state=1))

# 6. Correlation
corr = df.corr(numeric_only=True)
sns.heatmap(corr, annot=True, fmt=".2f", cmap="RdBu_r", vmin=-1, vmax=1, square=True)

# 7. Outliers (IQR per numeric column)
num = df.select_dtypes("number")
q1, q3 = num.quantile([.25, .75]); iqr = q3 - q1
outlier_flags = (num < q1 - 1.5*iqr) | (num > q3 + 1.5*iqr)
print("outlier rows:", outlier_flags.any(axis=1).sum(), "of", len(df))

# 8. Distribution
print(df.skew(numeric_only=True).sort_values(ascending=False))

# 9. Insights — your turn: write down what you found.

A quick reference table

The order matters: overview → missing → descriptive gives you the lay of the land; univariate → bivariate → correlation builds the picture variable by variable; outliers → distribution tells you what to fix and how; insights closes the loop. Ten minutes of structured EDA saves hours of debugging a model built on data you never looked at.