Probabilistic Methods — Interactive Solver

List 7 — Data Analysis

Descriptive statistics, distributions in data, correlation

Biased Die: Empirical vs Theoretical Distribution

A six-sided die is suspected to be biased. Its true probabilities are $p=[0.14, 0.16, 0.17, 0.18, 0.16, 0.19]$ . Roll it $n=600$ times, record results, compute empirical relative frequencies, and compare them to the theoretical probabilities. Discuss the chi-square intuition.

a)Generate a dataset of 600 rolls using the given probability vector.
b)Compute the empirical relative frequency of each face.
c)Compare empirical frequencies to theoretical probabilities in a table.
d)Compute the chi-square statistic $\chi^2 = \sum_{k=1}^{6} \frac{(O_k - E_k)^2}{E_k}$ and interpret it.
e)Use the live simulator below to watch convergence as $n$ grows.

Step 1 — Generate Dataset

Use NumPy to simulate 600 rolls of a biased die. Each row stores a single roll outcome (1–6).

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(42)
probs = [0.14, 0.16, 0.17, 0.18, 0.16, 0.19]
faces = [1, 2, 3, 4, 5, 6]
rolls = np.random.choice(faces, size=600, p=probs)
df = pd.DataFrame({'roll': rolls})

counts = df['roll'].value_counts().sort_index()
df_summary = pd.DataFrame({
    'face': faces,
    'observed': counts.values,
    'empirical_freq': (counts.values / 600).round(4),
    'theoretical_prob': probs
})
print(df_summary.to_string(index=False))

Step 2 — Empirical Relative Frequencies

For each face $k$ , the empirical relative frequency is $\hat{p}_k = O_k / n$ . With $n=600$ rolls the law of large numbers guarantees $\hat{p}_k \to p_k$ .

\hat{p}_k = \frac{O_k}{n}, \quad k = 1,\ldots,6

Step 3 — Chi-Square Calculation

The expected count for face $k$ is $E_k = n \cdot p_k$ . Compute the chi-square statistic and compare to the critical value $\chi^2_{0.05,5} \approx 11.07$ .

\chi^2 = \sum_{k=1}^{6} \frac{(O_k - E_k)^2}{E_k}

dataset generator · Pythonread-only

import numpy as np
from scipy.stats import chisquare

probs = [0.14, 0.16, 0.17, 0.18, 0.16, 0.19]
n = 600
expected = [p * n for p in probs]

np.random.seed(42)
faces = [1, 2, 3, 4, 5, 6]
rolls = np.random.choice(faces, size=n, p=probs)
observed = [np.sum(rolls == k) for k in faces]

chi2, p_val = chisquare(observed, f_exp=expected)
print(f"Observed: {observed}")
print(f"Expected: {expected}")
print(f"Chi-square: {chi2:.4f}, p-value: {p_val:.4f}")
print("Conclusion: small chi2 => empirical matches theoretical.")

Step 4 — Live Simulator

Use the die simulator below. Increase the number of rolls and observe how the bar of each face converges toward its theoretical probability. This illustrates the Law of Large Numbers.

Die rolls — empirical vs theoretical

rolls: 200

Answer

With

n=600

rolls the empirical frequencies should lie close to the theoretical probabilities. The chi-square statistic will typically be well below the critical value of 11.07, indicating no statistically significant deviation from the declared bias.

Explanation tip: Face 6 has the highest probability (0.19); face 1 has the lowest (0.14). With small

n

individual runs can look quite different — the live simulator makes convergence visible.

Biased Coin: Law of Large Numbers in Action

A coin has $P(\text{heads}) = 0.52$ . Flip it $n=1000$ times and track the running relative frequency of heads. Observe how the empirical proportion converges to the true probability as $n$ increases.

a)Simulate 1000 coin flips with $P(H)=0.52$ .
b)Compute the running (cumulative) relative frequency after each flip.
c)At what $n$ does the running frequency stay within $\pm 0.02$ of 0.52?
d)Use the live LLN simulator below to explore convergence interactively.

Step 1 — Simulate Flips

Encode heads as 1, tails as 0. The running mean after $k$ flips is the empirical estimate of $P(H)$ .

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(7)
flips = np.random.binomial(1, 0.52, size=1000)
running_freq = np.cumsum(flips) / np.arange(1, 1001)

df = pd.DataFrame({
    'flip_number': range(1, 1001),
    'outcome': flips,
    'running_freq': running_freq.round(4)
})
print(df.head(10).to_string(index=False))
print("...")
print(df.tail(5).to_string(index=False))

Step 2 — Law of Large Numbers Formula

By the (weak) Law of Large Numbers, for any $\varepsilon>0$ : as $n\to\infty$ the probability that $|\hat{p}_n - p|>\varepsilon$ goes to 0.

\hat{p}_n = \frac{1}{n}\sum_{i=1}^{n} X_i \xrightarrow{\,p\,} p = 0.52

Step 3 — First Stable Crossing

Find the smallest $n$ after which the running frequency stays within $\pm 0.02$ of the true probability for the remaining flips.

dataset generator · Pythonread-only

import numpy as np

np.random.seed(7)
flips = np.random.binomial(1, 0.52, size=1000)
running_freq = np.cumsum(flips) / np.arange(1, 1001)

eps = 0.02
within = np.abs(running_freq - 0.52) <= eps
stable_n = None
for i in range(len(within)):
    if np.all(within[i:]):
        stable_n = i + 1
        break
print(f"Stable within +-{eps} of 0.52 from flip n={stable_n}")
print(f"Final running frequency: {running_freq[-1]:.4f}")

Step 4 — Live LLN Simulator

Use the simulator below to flip the coin thousands of times interactively. Watch how the running relative frequency approaches the red line at $p=0.52$ .

Law of Large Numbers — coin

flips: 200

Relative frequency of heads converges to p = 0.52 as the number of flips grows.

Answer

The running relative frequency converges to 0.52 as

n

grows. By around

n=300

–

500

the estimate typically stabilizes within

\pm 0.02

. The key takeaway is that randomness averages out over many trials.

Explanation tip: The spread

P(H)=0.52

is deliberately close to 0.50 to illustrate that detecting a small bias requires a large sample.

Exam Scores: Comparing Two Groups

Two teaching groups took the same exam. Group A scores follow $N(68, 12^2)$ and Group B scores follow $N(74, 18^2)$ , both with $n=200$ students. Compare the distributions on mean, median, standard deviation, and pass rate ( $\text{score} \geq 60$ ).

a)Generate exam scores for both groups.
b)Compute mean, median, and standard deviation for each group.
c)Compute the pass rate ( $\text{score}\geq 60$ ) for each group.
d)Explain why Group B has a higher mean yet a lower pass rate is possible in some simulations.

Step 1 — Generate Scores

Simulate exam scores. Scores are clipped to [0, 100] since the exam is out of 100 marks.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(21)
scores_a = np.clip(np.random.normal(68, 12, 200), 0, 100).round(1)
scores_b = np.clip(np.random.normal(74, 18, 200), 0, 100).round(1)

df = pd.DataFrame({
    'student': list(range(1, 201)) * 2,
    'group': ['A'] * 200 + ['B'] * 200,
    'score': np.concatenate([scores_a, scores_b])
})
print(df.head(5).to_string(index=False))
print(df.groupby('group')['score'].describe().round(2))

Step 2 — Descriptive Statistics

Compute the core descriptive statistics. Sample mean and sample standard deviation use the unbiased estimators.

\bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i, \quad s = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2}

dataset generator · Pythonread-only

import numpy as np

np.random.seed(21)
scores_a = np.clip(np.random.normal(68, 12, 200), 0, 100)
scores_b = np.clip(np.random.normal(74, 18, 200), 0, 100)

for label, scores in [('A', scores_a), ('B', scores_b)]:
    print(f"Group {label}: mean={scores.mean():.2f}, "
          f"median={np.median(scores):.2f}, "
          f"std={scores.std(ddof=1):.2f}, "
          f"pass_rate={(scores >= 60).mean()*100:.1f}%")

Step 3 — Interpretation

Group B has a higher theoretical mean (74 vs 68) but also much greater spread ( $\sigma=18$ vs $\sigma=12$ ). More Group B students score both very high and very low. The heavier lower tail can erode Group B's pass rate despite its higher mean.

P(X_B \geq 60) = P\!\left(Z \geq \frac{60-74}{18}\right) = P(Z \geq -0.78) \approx 0.782

Answer

Group A: mean

\approx 68

, std

\approx 12

, pass rate

\approx 84\%

. Group B: mean

\approx 74

, std

\approx 18

, pass rate

\approx 78\%

. Higher mean does not guarantee higher pass rate when variance differs substantially.

Explanation tip: Always report variability alongside the mean. Two distributions with the same mean but different spread have very different practical implications.

Delivery Times: Outliers and Robust Statistics

Delivery times (in minutes) follow a lognormal distribution with $\mu=3.5, \sigma=0.4$ . Five outliers representing "lost packages" (times 300–500 min) are injected. Detect outliers using the $1.5 \cdot \text{IQR}$ rule and compare mean/median and std/IQR before and after removal.

a)Generate 200 delivery times from the lognormal distribution and inject 5 outliers.
b)Compute Q1, Q3, IQR, and the outlier fences.
c)Identify outliers using the $1.5 \cdot \text{IQR}$ rule.
d)Compare mean vs median and std vs IQR with and without outliers.

Step 1 — Generate Data with Outliers

The lognormal distribution is right-skewed — a natural model for delivery times. We inject five extreme values to represent failed deliveries.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(33)
normal_times = np.random.lognormal(mean=3.5, sigma=0.4, size=200)
outlier_times = np.random.uniform(300, 500, size=5)
times = np.concatenate([normal_times, outlier_times])
np.random.shuffle(times)

df = pd.DataFrame({'delivery_time': times.round(2)})
print(df['delivery_time'].describe().round(2))
print(f"Total records: {len(df)}")

Step 2 — IQR Outlier Detection

The IQR rule flags any value below $Q_1 - 1.5 \cdot \text{IQR}$ or above $Q_3 + 1.5 \cdot \text{IQR}$ as a potential outlier.

\text{IQR} = Q_3 - Q_1, \quad \text{lower} = Q_1 - 1.5\,\text{IQR}, \quad \text{upper} = Q_3 + 1.5\,\text{IQR}

dataset generator · Pythonread-only

import numpy as np

np.random.seed(33)
times = np.concatenate([
    np.random.lognormal(3.5, 0.4, 200),
    np.random.uniform(300, 500, 5)
])

q1, q3 = np.percentile(times, [25, 75])
iqr = q3 - q1
lower, upper = q1 - 1.5 * iqr, q3 + 1.5 * iqr
outliers = times[(times < lower) | (times > upper)]
print(f"Q1={q1:.2f}, Q3={q3:.2f}, IQR={iqr:.2f}")
print(f"Fences: [{lower:.2f}, {upper:.2f}]")
print(f"Outliers detected: {len(outliers)}")
print(f"Outlier values: {sorted(outliers.round(1))}")

Step 3 — Robust vs Non-Robust Statistics

Outliers inflate the mean and standard deviation dramatically. The median and IQR are resistant to extreme values.

dataset generator · Pythonread-only

import numpy as np

np.random.seed(33)
times = np.concatenate([
    np.random.lognormal(3.5, 0.4, 200),
    np.random.uniform(300, 500, 5)
])
q1, q3 = np.percentile(times, [25, 75])
iqr = q3 - q1
clean = times[(times >= q1 - 1.5*iqr) & (times <= q3 + 1.5*iqr)]

print("With outliers:   mean={:.1f}, median={:.1f}, std={:.1f}, IQR={:.1f}".format(
    times.mean(), np.median(times), times.std(ddof=1), iqr))
print("Without outliers: mean={:.1f}, median={:.1f}, std={:.1f}, IQR={:.1f}".format(
    clean.mean(), np.median(clean), clean.std(ddof=1),
    np.percentile(clean,75)-np.percentile(clean,25)))

Answer

The 5 injected outliers raise the mean by roughly 30–50 minutes but barely shift the median. IQR is similarly stable. Use median and IQR for skewed or outlier-prone data; reserve mean and std for roughly symmetric distributions.

Customer Survey: Conditional Frequencies and Cross-Tabs

A customer survey records: age group (18–30, 31–50, 51+), acquisition channel (Online, Store, Referral), satisfaction score (1–5), and renewal flag (Yes/No). Analyze renewal rates by channel and satisfaction level.

a)Generate 500 synthetic survey records.
b)Build a cross-tabulation of channel vs renewal.
c)Compute conditional renewal rate by satisfaction score.
d)Identify which channel and satisfaction combination has the highest renewal rate.

Step 1 — Generate Survey Data

Higher satisfaction scores are linked to higher renewal probability. Referral customers have a slightly better renewal baseline.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(55)
n = 500
channels = np.random.choice(['Online','Store','Referral'], n, p=[0.5,0.3,0.2])
ages = np.random.choice(['18-30','31-50','51+'], n, p=[0.35,0.45,0.20])
satisfaction = np.random.choice([1,2,3,4,5], n, p=[0.08,0.12,0.25,0.35,0.20])
renew_prob = 0.3 + 0.12 * (satisfaction - 1)
renew_prob += np.where(channels == 'Referral', 0.05, 0)
renewal = (np.random.rand(n) < renew_prob).astype(int)

df = pd.DataFrame({'age':ages,'channel':channels,
                   'satisfaction':satisfaction,'renewal':renewal})
print(df.head(8).to_string(index=False))

Step 2 — Cross-Tabulation: Channel vs Renewal

The cross-tab shows absolute counts. Dividing each row by its total gives the conditional renewal rate per channel.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(55)
n = 500
channels = np.random.choice(['Online','Store','Referral'], n, p=[0.5,0.3,0.2])
satisfaction = np.random.choice([1,2,3,4,5], n, p=[0.08,0.12,0.25,0.35,0.20])
renew_prob = 0.3 + 0.12*(satisfaction-1)
renew_prob += np.where(channels=='Referral', 0.05, 0)
renewal = (np.random.rand(n) < renew_prob).astype(int)
df = pd.DataFrame({'channel':channels,'satisfaction':satisfaction,'renewal':renewal})

crosstab = pd.crosstab(df['channel'], df['renewal'],
                        margins=True, margins_name='Total')
print(crosstab)
rate = df.groupby('channel')['renewal'].mean().round(3)
print("
Renewal rate by channel:")
print(rate)

Step 3 — Renewal Rate by Satisfaction Score

Conditional renewal rate: $P(\text{renewal} \mid \text{satisfaction}=k)$ . This reveals the monotone relationship between satisfaction and loyalty.

P(\text{renewal} \mid S=k) = \frac{|\{i : \text{renewal}_i=1,\, S_i=k\}|}{|\{i : S_i=k\}|}

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(55)
n = 500
channels = np.random.choice(['Online','Store','Referral'], n, p=[0.5,0.3,0.2])
satisfaction = np.random.choice([1,2,3,4,5], n, p=[0.08,0.12,0.25,0.35,0.20])
renew_prob = 0.3 + 0.12*(satisfaction-1)
renew_prob += np.where(channels=='Referral', 0.05, 0)
renewal = (np.random.rand(n) < renew_prob).astype(int)
df = pd.DataFrame({'channel':channels,'satisfaction':satisfaction,'renewal':renewal})

rate_s = df.groupby('satisfaction')['renewal'].mean().round(3)
print("Renewal rate by satisfaction:")
print(rate_s)
best = df.groupby(['channel','satisfaction'])['renewal'].mean().idxmax()
print(f"Highest renewal: channel={best[0]}, satisfaction={best[1]}")

Answer

Referral customers with satisfaction

\geq 4

show the highest renewal rates (

\approx 75\%

). Renewal rate increases monotonically with satisfaction score. Online channel has the most customers but not the best loyalty.

Explanation tip: Conditional frequencies reveal relationships hidden in marginal totals. Always segment by relevant covariates before reporting summary statistics.

Factory Measurements: Machine Comparison and Spec Limits

Three factory machines (M1, M2, M3) produce parts with a target diameter of $50\,\text{mm}$ . Specification limits are $[49, 51]\,\text{mm}$ . Each machine produces 150 parts. Compare mean offset from target, standard deviation, and fraction within spec.

a)Generate measurements: M1 ~ N(50.1, 0.3²), M2 ~ N(49.8, 0.6²), M3 ~ N(50.0, 0.2²).
b)Compute mean, std, and bias (mean − 50) for each machine.
c)Compute the fraction of parts within spec [49, 51].
d)Rank the machines by quality; explain the bias-variance trade-off.

Step 1 — Generate Measurement Data

M1 is slightly high-biased but precise. M2 is low-biased and imprecise. M3 is centred and most precise.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(99)
m1 = np.random.normal(50.1, 0.3, 150)
m2 = np.random.normal(49.8, 0.6, 150)
m3 = np.random.normal(50.0, 0.2, 150)

df = pd.DataFrame({
    'measurement': np.concatenate([m1, m2, m3]),
    'machine': ['M1']*150 + ['M2']*150 + ['M3']*150
})
print(df.groupby('machine')['measurement'].describe().round(4))

Step 2 — Bias and Variability

Bias measures systematic error (offset from target). Variance measures random error (scatter around the machine mean).

\text{bias}_j = \bar{x}_j - 50, \quad s_j = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_{ij}-\bar{x}_j)^2}

dataset generator · Pythonread-only

import numpy as np

np.random.seed(99)
machines = {
    'M1': np.random.normal(50.1, 0.3, 150),
    'M2': np.random.normal(49.8, 0.6, 150),
    'M3': np.random.normal(50.0, 0.2, 150),
}
for name, data in machines.items():
    bias = data.mean() - 50
    std = data.std(ddof=1)
    in_spec = ((data >= 49) & (data <= 51)).mean()
    print(f"{name}: bias={bias:+.4f} mm, std={std:.4f} mm, "
          f"within_spec={in_spec*100:.1f}%")

Step 3 — Ranking and Interpretation

M3 dominates: it is centred on target and has the smallest spread, yielding the highest fraction within spec. M2 is worst due to combined bias and high variance. M1 has small bias but is less precise than M3.

Answer

M3 (mean

\approx 50.00

, std

\approx 0.20

\approx 99.7\%

within spec) is best. M1 (std

\approx 0.30

, slight positive bias) is second. M2 (std

\approx 0.60

, negative bias) fails spec most often.

Explanation tip: In manufacturing, both bias and variance matter. A centred machine with low spread minimises rejects. The process capability index

C_{pk}

formally combines both.

Call-Center Requests: Poisson Distribution by Hour

A call-center logs hourly request counts over 60 days. The arrival rate $\lambda$ varies by time of day: $\lambda=8$ (off-peak), $\lambda=25$ (peak 10–12, 14–16), $\lambda=15$ (shoulder 8–10, 12–14, 16–18). Verify that empirical mean $\approx$ variance (Poisson property).

a)Generate 60 days × 10 hours of Poisson-distributed counts.
b)Aggregate to get the empirical mean and variance per hour.
c)Verify mean ≈ variance for each hour band (Poisson property).
d)Identify the peak hour and compare its distribution to the off-peak hour.

Step 1 — Generate Hourly Request Counts

Each cell is an independent Poisson draw. The resulting dataset has 600 rows (day × hour).

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(77)
hours = list(range(8, 18))
lambda_map = {8:15, 9:15, 10:25, 11:25, 12:15,
              13:15, 14:25, 15:25, 16:15, 17:8}
records = []
for day in range(1, 61):
    for hour in hours:
        lam = lambda_map[hour]
        count = np.random.poisson(lam)
        records.append({'day': day, 'hour': hour,
                        'lambda': lam, 'requests': count})
df = pd.DataFrame(records)
print(df.head(10).to_string(index=False))

Step 2 — Verify Poisson Property (Mean ≈ Variance)

For a Poisson random variable $X \sim \text{Pois}(\lambda)$ , both the mean and variance equal $\lambda$ . Empirically, $\bar{x} \approx s^2$ for each hour.

E[X] = \text{Var}(X) = \lambda

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(77)
hours = list(range(8, 18))
lambda_map = {8:15, 9:15, 10:25, 11:25, 12:15,
              13:15, 14:25, 15:25, 16:15, 17:8}
records = []
for day in range(1, 61):
    for hour in hours:
        records.append({'hour': hour, 'lambda': lambda_map[hour],
                        'requests': np.random.poisson(lambda_map[hour])})
df = pd.DataFrame(records)

summary = df.groupby('hour')['requests'].agg(
    mean='mean', variance='var', true_lambda='first'
)
summary['true_lambda'] = [lambda_map[h] for h in summary.index]
summary = summary.round(2)
print(summary)

Step 3 — Interpretation

For each hour band the sample mean and sample variance are close to the true $\lambda$ . The ratio $s^2 / \bar{x}$ (dispersion index) should hover near 1.0 for a well-fitting Poisson model.

\text{Dispersion index} = \frac{s^2}{\bar{x}} \approx 1 \quad (\text{Poisson})

Answer

Peak hours 10–11 and 14–15 have

\lambda=25

; off-peak hour 17 has

\lambda=8

. In each band the empirical mean and variance are within

\sim 5\%

\lambda

, confirming the Poisson model.

Waiting Times: Gamma Distribution and Queue Comparison

A service centre has two queues. Standard queue waiting times follow $\text{Gamma}(k=2, \theta=8)$ (mean 16 min). Priority queue follows $\text{Gamma}(k=2, \theta=3)$ (mean 6 min). Compare distributions via empirical CDF, quantiles, and fraction served within 10 minutes.

a)Generate 300 waiting times for each queue.
b)Compute mean, median, Q1, Q3, and the 90th percentile.
c)Compute the fraction served within 10 minutes for each queue.
d)Compare the shapes: same $k$ means same coefficient of variation.

Step 1 — Generate Waiting Times

The gamma distribution with shape $k=2$ is unimodal and right-skewed. Scaling $\theta$ shifts the whole distribution.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(11)
standard = np.random.gamma(shape=2, scale=8, size=300)
priority = np.random.gamma(shape=2, scale=3, size=300)

df = pd.DataFrame({
    'wait_time': np.concatenate([standard, priority]),
    'queue': ['Standard']*300 + ['Priority']*300
})
print(df.groupby('queue')['wait_time'].describe().round(2))

Step 2 — Quantile Comparison

Quantiles directly answer "what fraction of customers wait more than $t$ minutes?" — more useful than the mean alone for service-level agreements.

Q_p = F^{-1}(p), \quad P(X \leq Q_p) = p

dataset generator · Pythonread-only

import numpy as np

np.random.seed(11)
standard = np.random.gamma(2, 8, 300)
priority = np.random.gamma(2, 3, 300)

for label, data in [('Standard', standard), ('Priority', priority)]:
    q = np.percentile(data, [25, 50, 75, 90])
    within10 = (data <= 10).mean()
    print(f"{label}: Q1={q[0]:.1f}, median={q[1]:.1f}, "
          f"Q3={q[2]:.1f}, P90={q[3]:.1f}, "
          f"served<=10min: {within10*100:.1f}%")

Step 3 — Empirical CDF Interpretation

The empirical CDF $\hat{F}(t) = \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}[x_i \leq t]$ estimates the probability of being served within $t$ minutes. The gap between the two CDFs quantifies the priority queue advantage.

\hat{F}(t) = \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}[x_i \leq t]

Answer

Priority queue: median

\approx 4.2

min,

\approx 86\%

served within 10 min. Standard queue: median

\approx 11

min,

\approx 47\%

served within 10 min. Same shape (

k=2

) means

\text{CV} = 1/\sqrt{k} \approx 0.71

for both — relative variability is identical.

Explanation tip: When distributions have the same shape but different scales, comparing medians and quantiles is more informative than comparing means, especially for right-skewed data.

Warehouse Orders: Demand and Return Rates Across Cells

Five warehouses (W1–W5) each stock five product groups (Electronics, Clothing, Food, Furniture, Books). For each warehouse–product cell, we observe total orders and total returns over one quarter. Find the cell with the highest return rate and the product group with the most variable demand.

a)Generate a 5×5 demand matrix and a 5×5 return matrix.
b)Compute the return rate (returns / orders) for each cell.
c)Identify the cell with the highest return rate.
d)Identify the product group with the highest coefficient of variation in demand across warehouses.

Step 1 — Generate Demand and Returns

Demand and return counts are drawn from realistic distributions. Electronics has higher return rates; Food has the most stable demand.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(42)
warehouses = ['W1','W2','W3','W4','W5']
products = ['Electronics','Clothing','Food','Furniture','Books']

base_demand = np.array([
    [500,300,800,200,150],
    [480,320,750,210,160],
    [510,290,820,190,140],
    [490,310,780,205,155],
    [520,280,810,195,145],
], dtype=float)
noise = np.random.normal(0, 20, base_demand.shape)
demand = (base_demand + noise).clip(50).round().astype(int)

return_rate_base = np.array([0.08,0.05,0.02,0.07,0.03])
returns = np.round(demand * (return_rate_base +
          np.random.normal(0, 0.01, demand.shape))).clip(0).astype(int)

df_d = pd.DataFrame(demand, index=warehouses, columns=products)
df_r = pd.DataFrame(returns, index=warehouses, columns=products)
print("Demand:"); print(df_d)
print("Returns:"); print(df_r)

Step 2 — Return Rate Matrix

Return rate per cell: $r_{ij} = \text{returns}_{ij} / \text{orders}_{ij}$ . A high return rate signals quality issues, sizing problems, or unmet expectations.

r_{ij} = \frac{\text{returns}_{ij}}{\text{orders}_{ij}}

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(42)
warehouses = ['W1','W2','W3','W4','W5']
products = ['Electronics','Clothing','Food','Furniture','Books']
base_demand = np.array([[500,300,800,200,150],[480,320,750,210,160],
    [510,290,820,190,140],[490,310,780,205,155],[520,280,810,195,145]],float)
demand = (base_demand + np.random.normal(0,20,base_demand.shape)).clip(50).round().astype(int)
return_rate_base = np.array([0.08,0.05,0.02,0.07,0.03])
returns = np.round(demand*(return_rate_base+np.random.normal(0,0.01,demand.shape))).clip(0).astype(int)
rate = returns / demand
df_rate = pd.DataFrame(rate.round(4), index=warehouses, columns=products)
print("Return rates:"); print(df_rate)
idx = np.unravel_index(rate.argmax(), rate.shape)
print(f"Highest return rate: {warehouses[idx[0]]}/{products[idx[1]]} = {rate[idx]:.3f}")

Step 3 — Coefficient of Variation by Product Group

The coefficient of variation $\text{CV} = s/\bar{x}$ measures relative variability in demand. A high CV means the product group has inconsistent demand across warehouses.

\text{CV} = \frac{s}{\bar{x}}

dataset generator · Pythonread-only

import numpy as np

np.random.seed(42)
products = ['Electronics','Clothing','Food','Furniture','Books']
base = np.array([[500,300,800,200,150],[480,320,750,210,160],
    [510,290,820,190,140],[490,310,780,205,155],[520,280,810,195,145]],float)
demand = (base + np.random.normal(0,20,base.shape)).clip(50).round().astype(int)

for j, prod in enumerate(products):
    col = demand[:, j]
    cv = col.std(ddof=1) / col.mean()
    print(f"{prod}: mean={col.mean():.1f}, std={col.std(ddof=1):.1f}, CV={cv:.4f}")

Answer

Electronics consistently has the highest return rate (

\approx 8\%

). Food has the most stable demand (lowest CV across warehouses). Furniture shows moderate returns; Books have the lowest absolute volume.

Correlation Traps: Why r ≈ 0 Does Not Mean No Relationship

A dataset contains $x \in [-3, 3]$ and $y = x^2 + \varepsilon$ (a perfect quadratic relationship). Compute Pearson $r$ . Observe that $r \approx 0$ despite a clear non-linear relationship. Add three extreme outliers that can flip the sign of $r$ .

a)Generate $x$ uniformly on $[-3, 3]$ and set $y = x^2 + \varepsilon$ .
b)Compute Pearson $r$ and note it is near zero.
c)Add three outliers at $(x, y) = (10, 5), (11, 6), (12, 4)$ and recompute $r$ .
d)Explain why Pearson $r$ only measures LINEAR association.

Step 1 — Generate Quadratic Data

The relationship $y = x^2$ is deterministic and strong, but symmetric around $x=0$ . Positive and negative $x$ values cancel each other in the cross-product sum, driving $r$ toward zero.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(8)
x = np.random.uniform(-3, 3, 200)
y = x**2 + np.random.normal(0, 0.5, 200)

df = pd.DataFrame({'x': x.round(3), 'y': y.round(3)})
print(df.head(8).to_string(index=False))
print(f"Pearson r = {np.corrcoef(x, y)[0,1]:.4f}")
print("Note: |r| near 0 despite strong quadratic pattern!")

Step 2 — Pearson r Formula

Pearson $r$ measures the strength and direction of the LINEAR relationship. For $y=x^2$ , the positive right-half and negative left-half contributions cancel.

r = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2} \cdot \sqrt{\sum_{i=1}^{n}(y_i-\bar{y})^2}}

Step 3 — Effect of Outliers

Adding a few high-leverage outliers can dramatically change $r$ — even reversing its sign. This shows that $r$ is sensitive to outliers.

dataset generator · Pythonread-only

import numpy as np

np.random.seed(8)
x = np.random.uniform(-3, 3, 200)
y = x**2 + np.random.normal(0, 0.5, 200)

x_out = np.append(x, [10, 11, 12])
y_out = np.append(y, [5, 6, 4])

r_before = np.corrcoef(x, y)[0, 1]
r_after = np.corrcoef(x_out, y_out)[0, 1]
print(f"r without outliers: {r_before:.4f}")
print(f"r with outliers:    {r_after:.4f}")
print("Three outliers can shift r by more than 0.3!")

Step 4 — Spearman and Visual Checks

Spearman rank correlation $r_s$ detects monotone non-linear relationships. Always inspect a scatter plot; $r$ alone is insufficient. For the quadratic dataset, $r \approx 0$ but $r_s \approx 0$ too — here even rank correlation fails because the relationship is non-monotone. Use mutual information or polynomial regression instead.

dataset generator · Pythonread-only

import numpy as np
from scipy.stats import spearmanr, pearsonr

np.random.seed(8)
x = np.random.uniform(-3, 3, 200)
y = x**2 + np.random.normal(0, 0.5, 200)

r_p, _ = pearsonr(x, y)
r_s, _ = spearmanr(x, y)
print(f"Pearson r  = {r_p:.4f}  (near 0, misses quadratic)")
print(f"Spearman r = {r_s:.4f}  (also near 0: non-monotone)")
print("Conclusion: ALWAYS plot your data before trusting r.")

Answer

r \approx 0

for

y = x^2 + \varepsilon

on a symmetric domain, despite a deterministic relationship. Pearson

r

only captures linear association. Outliers can shift

r

0.3

or more. Best practice: plot first, then compute correlation as a supplementary measure.

Explanation tip: "Correlation near zero" means "no linear relationship", NOT "no relationship". Always visualise the data and consider non-linear measures for non-linear patterns.

Final Report: The Complete Data Analysis Workflow

Summarise the end-to-end data analysis workflow practiced in this section: load → clean → describe → visualise → model → interpret. For each stage, list the key tools and common pitfalls encountered across Problems 1–10.

Stage 1 — Load Data

Read data from CSV, database, or simulation. Check dtypes, encoding, and missing values immediately after loading. Use `df.info()` and `df.head()` as sanity checks.

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

# Simulate loading a mixed dataset (as in Problems 3-10)
np.random.seed(0)
df = pd.DataFrame({
    'id': range(1, 101),
    'value': np.random.normal(50, 10, 100),
    'group': np.random.choice(['A','B','C'], 100),
    'flag': np.random.randint(0, 2, 100)
})
print(df.info())
print(df.head(5).to_string(index=False))
print(f"Missing values: {df.isnull().sum().sum()}")

Stage 2 — Clean Data

Handle missing values, fix dtypes, remove or flag duplicates, and apply the outlier detection rule from Problem 4.

x_i \text{ is an outlier if } x_i < Q_1 - 1.5\,\text{IQR} \text{ or } x_i > Q_3 + 1.5\,\text{IQR}

dataset generator · Pythonread-only

import numpy as np
import pandas as pd

np.random.seed(0)
values = np.random.normal(50, 10, 100)
values[::25] = np.nan   # inject missing
values[5] = 999          # inject outlier

s = pd.Series(values).dropna()
q1, q3 = s.quantile(0.25), s.quantile(0.75)
iqr = q3 - q1
clean = s[(s >= q1 - 1.5*iqr) & (s <= q3 + 1.5*iqr)]
print(f"Original: n={len(values)}, missing={np.isnan(values).sum()}")
print(f"After cleaning: n={len(clean)}, "
      f"mean={clean.mean():.2f}, std={clean.std():.2f}")

Stage 3 — Describe

Compute the standard descriptive statistics. Use all five: mean, median, std, IQR, and range. Do not rely on a single number.

\bar{x} = \frac{1}{n}\sum x_i, \quad s^2 = \frac{1}{n-1}\sum(x_i-\bar{x})^2, \quad \text{IQR} = Q_3 - Q_1

Stage 4 — Visualise

A histogram reveals shape; a box-plot shows outliers and IQR; a scatter plot reveals relationships. As Problem 10 showed: always plot before computing $r$ . As Problems 3 and 8 showed: overlay distributions when comparing groups.

Stage 5 — Model and Test

Fit the appropriate model: Poisson for counts (Problem 7), lognormal for skewed positives (Problem 4), gamma for waiting times (Problem 8), normal for exam scores (Problem 3). Validate with a goodness-of-fit test (Problem 1 used the chi-square intuition).

\chi^2 = \sum_{k} \frac{(O_k - E_k)^2}{E_k}

Stage 6 — Interpret with Caveats

State conclusions in plain language with uncertainty bounds. Highlight key caveats: correlation $\neq$ causation; $r \approx 0$ $\neq$ no relationship (Problem 10); higher mean $\neq$ better outcome when variance differs (Problem 3); mean $\neq$ typical value for skewed distributions (Problem 4).

dataset generator · Pythonread-only

# Summary template: print a mini-report for any numeric column
import numpy as np
import pandas as pd

np.random.seed(1)
data = np.concatenate([
    np.random.lognormal(3.5, 0.4, 195),
    np.random.uniform(300, 500, 5)
])
s = pd.Series(data, name='metric')
q1, q3 = s.quantile(0.25), s.quantile(0.75)
print("=== Mini Analysis Report ===")
print(f"n         : {len(s)}")
print(f"Mean      : {s.mean():.2f}")
print(f"Median    : {s.median():.2f}")
print(f"Std       : {s.std():.2f}")
print(f"IQR       : {q3-q1:.2f}")
print(f"Skewness  : {s.skew():.3f}")
print("Caveats: distribution is right-skewed; "
      "median is a better centre measure than mean.")

Answer

A rigorous analysis always follows all six stages. Skipping any stage risks misleading conclusions. The ten problems in this section each illustrated a different failure mode: convergence illusions (P1–P2), group heterogeneity (P3), outlier sensitivity (P4), conditional vs marginal rates (P5), bias-variance in manufacturing (P6), distributional model validation (P7–P8), multi-dimensional aggregation (P9), and the limits of Pearson correlation (P10).

Explanation tip: Data analysis is iterative. After interpreting results, revisit earlier stages as new questions arise. Document every transformation so the analysis is reproducible.