Lab — Statistical Foundations for Data Science

15 experiments, done in five Python programs, each set out as 1. Question, 2. Aim, 3. Steps,

4. Programme, 5. Execution and Results.

The prescribed lab is headed "Advanced Spreadsheets/Excel Lab/PSPP Open Source". Every experiment is a spreadsheet exercise, and the practical exam tests it in a spreadsheet.

Do each experiment twice

Once in Excel — that is what the exam marks. Once in Python — that is what the course is for.

The prescribed lab never uses Python, even though Python Programming and Data Structures teaches the language; Python-based analysis is taught in Python for Data Analysis and Visualization. That gap is worth closing yourself. See SYLLABUS-REVIEW.md finding D8.

bash tools/data-science/run_stats_labs.sh                       # verify everything
cd labs/course-4-stats/python && python3 test_statlib.py

The Python programs share one file, so the experiments are set out below by program, in the order of their first experiment; each program's steps follow its experiments.

The experiments

Experiment 2 is reconstructed

The official text of experiment 2 survives only as the fragment "a positive result." — the question stem is missing from the PDF. See review findings D1 and D3.

Reconstructed question:

Answer: about 16.7%, not 99%. Of 10,000 people, 99 true positives are outnumbered by 495 false positives. Build that 10,000-person table in your sheet — it makes the result obvious and it earns marks.

This is the base rate fallacy, and Bayes' theorem is examined despite being absent from the syllabus units. Do not skip it.

Experiments 1 and 2 — Contingency tables and Bayes' theorem

1. Question

Unit 1. Experiment 1: from sales data, build the contingency table of region against purchase type, and find the joint, marginal and conditional probabilities; are region and purchase independent? Experiment 2 (reconstructed): a disease affects 1% of a population, and a test is 99% sensitive and 95% specific; given a positive result, what is the probability that the person has the disease?

2. Aim

Read probabilities off a contingency table, and turn a test's accuracy into the chance that a positive result is right.

3. Steps

1. Experiment 1: the contingency table.
2. Joint probabilities.
3. Marginal probabilities.
4. Conditional probabilities.
5. The independence check.
6. Experiment 2: the medical test.
7. The law of total probability.
8. Bayes' theorem.

4. Programme

"""Course 4 Lab, experiments 1-2: contingency tables, conditional probability,
independence, and Bayes' theorem.

Experiment 2 is reconstructed -- the official text survives only as the
fragment "a positive result." See SYLLABUS-REVIEW.md findings D1 and D3.

Standard library only: no numpy, no pandas, no scipy.
"""

from fractions import Fraction

# ---------------------------------------------------------------------
# EXPERIMENT 1: contingency table, conditional probabilities, independence
# ---------------------------------------------------------------------
# Step 1: Experiment 1: the contingency table
print("=" * 66)
print("EXPERIMENT 1: Contingency table from sales data")
print("=" * 66)

# Rows = region, columns = whether the customer bought the premium product.
table = {
    "North": {"Premium": 30, "Standard": 70},
    "South": {"Premium": 45, "Standard": 55},
    "East":  {"Premium": 25, "Standard": 75},
}
columns = ["Premium", "Standard"]

row_totals = {r: sum(cells.values()) for r, cells in table.items()}
col_totals = {c: sum(table[r][c] for r in table) for c in columns}
grand_total = sum(row_totals.values())

print(f"\n{'Region':<10}" + "".join(f"{c:>12}" for c in columns) + f"{'Total':>10}")
print("-" * 46)
for region, cells in table.items():
    print(f"{region:<10}" + "".join(f"{cells[c]:>12}" for c in columns)
          + f"{row_totals[region]:>10}")
print("-" * 46)
print(f"{'Total':<10}" + "".join(f"{col_totals[c]:>12}" for c in columns)
      + f"{grand_total:>10}")

# Step 2: Joint probabilities
print("\nJOINT probabilities  P(Region and Purchase) = cell / grand total")
for region, cells in table.items():
    for col in columns:
        print(f"  P({region} and {col:<8}) = {cells[col]:>3}/{grand_total} "
              f"= {cells[col] / grand_total:.4f}")

# Step 3: Marginal probabilities
print("\nMARGINAL probabilities")
for region in table:
    print(f"  P({region:<6}) = {row_totals[region]}/{grand_total} "
          f"= {row_totals[region] / grand_total:.4f}")
for col in columns:
    print(f"  P({col:<8}) = {col_totals[col]}/{grand_total} "
          f"= {col_totals[col] / grand_total:.4f}")

# Step 4: Conditional probabilities
print("\nCONDITIONAL probabilities  P(A|B) = P(A and B) / P(B)")
for region in table:
    joint = table[region]["Premium"] / grand_total
    marginal = row_totals[region] / grand_total
    print(f"  P(Premium | {region:<6}) = {joint:.4f} / {marginal:.4f} "
          f"= {joint / marginal:.4f}")

# Step 5: The independence check
print("\nINDEPENDENCE CHECK")
print("  A and B are independent if P(A and B) = P(A) x P(B)")
for region in table:
    joint = table[region]["Premium"] / grand_total
    product = (row_totals[region] / grand_total) * (col_totals["Premium"] / grand_total)
    verdict = "independent" if abs(joint - product) < 1e-9 else "NOT independent"
    print(f"  {region:<6}: P(joint) = {joint:.4f}, P(A)xP(B) = {product:.4f}"
          f"  -> {verdict}")
print("\n  Conclusion: region and purchase type are dependent -- knowing the")
print("  region changes the probability of a premium purchase.")

# ---------------------------------------------------------------------
# EXPERIMENT 2 [RECONSTRUCTED]: Bayes' theorem, medical testing
# ---------------------------------------------------------------------
# Step 6: Experiment 2: the medical test
print("\n" + "=" * 66)
print("EXPERIMENT 2 [RECONSTRUCTED]: Bayes' theorem")
print("=" * 66)
print("""
Question (reconstructed from the surviving fragment "a positive result."):

  A disease affects 1% of a population. A test for it is 99% sensitive
  (it correctly flags 99% of people who have the disease) and 95% specific
  (it correctly clears 95% of people who do not).

  A randomly chosen person tests positive. What is the probability that
  they actually have the disease, given a positive result?
""")

p_disease = 0.01                 # prior P(D)
p_no_disease = 1 - p_disease
sensitivity = 0.99               # P(+ | D)
specificity = 0.95               # P(- | not D)
false_positive = 1 - specificity  # P(+ | not D)

# Step 7: The law of total probability
# Total probability of a positive result (the denominator).
p_positive = sensitivity * p_disease + false_positive * p_no_disease
posterior = sensitivity * p_disease / p_positive

print("  Given:")
print(f"    P(D)          = {p_disease}          prior probability of disease")
print(f"    P(+ | D)      = {sensitivity}          sensitivity")
print(f"    P(- | not D)  = {specificity}          specificity")
print(f"    P(+ | not D)  = {false_positive:.2f}          false positive rate")

print("\n  Law of total probability -- P(+):")
print(f"    P(+) = P(+|D)P(D) + P(+|not D)P(not D)")
print(f"         = {sensitivity} x {p_disease} + {false_positive:.2f} x {p_no_disease}")
print(f"         = {sensitivity * p_disease:.4f} + {false_positive * p_no_disease:.4f}"
      f" = {p_positive:.4f}")

# Step 8: Bayes' theorem
print("\n  Bayes' theorem -- P(D | +):")
print(f"    P(D|+) = P(+|D)P(D) / P(+)")
print(f"           = {sensitivity * p_disease:.4f} / {p_positive:.4f}")
print(f"           = {posterior:.4f}  ->  about {posterior * 100:.1f}%")

print(f"""
  THE POINT OF THIS QUESTION: the test is 99% sensitive, yet a positive
  result means only a {posterior * 100:.0f}% chance of having the disease. Because the
  disease is rare, the 5% of false positives drawn from the 99% who are
  healthy vastly outnumber the true positives.

  Out of 10,000 people:
    {int(10000 * p_disease)} have the disease, of whom {10000 * p_disease * sensitivity:.0f} test positive
    {int(10000 * p_no_disease)} do not, of whom {10000 * p_no_disease * false_positive:.0f} still test positive
    so {10000 * p_positive:.0f} positives in total, of which only {10000 * p_disease * sensitivity:.0f} are real.

  Confusing P(D|+) with P(+|D) is called the base rate fallacy, and it is
  the single most examined idea in this part of the syllabus.
""")

print("  Exact fraction:", Fraction(posterior).limit_denominator(10000))

5. Execution and Results

==================================================================
EXPERIMENT 1: Contingency table from sales data
==================================================================

Region         Premium    Standard     Total
----------------------------------------------
North               30          70       100
South               45          55       100
East                25          75       100
----------------------------------------------
Total              100         200       300

JOINT probabilities  P(Region and Purchase) = cell / grand total
  P(North and Premium ) =  30/300 = 0.1000
  P(North and Standard) =  70/300 = 0.2333
  P(South and Premium ) =  45/300 = 0.1500
  P(South and Standard) =  55/300 = 0.1833
  P(East and Premium ) =  25/300 = 0.0833
  P(East and Standard) =  75/300 = 0.2500

MARGINAL probabilities
  P(North ) = 100/300 = 0.3333
  P(South ) = 100/300 = 0.3333
  P(East  ) = 100/300 = 0.3333
  P(Premium ) = 100/300 = 0.3333
  P(Standard) = 200/300 = 0.6667

CONDITIONAL probabilities  P(A|B) = P(A and B) / P(B)
  P(Premium | North ) = 0.1000 / 0.3333 = 0.3000
  P(Premium | South ) = 0.1500 / 0.3333 = 0.4500
  P(Premium | East  ) = 0.0833 / 0.3333 = 0.2500

INDEPENDENCE CHECK
  A and B are independent if P(A and B) = P(A) x P(B)
  North : P(joint) = 0.1000, P(A)xP(B) = 0.1111  -> NOT independent
  South : P(joint) = 0.1500, P(A)xP(B) = 0.1111  -> NOT independent
  East  : P(joint) = 0.0833, P(A)xP(B) = 0.1111  -> NOT independent

  Conclusion: region and purchase type are dependent -- knowing the
  region changes the probability of a premium purchase.

==================================================================
EXPERIMENT 2 [RECONSTRUCTED]: Bayes' theorem
==================================================================

Question (reconstructed from the surviving fragment "a positive result."):

  A disease affects 1% of a population. A test for it is 99% sensitive
  (it correctly flags 99% of people who have the disease) and 95% specific
  (it correctly clears 95% of people who do not).

  A randomly chosen person tests positive. What is the probability that
  they actually have the disease, given a positive result?

  Given:
    P(D)          = 0.01          prior probability of disease
    P(+ | D)      = 0.99          sensitivity
    P(- | not D)  = 0.95          specificity
    P(+ | not D)  = 0.05          false positive rate

  Law of total probability -- P(+):
    P(+) = P(+|D)P(D) + P(+|not D)P(not D)
         = 0.99 x 0.01 + 0.05 x 0.99
         = 0.0099 + 0.0495 = 0.0594

  Bayes' theorem -- P(D | +):
    P(D|+) = P(+|D)P(D) / P(+)
           = 0.0099 / 0.0594
           = 0.1667  ->  about 16.7%

  THE POINT OF THIS QUESTION: the test is 99% sensitive, yet a positive
  result means only a 17% chance of having the disease. Because the
  disease is rare, the 5% of false positives drawn from the 99% who are
  healthy vastly outnumber the true positives.

  Out of 10,000 people:
    100 have the disease, of whom 99 test positive
    9900 do not, of whom 495 still test positive
    so 594 positives in total, of which only 99 are real.

  Confusing P(D|+) with P(+|D) is called the base rate fallacy, and it is
  the single most examined idea in this part of the syllabus.

  Exact fraction: 1/6

The dataset used in the Python version gives χ² = 9.75 on 2 df (p = 0.0076), so region and purchase type are genuinely associated. Experiments 1 and 15 are answering the same question at two levels of rigour.

Experiments 3 to 6 — Central tendency, dispersion, the histogram and the bar chart

1. Question

Unit 1. For a set of marks: 3, find the mean, median and mode; 4, the range, quartiles, variance, standard deviation and coefficient of variation; 5, draw a histogram and describe the shape of the distribution; 6, draw a bar chart of categorical data.

2. Aim

Summarise a data set by its centre, its spread and its shape, and know which summary to quote.

3. Steps

1. Experiment 3: the mean.
2. The median.
3. The mode.
4. Which one to use.
5. Experiment 4: range and quartiles.
6. Variance and standard deviation.
7. Coefficient of variation, and outliers.
8. Experiment 5: the histogram and its shape.
9. Experiment 6: a bar chart of categories.

4. Programme

"""Course 4 Lab, experiments 3-6: measures of central tendency and dispersion,
histograms and bar charts.

Standard library only -- `statistics` is part of Python, no install needed.
Charts are drawn as text so they render anywhere.
"""

import statistics
from collections import Counter

marks = [45, 67, 78, 52, 89, 91, 73, 64, 58, 82,
         76, 69, 71, 85, 60, 55, 93, 48, 79, 66]

# Step 1: Experiment 3: the mean
print("=" * 62)
print("EXPERIMENT 3: Measures of central tendency")
print("=" * 62)
print(f"Dataset (n = {len(marks)}): {sorted(marks)}\n")

n = len(marks)
mean = sum(marks) / n
ordered = sorted(marks)
if n % 2 == 0:
    median = (ordered[n // 2 - 1] + ordered[n // 2]) / 2
else:
    median = ordered[n // 2]

counts = Counter(marks)
top = max(counts.values())
modes = [v for v, c in counts.items() if c == top]

print("MEAN -- the balance point. Add everything, divide by how many.")
print(f"  mean = {sum(marks)} / {n} = {mean:.2f}")
print(f"  cross-check statistics.mean() = {statistics.mean(marks):.2f}")

# Step 2: The median
print("\nMEDIAN -- the middle value once sorted. Half are below, half above.")
print(f"  n = {n} is even, so average the {n//2}th and {n//2+1}th values:")
print(f"  ({ordered[n//2-1]} + {ordered[n//2]}) / 2 = {median}")
print(f"  cross-check statistics.median() = {statistics.median(marks)}")

# Step 3: The mode
print("\nMODE -- the most frequent value.")
if top == 1:
    print("  every value appears once, so there is no mode")
else:
    print(f"  {modes} (appearing {top} times)")

# Step 4: Which one to use
print("\nWHICH ONE TO USE")
print("  The mean uses every value, so one extreme value drags it. Add a")
print("  single mark of 500 to this dataset:")
skewed = marks + [500]
print(f"    mean   {mean:.2f} -> {sum(skewed)/len(skewed):.2f}   moved a lot")
print(f"    median {median} -> {statistics.median(skewed)}   barely moved")
print("  That resistance is why income and house prices are quoted as medians.")

# Step 5: Experiment 4: range and quartiles
print("\n" + "=" * 62)
print("EXPERIMENT 4: Measures of dispersion")
print("=" * 62)

data_range = max(marks) - min(marks)


def quartile(sorted_data, q):
    """Linear-interpolation quartile, matching Excel's QUARTILE.INC."""
    pos = (len(sorted_data) - 1) * q
    lower = int(pos)
    upper = min(lower + 1, len(sorted_data) - 1)
    return sorted_data[lower] + (pos - lower) * (sorted_data[upper] - sorted_data[lower])


q1, q2, q3 = quartile(ordered, 0.25), quartile(ordered, 0.5), quartile(ordered, 0.75)
iqr = q3 - q1

# Population vs sample: the divisor is n for a population, n-1 for a sample.
pop_var = sum((x - mean) ** 2 for x in marks) / n
sam_var = sum((x - mean) ** 2 for x in marks) / (n - 1)

print(f"RANGE = max - min = {max(marks)} - {min(marks)} = {data_range}")
print("  Uses only two values, so a single outlier defines it entirely.")

print(f"\nQUARTILES and IQR")
print(f"  Q1 = {q1:.2f}   Q2 (median) = {q2:.2f}   Q3 = {q3:.2f}")
print(f"  IQR = Q3 - Q1 = {q3:.2f} - {q1:.2f} = {iqr:.2f}")
print("  The IQR is the spread of the middle 50%, so outliers cannot inflate it.")

# Step 6: Variance and standard deviation
print(f"\nVARIANCE -- the mean squared distance from the mean")
print(f"  population variance (divide by n)     = {pop_var:.2f}")
print(f"  sample variance     (divide by n - 1) = {sam_var:.2f}")
print("  Use n-1 when the data is a SAMPLE. Dividing by n underestimates the")
print("  spread, because deviations are measured from the sample's own mean.")
print("  This is Bessel's correction, and choosing wrongly costs marks.")
print(f"  cross-check statistics.pvariance() = {statistics.pvariance(marks):.2f}")
print(f"  cross-check statistics.variance()  = {statistics.variance(marks):.2f}")

print(f"\nSTANDARD DEVIATION -- the square root of the variance")
print(f"  population sd = {pop_var ** 0.5:.2f}")
print(f"  sample sd     = {sam_var ** 0.5:.2f}")
print("  Back in the original units (marks), unlike variance (marks squared).")

# Step 7: Coefficient of variation, and outliers
print(f"\nCOEFFICIENT OF VARIATION -- relative spread, unit-free")
print(f"  CV = sd / mean x 100 = {sam_var ** 0.5 / mean * 100:.2f}%")

print("\nOUTLIER RULE: anything below Q1 - 1.5xIQR or above Q3 + 1.5xIQR")
low_fence, high_fence = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print(f"  fences: [{low_fence:.2f}, {high_fence:.2f}]")
outliers = [x for x in marks if x < low_fence or x > high_fence]
print(f"  outliers: {outliers if outliers else 'none'}")

# Step 8: Experiment 5: the histogram and its shape
print("\n" + "=" * 62)
print("EXPERIMENT 5: Histogram and the shape of the distribution")
print("=" * 62)

bins = [(40, 50), (50, 60), (60, 70), (70, 80), (80, 90), (90, 100)]
print(f"\n{'Class':<12}{'Frequency':<12}Histogram")
for low, high in bins:
    freq = sum(1 for m in marks if low <= m < high)
    print(f"{low}-{high:<8} {freq:<12}{'#' * freq * 3}")

# Pearson's second coefficient of skewness.
skew = 3 * (mean - median) / (sam_var ** 0.5)
print(f"\nSHAPE")
print(f"  mean = {mean:.2f}, median = {median}")
print(f"  Pearson skewness = 3(mean - median)/sd = {skew:.3f}")
if abs(skew) < 0.5:
    shape = "roughly symmetric"
elif skew > 0:
    shape = "positively skewed -- a tail to the right"
else:
    shape = "negatively skewed -- a tail to the left"
print(f"  The distribution is {shape}.")
print("  Rule of thumb: mean > median suggests a right tail; mean < median a")
print("  left tail; mean = median = mode means perfectly symmetric.")

# Step 9: Experiment 6: a bar chart of categories
print("\n" + "=" * 62)
print("EXPERIMENT 6: Bar chart of categorical data")
print("=" * 62)

survey = [("Male", "A"), ("Female", "A"), ("Male", "B"), ("Female", "B"),
          ("Male", "A"), ("Female", "C"), ("Male", "C"), ("Female", "A"),
          ("Male", "B"), ("Female", "B"), ("Male", "A"), ("Female", "A"),
          ("Male", "C"), ("Female", "B"), ("Male", "B")]

sections = sorted({s for _, s in survey})
genders = sorted({g for g, _ in survey})

print(f"\n{'Section':<10}" + "".join(f"{g:>10}" for g in genders) + f"{'Total':>10}")
print("-" * 40)
for section in sections:
    row = [sum(1 for g, s in survey if s == section and g == gender)
           for gender in genders]
    print(f"{section:<10}" + "".join(f"{v:>10}" for v in row) + f"{sum(row):>10}")

print("\nGrouped bar chart")
for section in sections:
    for gender in genders:
        count = sum(1 for g, s in survey if s == section and g == gender)
        print(f"  {section}-{gender:<8} {'|' * count * 4} {count}")

print("\nINTERPRETATION")
print("  A bar chart is for CATEGORIES -- the bars are separated, and their")
print("  order carries no meaning. A histogram is for CONTINUOUS data -- the")
print("  bars touch, because the classes are adjacent intervals. Drawing one")
print("  when the question asks for the other is a common way to lose marks.")

5. Execution and Results

==============================================================
EXPERIMENT 3: Measures of central tendency
==============================================================
Dataset (n = 20): [45, 48, 52, 55, 58, 60, 64, 66, 67, 69, 71, 73, 76, 78, 79, 82, 85, 89, 91, 93]

MEAN -- the balance point. Add everything, divide by how many.
  mean = 1401 / 20 = 70.05
  cross-check statistics.mean() = 70.05

MEDIAN -- the middle value once sorted. Half are below, half above.
  n = 20 is even, so average the 10th and 11th values:
  (69 + 71) / 2 = 70.0
  cross-check statistics.median() = 70.0

MODE -- the most frequent value.
  every value appears once, so there is no mode

WHICH ONE TO USE
  The mean uses every value, so one extreme value drags it. Add a
  single mark of 500 to this dataset:
    mean   70.05 -> 90.52   moved a lot
    median 70.0 -> 71   barely moved
  That resistance is why income and house prices are quoted as medians.

==============================================================
EXPERIMENT 4: Measures of dispersion
==============================================================
RANGE = max - min = 93 - 45 = 48
  Uses only two values, so a single outlier defines it entirely.

QUARTILES and IQR
  Q1 = 59.50   Q2 (median) = 70.00   Q3 = 79.75
  IQR = Q3 - Q1 = 79.75 - 59.50 = 20.25
  The IQR is the spread of the middle 50%, so outliers cannot inflate it.

VARIANCE -- the mean squared distance from the mean
  population variance (divide by n)     = 192.75
  sample variance     (divide by n - 1) = 202.89
  Use n-1 when the data is a SAMPLE. Dividing by n underestimates the
  spread, because deviations are measured from the sample's own mean.
  This is Bessel's correction, and choosing wrongly costs marks.
  cross-check statistics.pvariance() = 192.75
  cross-check statistics.variance()  = 202.89

STANDARD DEVIATION -- the square root of the variance
  population sd = 13.88
  sample sd     = 14.24
  Back in the original units (marks), unlike variance (marks squared).

COEFFICIENT OF VARIATION -- relative spread, unit-free
  CV = sd / mean x 100 = 20.33%

OUTLIER RULE: anything below Q1 - 1.5xIQR or above Q3 + 1.5xIQR
  fences: [29.12, 110.12]
  outliers: none

==============================================================
EXPERIMENT 5: Histogram and the shape of the distribution
==============================================================

Class       Frequency   Histogram
40-50       2           ######
50-60       3           #########
60-70       5           ###############
70-80       5           ###############
80-90       3           #########
90-100      2           ######

SHAPE
  mean = 70.05, median = 70.0
  Pearson skewness = 3(mean - median)/sd = 0.011
  The distribution is roughly symmetric.
  Rule of thumb: mean > median suggests a right tail; mean < median a
  left tail; mean = median = mode means perfectly symmetric.

==============================================================
EXPERIMENT 6: Bar chart of categorical data
==============================================================

Section       Female      Male     Total
----------------------------------------
A                  3         3         6
B                  3         3         6
C                  1         2         3

Grouped bar chart
  A-Female   |||||||||||| 3
  A-Male     |||||||||||| 3
  B-Female   |||||||||||| 3
  B-Male     |||||||||||| 3
  C-Female   |||| 1
  C-Male     |||||||| 2

INTERPRETATION
  A bar chart is for CATEGORIES -- the bars are separated, and their
  order carries no meaning. A histogram is for CONTINUOUS data -- the
  bars touch, because the classes are adjacent intervals. Drawing one
  when the question asks for the other is a common way to lose marks.

Experiments 7, 12 and 13 — Covariance, correlation and regression

1. Question

Units 1 and 4. For hours studied against exam score: 7, draw the scatter plot and find the covariance and correlation; 12, find Pearson's and Spearman's correlation; 13, fit the simple linear regression of score on hours, and test it.

2. Aim

Measure how closely two variables move together, and fit the line that predicts one from the other.

3. Steps

1. The data, and the deviations from the means.
2. Covariance.
3. Pearson's r.
4. Spearman's rank correlation.
5. Experiment 13: the regression line.
6. Residuals.
7. The analysis of variance, and R squared.
8. Testing the slope.
9. Prediction.

4. Programme

"""Course 4 Lab, experiments 7, 12 and 13: covariance, Pearson and Spearman
correlation, and simple linear regression with ANOVA.
"""

import statlib as S

# Step 1: The data, and the deviations from the means
# The classic paired dataset: hours studied against exam score.
hours = [2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
scores = [52, 55, 61, 64, 70, 72, 78, 82, 85, 91]
n = len(hours)

print("=" * 66)
print("EXPERIMENT 7 and 12: Covariance and correlation")
print("=" * 66)

mean_x = sum(hours) / n
mean_y = sum(scores) / n

print(f"\n{'Hours (x)':<12}{'Score (y)':<12}{'x-mx':<10}{'y-my':<10}"
      f"{'(x-mx)(y-my)':<15}{'(x-mx)^2':<12}{'(y-my)^2'}")
print("-" * 82)
sum_xy = sum_xx = sum_yy = 0.0
for x, y in zip(hours, scores):
    dx, dy = x - mean_x, y - mean_y
    sum_xy += dx * dy
    sum_xx += dx * dx
    sum_yy += dy * dy
    print(f"{x:<12}{y:<12}{dx:<10.2f}{dy:<10.2f}{dx * dy:<15.2f}"
          f"{dx * dx:<12.2f}{dy * dy:.2f}")
print("-" * 82)
print(f"{'Mean':<12}{mean_x:<12.2f}{'':<10}{'':<10}"
      f"{sum_xy:<15.2f}{sum_xx:<12.2f}{sum_yy:.2f}")
print(f"  mean of y = {mean_y:.2f}")

cov_sample = sum_xy / (n - 1)
cov_pop = sum_xy / n
# Step 2: Covariance
print(f"\nCOVARIANCE")
print(f"  sample     (divide by n-1) = {sum_xy:.2f} / {n-1} = {cov_sample:.3f}")
print(f"  population (divide by n)   = {sum_xy:.2f} / {n} = {cov_pop:.3f}")
print("  Sign tells direction, magnitude is meaningless -- it depends on the")
print("  units. Change hours to minutes and the covariance multiplies by 60.")

r = sum_xy / (sum_xx * sum_yy) ** 0.5
# Step 3: Pearson's r
print(f"\nPEARSON CORRELATION r")
print(f"  r = sum((x-mx)(y-my)) / sqrt(sum(x-mx)^2 x sum(y-my)^2)")
print(f"    = {sum_xy:.2f} / sqrt({sum_xx:.2f} x {sum_yy:.2f})")
print(f"    = {sum_xy:.2f} / {(sum_xx * sum_yy) ** 0.5:.2f}")
print(f"    = {r:.4f}")
print("  r is unit-free and always lies between -1 and +1. Excel: =CORREL(x,y)")

strength = ("very strong" if abs(r) >= 0.9 else "strong" if abs(r) >= 0.7
            else "moderate" if abs(r) >= 0.4 else "weak")
direction = "positive" if r > 0 else "negative"
print(f"  Interpretation: a {strength} {direction} linear relationship.")


# Step 4: Spearman's rank correlation
def rank(values):
    """Ranks with ties averaged -- what Spearman requires."""
    ordered = sorted(range(len(values)), key=lambda i: values[i])
    ranks = [0.0] * len(values)
    i = 0
    while i < len(ordered):
        j = i
        while j + 1 < len(ordered) and values[ordered[j + 1]] == values[ordered[i]]:
            j += 1
        average_rank = (i + j) / 2 + 1
        for k in range(i, j + 1):
            ranks[ordered[k]] = average_rank
        i = j + 1
    return ranks


rx, ry = rank(hours), rank(scores)
d_squared = sum((a - b) ** 2 for a, b in zip(rx, ry))
rho = 1 - (6 * d_squared) / (n * (n * n - 1))

print(f"\nSPEARMAN RANK CORRELATION rho")
print(f"  {'x':<8}{'y':<8}{'rank x':<10}{'rank y':<10}{'d':<8}{'d^2'}")
for x, y, a, b in zip(hours, scores, rx, ry):
    print(f"  {x:<8}{y:<8}{a:<10.1f}{b:<10.1f}{a - b:<8.1f}{(a - b) ** 2:.2f}")
print(f"  sum of d^2 = {d_squared:.2f}")
print(f"  rho = 1 - 6.sum(d^2) / n(n^2-1) = 1 - 6({d_squared:.0f}) / "
      f"{n}({n * n - 1}) = {rho:.4f}")
print("  Spearman works on ranks, so it detects any MONOTONIC relationship,")
print("  not only a straight-line one, and is unaffected by outliers.")

# Step 5: Experiment 13: the regression line
print("\n" + "=" * 66)
print("EXPERIMENT 13: Simple linear regression")
print("=" * 66)

b1 = sum_xy / sum_xx                 # slope
b0 = mean_y - b1 * mean_x            # intercept

print(f"\n  Model: y = b0 + b1.x")
print(f"  b1 (slope)     = sum((x-mx)(y-my)) / sum((x-mx)^2)")
print(f"                 = {sum_xy:.2f} / {sum_xx:.2f} = {b1:.4f}")
print(f"  b0 (intercept) = my - b1.mx = {mean_y:.2f} - {b1:.4f} x {mean_x:.2f}"
      f" = {b0:.4f}")
print(f"\n  FITTED LINE:  y = {b0:.4f} + {b1:.4f}x")
print(f"  Meaning: each extra hour of study is associated with about "
      f"{b1:.2f} more marks.")
print(f"  The intercept {b0:.2f} is the predicted score at zero hours -- treat")
print("  it cautiously, since x = 0 lies outside the observed range.")

# Step 6: Residuals
print(f"\n  RESIDUALS")
print(f"  {'x':<8}{'observed y':<14}{'fitted y':<14}{'residual':<12}{'residual^2'}")
ss_res = ss_tot = 0.0
for x, y in zip(hours, scores):
    fitted = b0 + b1 * x
    residual = y - fitted
    ss_res += residual ** 2
    ss_tot += (y - mean_y) ** 2
    print(f"  {x:<8}{y:<14}{fitted:<14.3f}{residual:<12.3f}{residual ** 2:.4f}")

ss_reg = ss_tot - ss_res
r_squared = ss_reg / ss_tot

# Step 7: The analysis of variance, and R squared
print(f"\n  ANALYSIS OF VARIANCE")
print(f"  {'Source':<14}{'SS':<14}{'df':<8}{'MS':<14}{'F'}")
print("  " + "-" * 58)
df_reg, df_res = 1, n - 2
ms_reg, ms_res = ss_reg / df_reg, ss_res / df_res
f_stat = ms_reg / ms_res
print(f"  {'Regression':<14}{ss_reg:<14.4f}{df_reg:<8}{ms_reg:<14.4f}{f_stat:.4f}")
print(f"  {'Residual':<14}{ss_res:<14.4f}{df_res:<8}{ms_res:<14.4f}")
print(f"  {'Total':<14}{ss_tot:<14.4f}{n - 1:<8}")

p_value = S.f_sf(f_stat, df_reg, df_res)
print(f"\n  F = {f_stat:.4f} on ({df_reg}, {df_res}) df,  p = {p_value:.3e}")
print(f"  R^2 = SS_regression / SS_total = {ss_reg:.4f} / {ss_tot:.4f} = "
      f"{r_squared:.4f}")
print(f"  So {r_squared * 100:.2f}% of the variation in scores is explained by hours.")
print(f"  Note that R^2 = r^2: {r:.4f}^2 = {r ** 2:.4f}  -- true for SIMPLE")
print("  regression only, not for multiple regression.")

se_slope = (ms_res / sum_xx) ** 0.5
t_stat = b1 / se_slope
# Step 8: Testing the slope
print(f"\n  TESTING THE SLOPE   H0: b1 = 0   H1: b1 != 0")
print(f"  standard error of b1 = sqrt(MS_res / sum(x-mx)^2) = {se_slope:.4f}")
print(f"  t = b1 / se = {b1:.4f} / {se_slope:.4f} = {t_stat:.4f} on {df_res} df")
print(f"  two-tailed p = {S.t_sf_two_tailed(t_stat, df_res):.3e}")
print(f"  p < 0.05, so reject H0 -- the slope is significantly different from 0.")
print(f"  For simple regression t^2 = F: {t_stat:.4f}^2 = {t_stat ** 2:.4f}"
      f" = {f_stat:.4f}")

# Step 9: Prediction
print(f"\n  PREDICTION")
for x in (7.5, 12):
    note = "" if min(hours) <= x <= max(hours) else "   <- EXTRAPOLATION, unsafe"
    print(f"    x = {x:<5} -> predicted score {b0 + b1 * x:.2f}{note}")

print("\n  CORRELATION IS NOT CAUSATION. This data cannot show that studying")
print("  causes higher scores -- only that they move together. A confounder")
print("  (prior ability, motivation) could drive both.")

5. Execution and Results

==================================================================
EXPERIMENT 7 and 12: Covariance and correlation
==================================================================

Hours (x)   Score (y)   x-mx      y-my      (x-mx)(y-my)   (x-mx)^2    (y-my)^2
----------------------------------------------------------------------------------
2           52          -4.50     -19.00    85.50          20.25       361.00
3           55          -3.50     -16.00    56.00          12.25       256.00
4           61          -2.50     -10.00    25.00          6.25        100.00
5           64          -1.50     -7.00     10.50          2.25        49.00
6           70          -0.50     -1.00     0.50           0.25        1.00
7           72          0.50      1.00      0.50           0.25        1.00
8           78          1.50      7.00      10.50          2.25        49.00
9           82          2.50      11.00     27.50          6.25        121.00
10          85          3.50      14.00     49.00          12.25       196.00
11          91          4.50      20.00     90.00          20.25       400.00
----------------------------------------------------------------------------------
Mean        6.50                            355.00         82.50       1534.00
  mean of y = 71.00

COVARIANCE
  sample     (divide by n-1) = 355.00 / 9 = 39.444
  population (divide by n)   = 355.00 / 10 = 35.500
  Sign tells direction, magnitude is meaningless -- it depends on the
  units. Change hours to minutes and the covariance multiplies by 60.

PEARSON CORRELATION r
  r = sum((x-mx)(y-my)) / sqrt(sum(x-mx)^2 x sum(y-my)^2)
    = 355.00 / sqrt(82.50 x 1534.00)
    = 355.00 / 355.75
    = 0.9979
  r is unit-free and always lies between -1 and +1. Excel: =CORREL(x,y)
  Interpretation: a very strong positive linear relationship.

SPEARMAN RANK CORRELATION rho
  x       y       rank x    rank y    d       d^2
  2       52      1.0       1.0       0.0     0.00
  3       55      2.0       2.0       0.0     0.00
  4       61      3.0       3.0       0.0     0.00
  5       64      4.0       4.0       0.0     0.00
  6       70      5.0       5.0       0.0     0.00
  7       72      6.0       6.0       0.0     0.00
  8       78      7.0       7.0       0.0     0.00
  9       82      8.0       8.0       0.0     0.00
  10      85      9.0       9.0       0.0     0.00
  11      91      10.0      10.0      0.0     0.00
  sum of d^2 = 0.00
  rho = 1 - 6.sum(d^2) / n(n^2-1) = 1 - 6(0) / 10(99) = 1.0000
  Spearman works on ranks, so it detects any MONOTONIC relationship,
  not only a straight-line one, and is unaffected by outliers.

==================================================================
EXPERIMENT 13: Simple linear regression
==================================================================

  Model: y = b0 + b1.x
  b1 (slope)     = sum((x-mx)(y-my)) / sum((x-mx)^2)
                 = 355.00 / 82.50 = 4.3030
  b0 (intercept) = my - b1.mx = 71.00 - 4.3030 x 6.50 = 43.0303

  FITTED LINE:  y = 43.0303 + 4.3030x
  Meaning: each extra hour of study is associated with about 4.30 more marks.
  The intercept 43.03 is the predicted score at zero hours -- treat
  it cautiously, since x = 0 lies outside the observed range.

  RESIDUALS
  x       observed y    fitted y      residual    residual^2
  2       52            51.636        0.364       0.1322
  3       55            55.939        -0.939      0.8825
  4       61            60.242        0.758       0.5739
  5       64            64.545        -0.545      0.2975
  6       70            68.848        1.152       1.3260
  7       72            73.152        -1.152      1.3260
  8       78            77.455        0.545       0.2975
  9       82            81.758        0.242       0.0588
  10      85            86.061        -1.061      1.1249
  11      91            90.364        0.636       0.4050

  ANALYSIS OF VARIANCE
  Source        SS            df      MS            F
  ----------------------------------------------------------
  Regression    1527.5758     1       1527.5758     1902.2642
  Residual      6.4242        8       0.8030
  Total         1534.0000     9

  F = 1902.2642 on (1, 8) df,  p = 8.425e-11
  R^2 = SS_regression / SS_total = 1527.5758 / 1534.0000 = 0.9958
  So 99.58% of the variation in scores is explained by hours.
  Note that R^2 = r^2: 0.9979^2 = 0.9958  -- true for SIMPLE
  regression only, not for multiple regression.

  TESTING THE SLOPE   H0: b1 = 0   H1: b1 != 0
  standard error of b1 = sqrt(MS_res / sum(x-mx)^2) = 0.0987
  t = b1 / se = 4.3030 / 0.0987 = 43.6150 on 8 df
  two-tailed p = 8.425e-11
  p < 0.05, so reject H0 -- the slope is significantly different from 0.
  For simple regression t^2 = F: 43.6150^2 = 1902.2642 = 1902.2642

  PREDICTION
    x = 7.5   -> predicted score 75.30
    x = 12    -> predicted score 94.67   <- EXTRAPOLATION, unsafe

  CORRELATION IS NOT CAUSATION. This data cannot show that studying
  causes higher scores -- only that they move together. A confounder
  (prior ability, motivation) could drive both.

Experiments 8 to 11 — Random variables and the probability distributions

1. Question

Units 2 and 3. 8, simulate a discrete and a continuous random variable; 9, find the expectation and variance of a probability distribution; 10, work with the binomial and Poisson distributions; 11, with the normal and exponential distributions.

2. Aim

Generate random variables, and compute probabilities from the four distributions the syllabus names.

3. Steps

1. Experiment 8: a die rolled 1000 times.
2. 1000 draws from a normal distribution.
3. Experiment 9: expectation and variance.
4. Experiment 10: the binomial.
5. The Poisson.
6. Experiment 11: the normal.
7. The exponential.

4. Programme

"""Course 4 Lab, experiments 8-11: random variables, expectation and variance,
and the discrete and continuous probability distributions.

Excel equivalents are named against each section, since the prescribed lab is
a spreadsheet lab (see ../excel-walkthroughs.md).
"""

import random
import statlib as S

random.seed(42)          # reproducible output

# Step 1: Experiment 8: a die rolled 1000 times
print("=" * 66)
print("EXPERIMENT 8: Simulating discrete and continuous random variables")
print("=" * 66)

print("\nDISCRETE -- rolling a fair die 1000 times")
rolls = [random.randint(1, 6) for _ in range(1000)]
print(f"{'Face':<8}{'Observed':<12}{'Expected':<12}Frequency")
for face in range(1, 7):
    observed = rolls.count(face)
    print(f"{face:<8}{observed:<12}{1000/6:<12.1f}{'#' * (observed // 10)}")
print(f"  Excel: =RANDBETWEEN(1,6), then COUNTIF to tally")

# Step 2: 1000 draws from a normal distribution
print("\nCONTINUOUS -- 1000 draws from Normal(mean 100, sd 15)")
sample = [random.gauss(100, 15) for _ in range(1000)]
sample_mean = sum(sample) / len(sample)
sample_var = sum((x - sample_mean) ** 2 for x in sample) / (len(sample) - 1)
print(f"  sample mean = {sample_mean:.2f}   (population mean 100)")
print(f"  sample sd   = {sample_var ** 0.5:.2f}   (population sd 15)")
print("  Excel: =NORM.INV(RAND(), 100, 15)")

# Step 3: Experiment 9: expectation and variance
print("\n" + "=" * 66)
print("EXPERIMENT 9: Expectation and variance from a probability distribution")
print("=" * 66)

# A discrete random variable given by its PMF.
values = [0, 1, 2, 3, 4]
probs = [0.10, 0.25, 0.30, 0.25, 0.10]

print(f"\n{'x':<8}{'P(x)':<10}{'x.P(x)':<12}{'x^2.P(x)':<12}")
print("-" * 42)
expectation = 0.0
second_moment = 0.0
for x, p in zip(values, probs):
    expectation += x * p
    second_moment += x * x * p
    print(f"{x:<8}{p:<10.2f}{x * p:<12.3f}{x * x * p:<12.3f}")
print("-" * 42)
print(f"{'Sum':<8}{sum(probs):<10.2f}{expectation:<12.3f}{second_moment:<12.3f}")

variance = second_moment - expectation ** 2
print(f"\n  E(X)   = sum of x.P(x)          = {expectation:.3f}")
print(f"  E(X^2) = sum of x^2.P(x)        = {second_moment:.3f}")
print(f"  Var(X) = E(X^2) - [E(X)]^2      = {second_moment:.3f} - "
      f"{expectation:.3f}^2 = {variance:.3f}")
print(f"  SD(X)  = sqrt(Var(X))           = {variance ** 0.5:.3f}")
print("\n  The shortcut Var(X) = E(X^2) - [E(X)]^2 is quicker in an exam than")
print("  computing E[(X - mu)^2] term by term, and gives the same answer.")
print("  Check that the probabilities sum to 1 first -- if they do not, the")
print("  question is misread or mistyped.")

# Step 4: Experiment 10: the binomial
print("\n" + "=" * 66)
print("EXPERIMENT 10: Discrete distributions -- Binomial and Poisson")
print("=" * 66)

n, p = 10, 0.3
print(f"\nBINOMIAL(n={n}, p={p})")
print("  Use when: a fixed number of independent trials, each success or")
print("  failure, with a constant probability of success.")
print(f"\n{'k':<6}{'P(X=k)':<12}{'P(X<=k)':<12}Distribution")
for k in range(n + 1):
    pmf = S.binomial_pmf(k, n, p)
    cdf = S.binomial_cdf(k, n, p)
    print(f"{k:<6}{pmf:<12.5f}{cdf:<12.5f}{'#' * int(pmf * 150)}")
print(f"\n  Mean     = n.p        = {n} x {p} = {n * p:.2f}")
print(f"  Variance = n.p.(1-p)  = {n} x {p} x {1-p} = {n * p * (1 - p):.2f}")
print(f"  Excel: =BINOM.DIST(k, {n}, {p}, FALSE) for the PMF, TRUE for the CDF")

lam = 3
# Step 5: The Poisson
print(f"\nPOISSON(lambda={lam})")
print("  Use when: counting events in a fixed interval of time or space,")
print("  occurring independently at a constant average rate.")
print(f"\n{'k':<6}{'P(X=k)':<12}{'P(X<=k)':<12}Distribution")
for k in range(11):
    pmf = S.poisson_pmf(k, lam)
    cdf = S.poisson_cdf(k, lam)
    print(f"{k:<6}{pmf:<12.5f}{cdf:<12.5f}{'#' * int(pmf * 150)}")
print(f"\n  Mean = Variance = lambda = {lam}")
print("  That mean and variance are equal is the signature of a Poisson.")
print(f"  Excel: =POISSON.DIST(k, {lam}, FALSE)")

print("\n  POISSON AS A LIMIT OF THE BINOMIAL (n large, p small, np = lambda):")
print(f"  {'k':<6}{'Binomial(1000, 0.003)':<24}{'Poisson(3)':<14}difference")
for k in range(6):
    b = S.binomial_pmf(k, 1000, 0.003)
    po = S.poisson_pmf(k, 3)
    print(f"  {k:<6}{b:<24.6f}{po:<14.6f}{abs(b - po):.6f}")

# Step 6: Experiment 11: the normal
print("\n" + "=" * 66)
print("EXPERIMENT 11: Continuous distributions -- Normal and Exponential")
print("=" * 66)

mu, sigma = 100, 15
print(f"\nNORMAL(mu={mu}, sigma={sigma})  -- IQ scores, a standard example")
print("\n  THE EMPIRICAL RULE (68-95-99.7)")
for k in (1, 2, 3):
    low, high = mu - k * sigma, mu + k * sigma
    prob = S.normal_cdf(high, mu, sigma) - S.normal_cdf(low, mu, sigma)
    print(f"    within {k} sd  [{low:6.1f}, {high:6.1f}] -> {prob * 100:6.2f}%")

print("\n  SPECIFIC PROBABILITIES")
print(f"    P(X <= 115)      = {S.normal_cdf(115, mu, sigma):.4f}"
      f"   Excel =NORM.DIST(115,100,15,TRUE)")
print(f"    P(X > 130)       = {1 - S.normal_cdf(130, mu, sigma):.4f}"
      f"   Excel =1-NORM.DIST(130,100,15,TRUE)")
print(f"    P(85 <= X <= 115) = "
      f"{S.normal_cdf(115, mu, sigma) - S.normal_cdf(85, mu, sigma):.4f}")

print("\n  STANDARDISING -- z = (x - mu) / sigma")
for x in (85, 100, 115, 130):
    z = (x - mu) / sigma
    print(f"    x = {x:3d}  ->  z = {z:+.2f}  ->  P(X <= x) = "
          f"{S.normal_cdf(z):.4f}")

print("\n  PERCENTILES (inverse -- Excel's NORM.INV)")
for pct in (0.90, 0.95, 0.99):
    print(f"    {pct * 100:.0f}th percentile = {S.normal_ppf(pct, mu, sigma):.2f}")

rate = 0.5
# Step 7: The exponential
print(f"\nEXPONENTIAL(lambda={rate})  -- waiting time until the next event")
print(f"  Mean = 1/lambda = {1 / rate:.2f}, Variance = 1/lambda^2 = "
      f"{1 / rate ** 2:.2f}")
print(f"\n{'x':<8}{'PDF f(x)':<14}{'CDF P(X<=x)':<16}P(X>x)")
for x in (0, 1, 2, 3, 4, 5):
    print(f"{x:<8}{S.exponential_pdf(x, rate):<14.5f}"
          f"{S.exponential_cdf(x, rate):<16.5f}"
          f"{1 - S.exponential_cdf(x, rate):.5f}")
print(f"\n  Excel: =EXPON.DIST(x, {rate}, TRUE) for the CDF")
print("\n  MEMORYLESSNESS -- P(X > s+t | X > s) = P(X > t)")
s, t = 2, 3
joint = 1 - S.exponential_cdf(s + t, rate)
given = 1 - S.exponential_cdf(s, rate)
print(f"    P(X > {s+t} | X > {s}) = {joint:.5f} / {given:.5f} = {joint / given:.5f}")
print(f"    P(X > {t})           = {1 - S.exponential_cdf(t, rate):.5f}")
print("    Equal -- having waited 2 minutes tells you nothing about the next 3.")
print("    The exponential is the only continuous distribution with this property.")

5. Execution and Results

==================================================================
EXPERIMENT 8: Simulating discrete and continuous random variables
==================================================================

DISCRETE -- rolling a fair die 1000 times
Face    Observed    Expected    Frequency
1       166         166.7       ################
2       172         166.7       #################
3       162         166.7       ################
4       165         166.7       ################
5       160         166.7       ################
6       175         166.7       #################
  Excel: =RANDBETWEEN(1,6), then COUNTIF to tally

CONTINUOUS -- 1000 draws from Normal(mean 100, sd 15)
  sample mean = 100.12   (population mean 100)
  sample sd   = 14.92   (population sd 15)
  Excel: =NORM.INV(RAND(), 100, 15)

==================================================================
EXPERIMENT 9: Expectation and variance from a probability distribution
==================================================================

x       P(x)      x.P(x)      x^2.P(x)
------------------------------------------
0       0.10      0.000       0.000
1       0.25      0.250       0.250
2       0.30      0.600       1.200
3       0.25      0.750       2.250
4       0.10      0.400       1.600
------------------------------------------
Sum     1.00      2.000       5.300

  E(X)   = sum of x.P(x)          = 2.000
  E(X^2) = sum of x^2.P(x)        = 5.300
  Var(X) = E(X^2) - [E(X)]^2      = 5.300 - 2.000^2 = 1.300
  SD(X)  = sqrt(Var(X))           = 1.140

  The shortcut Var(X) = E(X^2) - [E(X)]^2 is quicker in an exam than
  computing E[(X - mu)^2] term by term, and gives the same answer.
  Check that the probabilities sum to 1 first -- if they do not, the
  question is misread or mistyped.

==================================================================
EXPERIMENT 10: Discrete distributions -- Binomial and Poisson
==================================================================

BINOMIAL(n=10, p=0.3)
  Use when: a fixed number of independent trials, each success or
  failure, with a constant probability of success.

k     P(X=k)      P(X<=k)     Distribution
0     0.02825     0.02825     ####
1     0.12106     0.14931     ##################
2     0.23347     0.38278     ###################################
3     0.26683     0.64961     ########################################
4     0.20012     0.84973     ##############################
5     0.10292     0.95265     ###############
6     0.03676     0.98941     #####
7     0.00900     0.99841     #
8     0.00145     0.99986
9     0.00014     0.99999
10    0.00001     1.00000

  Mean     = n.p        = 10 x 0.3 = 3.00
  Variance = n.p.(1-p)  = 10 x 0.3 x 0.7 = 2.10
  Excel: =BINOM.DIST(k, 10, 0.3, FALSE) for the PMF, TRUE for the CDF

POISSON(lambda=3)
  Use when: counting events in a fixed interval of time or space,
  occurring independently at a constant average rate.

k     P(X=k)      P(X<=k)     Distribution
0     0.04979     0.04979     #######
1     0.14936     0.19915     ######################
2     0.22404     0.42319     #################################
3     0.22404     0.64723     #################################
4     0.16803     0.81526     #########################
5     0.10082     0.91608     ###############
6     0.05041     0.96649     #######
7     0.02160     0.98810     ###
8     0.00810     0.99620     #
9     0.00270     0.99890
10    0.00081     0.99971

  Mean = Variance = lambda = 3
  That mean and variance are equal is the signature of a Poisson.
  Excel: =POISSON.DIST(k, 3, FALSE)

  POISSON AS A LIMIT OF THE BINOMIAL (n large, p small, np = lambda):
  k     Binomial(1000, 0.003)   Poisson(3)    difference
  0     0.049563                0.049787      0.000224
  1     0.149137                0.149361      0.000225
  2     0.224154                0.224042      0.000112
  3     0.224379                0.224042      0.000337
  4     0.168284                0.168031      0.000253
  5     0.100869                0.100819      0.000050

==================================================================
EXPERIMENT 11: Continuous distributions -- Normal and Exponential
==================================================================

NORMAL(mu=100, sigma=15)  -- IQ scores, a standard example

  THE EMPIRICAL RULE (68-95-99.7)
    within 1 sd  [  85.0,  115.0] ->  68.27%
    within 2 sd  [  70.0,  130.0] ->  95.45%
    within 3 sd  [  55.0,  145.0] ->  99.73%

  SPECIFIC PROBABILITIES
    P(X <= 115)      = 0.8413   Excel =NORM.DIST(115,100,15,TRUE)
    P(X > 130)       = 0.0228   Excel =1-NORM.DIST(130,100,15,TRUE)
    P(85 <= X <= 115) = 0.6827

  STANDARDISING -- z = (x - mu) / sigma
    x =  85  ->  z = -1.00  ->  P(X <= x) = 0.1587
    x = 100  ->  z = +0.00  ->  P(X <= x) = 0.5000
    x = 115  ->  z = +1.00  ->  P(X <= x) = 0.8413
    x = 130  ->  z = +2.00  ->  P(X <= x) = 0.9772

  PERCENTILES (inverse -- Excel's NORM.INV)
    90th percentile = 119.22
    95th percentile = 124.67
    99th percentile = 134.90

EXPONENTIAL(lambda=0.5)  -- waiting time until the next event
  Mean = 1/lambda = 2.00, Variance = 1/lambda^2 = 4.00

x       PDF f(x)      CDF P(X<=x)     P(X>x)
0       0.50000       0.00000         1.00000
1       0.30327       0.39347         0.60653
2       0.18394       0.63212         0.36788
3       0.11157       0.77687         0.22313
4       0.06767       0.86466         0.13534
5       0.04104       0.91792         0.08208

  Excel: =EXPON.DIST(x, 0.5, TRUE) for the CDF

  MEMORYLESSNESS -- P(X > s+t | X > s) = P(X > t)
    P(X > 5 | X > 2) = 0.08208 / 0.36788 = 0.22313
    P(X > 3)           = 0.22313
    Equal -- having waited 2 minutes tells you nothing about the next 3.
    The exponential is the only continuous distribution with this property.

Experiments 14 and 15 — Confidence intervals and hypothesis tests

1. Question

Unit 5. 14, from a sample of heights, find confidence intervals for the mean; 15, carry out a one-sample z-test, a two-sample t-test, a chi-square test of independence and an F-test for two variances.

2. Aim

Estimate a mean with a stated confidence, and test hypotheses with the test the situation calls for.

3. Steps

1. Experiment 14: the sample.
2. Confidence intervals, with t.
3. Experiment 15: the steps of every test.
4. One-sample z-test.
5. Two-sample t-test.
6. Chi-square test of independence.
7. F-test for two variances.
8. Type I and Type II errors, and power.

4. Programme

"""Course 4 Lab, experiments 14-15: confidence intervals and the four
hypothesis tests named in Unit 5 -- z, t, chi-square and F.
"""

import statlib as S

# Step 1: Experiment 14: the sample
print("=" * 68)
print("EXPERIMENT 14: Estimation and confidence intervals")
print("=" * 68)

sample = [68, 72, 75, 71, 69, 74, 73, 70, 76, 72,
          71, 73, 69, 75, 74, 72, 70, 73, 71, 74]
n = len(sample)
mean = sum(sample) / n
var = sum((x - mean) ** 2 for x in sample) / (n - 1)
sd = var ** 0.5
se = sd / n ** 0.5

print(f"\n  Sample of n = {n} heights (cm)")
print(f"  mean = {mean:.3f}   sample sd = {sd:.3f}   standard error = "
      f"sd/sqrt(n) = {se:.4f}")
print("\n  POINT ESTIMATE vs INTERVAL ESTIMATE")
print(f"    A point estimate is the single number {mean:.2f}. It is almost")
print("    certainly not exactly right. An interval estimate admits that.")

# Step 2: Confidence intervals, with t
print("\n  CONFIDENCE INTERVALS -- population sd unknown, so use t")
for level, alpha in ((0.90, 0.10), (0.95, 0.05), (0.99, 0.01)):
    # Critical t by bisection on the CDF.
    low, high = 0.0, 100.0
    for _ in range(200):
        mid = (low + high) / 2
        if S.t_cdf(mid, n - 1) < 1 - alpha / 2:
            low = mid
        else:
            high = mid
    t_crit = (low + high) / 2
    margin = t_crit * se
    print(f"    {level * 100:.0f}%: t({alpha/2:.3f}, {n-1}) = {t_crit:.4f}, "
          f"margin = {margin:.4f}  ->  ({mean - margin:.3f}, {mean + margin:.3f})")

print("\n  WHAT 95% CONFIDENCE ACTUALLY MEANS")
print("    If we repeated this sampling many times and built an interval each")
print("    time, about 95% of those intervals would contain the true")
print("    population mean. It does NOT mean there is a 95% probability that")
print("    the true mean lies in THIS interval -- the true mean is a fixed")
print("    number, not a random one. Saying otherwise loses marks.")
print("\n    Notice the interval widens as confidence rises: more certainty")
print("    costs precision. It narrows as n grows, in proportion to 1/sqrt(n)")
print("    -- to halve the width you need four times the data.")

# Step 3: Experiment 15: the steps of every test
print("\n" + "=" * 68)
print("EXPERIMENT 15: Hypothesis testing")
print("=" * 68)

print("""
  THE PROCEDURE, every time:
    1. State H0 (no effect) and H1 (the claim being tested)
    2. Choose the significance level alpha, usually 0.05
    3. Compute the test statistic
    4. Find the p-value, or compare against the critical value
    5. Decide: p < alpha means reject H0
    6. State the conclusion in the words of the original problem
""")

# Step 4: One-sample z-test
print("-" * 68)
print("1. ONE-SAMPLE z-TEST -- population sd KNOWN, large sample")
print("-" * 68)
mu0, sigma_known, n_z = 70, 3.0, 40
xbar = 71.2
z = (xbar - mu0) / (sigma_known / n_z ** 0.5)
p_two = 2 * (1 - S.normal_cdf(abs(z)))
print(f"""
  A machine should fill packets to a mean of {mu0} g, with a known population
  sd of {sigma_known} g. A sample of {n_z} packets averages {xbar} g.

  H0: mu = {mu0}        H1: mu != {mu0}       alpha = 0.05
  z = (xbar - mu0) / (sigma / sqrt(n))
    = ({xbar} - {mu0}) / ({sigma_known} / sqrt({n_z}))
    = {xbar - mu0:.2f} / {sigma_known / n_z ** 0.5:.4f} = {z:.4f}
  critical value: +/- 1.96      p-value = {p_two:.4f}
  Decision: {'reject H0' if p_two < 0.05 else 'fail to reject H0'}""")
print(f"  Conclusion: the mean fill weight {'differs' if p_two < 0.05 else 'does not differ'}"
      f" significantly from {mu0} g.")
print("  Excel: =2*(1-NORM.S.DIST(ABS(z),TRUE))")

# Step 5: Two-sample t-test
print("\n" + "-" * 68)
print("2. TWO-SAMPLE t-TEST -- population sd UNKNOWN")
print("-" * 68)
group_a = [78, 82, 75, 88, 79, 84, 80, 86, 77, 83]
group_b = [72, 75, 70, 78, 74, 71, 76, 73, 69, 77]


def describe(values):
    m = sum(values) / len(values)
    v = sum((x - m) ** 2 for x in values) / (len(values) - 1)
    return m, v


m_a, v_a = describe(group_a)
m_b, v_b = describe(group_b)
n_a, n_b = len(group_a), len(group_b)

pooled_var = ((n_a - 1) * v_a + (n_b - 1) * v_b) / (n_a + n_b - 2)
se_diff = (pooled_var * (1 / n_a + 1 / n_b)) ** 0.5
t = (m_a - m_b) / se_diff
df = n_a + n_b - 2
p_t = S.t_sf_two_tailed(t, df)

print(f"""
  Do two teaching methods give different mean scores?

  Group A: n = {n_a}, mean = {m_a:.2f}, variance = {v_a:.3f}
  Group B: n = {n_b}, mean = {m_b:.2f}, variance = {v_b:.3f}

  H0: mu_A = mu_B      H1: mu_A != mu_B      alpha = 0.05

  pooled variance = [(n_A-1)s_A^2 + (n_B-1)s_B^2] / (n_A + n_B - 2)
                  = [{n_a-1} x {v_a:.3f} + {n_b-1} x {v_b:.3f}] / {df}
                  = {pooled_var:.4f}
  standard error  = sqrt(pooled x (1/n_A + 1/n_B)) = {se_diff:.4f}
  t = (mean_A - mean_B) / se = {m_a - m_b:.2f} / {se_diff:.4f} = {t:.4f}
  df = {df}      two-tailed p = {p_t:.6f}
  Decision: {'reject H0' if p_t < 0.05 else 'fail to reject H0'}""")
print(f"  Conclusion: the two methods {'do' if p_t < 0.05 else 'do not'} differ "
      f"significantly in mean score.")
print("  Excel: =T.TEST(range_A, range_B, 2, 2)")

# Step 6: Chi-square test of independence
print("\n" + "-" * 68)
print("3. CHI-SQUARE TEST OF INDEPENDENCE")
print("-" * 68)
observed = [[30, 70], [45, 55], [25, 75]]
row_labels = ["North", "South", "East"]
col_labels = ["Premium", "Standard"]

row_sums = [sum(row) for row in observed]
col_sums = [sum(observed[i][j] for i in range(len(observed)))
            for j in range(len(observed[0]))]
total = sum(row_sums)

print(f"\n  Is region independent of purchase type?")
print(f"  H0: they are independent    H1: they are associated    alpha = 0.05")
print(f"\n  {'':<10}{'Observed':<24}{'Expected':<24}")
print(f"  {'Region':<10}" + "".join(f"{c:>11}" for c in col_labels)
      + "  " + "".join(f"{c:>11}" for c in col_labels))

chi2 = 0.0
for i, row in enumerate(observed):
    expected_row = []
    for j, obs in enumerate(row):
        exp = row_sums[i] * col_sums[j] / total
        expected_row.append(exp)
        chi2 += (obs - exp) ** 2 / exp
    print(f"  {row_labels[i]:<10}" + "".join(f"{v:>11}" for v in row)
          + "  " + "".join(f"{v:>11.2f}" for v in expected_row))

df_chi = (len(observed) - 1) * (len(observed[0]) - 1)
p_chi = S.chi2_sf(chi2, df_chi)
print(f"""
  Expected = (row total x column total) / grand total
  chi-square = sum (O - E)^2 / E = {chi2:.4f}
  df = (rows - 1)(columns - 1) = ({len(observed)}-1)({len(observed[0])}-1) = {df_chi}
  critical value at 0.05 with {df_chi} df = 5.991      p = {p_chi:.4f}
  Decision: {'reject H0' if p_chi < 0.05 else 'fail to reject H0'}""")
print(f"  Conclusion: region and purchase type are "
      f"{'associated' if p_chi < 0.05 else 'independent'}.")
print("\n  ASSUMPTION: every expected frequency should be at least 5.")
print(f"  Smallest expected here = "
      f"{min(row_sums[i] * col_sums[j] / total for i in range(len(observed)) for j in range(len(col_sums))):.2f} -- satisfied.")
print("  Excel: =CHISQ.TEST(observed_range, expected_range)")

# Step 7: F-test for two variances
print("\n" + "-" * 68)
print("4. F-TEST FOR EQUALITY OF TWO VARIANCES")
print("-" * 68)
# Convention: put the larger variance on top, giving a right-tailed test.
if v_a >= v_b:
    f_stat, df1, df2, top, bottom = v_a / v_b, n_a - 1, n_b - 1, "A", "B"
else:
    f_stat, df1, df2, top, bottom = v_b / v_a, n_b - 1, n_a - 1, "B", "A"
p_f = 2 * S.f_sf(f_stat, df1, df2)      # doubled for a two-tailed test

print(f"""
  Do the two groups have equal variances? (This is the assumption the pooled
  t-test above relies on, so it is worth checking.)

  H0: sigma_A^2 = sigma_B^2      H1: they differ      alpha = 0.05

  F = larger variance / smaller variance = s_{top}^2 / s_{bottom}^2
    = {max(v_a, v_b):.4f} / {min(v_a, v_b):.4f} = {f_stat:.4f}
  df = ({df1}, {df2})      two-tailed p = {p_f:.4f}
  Decision: {'reject H0' if p_f < 0.05 else 'fail to reject H0'}""")
print(f"  Conclusion: the variances {'differ' if p_f < 0.05 else 'are not significantly different'}"
      f", so the pooled t-test above {'was not appropriate' if p_f < 0.05 else 'was appropriate'}.")
print("  Excel: =F.TEST(range_A, range_B)")

# Step 8: Type I and Type II errors, and power
print("\n" + "=" * 68)
print("TYPE I and TYPE II ERRORS, and POWER")
print("=" * 68)
print("""
                        | H0 is actually TRUE    | H0 is actually FALSE
  ----------------------|------------------------|----------------------
  We REJECT H0          | Type I error (alpha)   | Correct  (power)
  We FAIL TO REJECT H0  | Correct                | Type II error (beta)

  alpha = P(Type I)  -- convicting an innocent person. We choose this, 0.05.
  beta  = P(Type II) -- acquitting a guilty one. Follows from the design.
  Power = 1 - beta   -- the chance of detecting a real effect.

  Lowering alpha to 0.01 makes a Type I error rarer but a Type II error more
  likely. The only way to reduce both at once is a larger sample.
""")

# Power of the z-test above against a true mean of 71.5.
mu_true, alpha = 71.5, 0.05
crit = 1.959964
se_z = sigma_known / n_z ** 0.5
shift = (mu_true - mu0) / se_z
power = (1 - S.normal_cdf(crit - shift)) + S.normal_cdf(-crit - shift)
print(f"  Worked example -- power of test 1 if the true mean were {mu_true}:")
print(f"    shift = (mu_true - mu0)/se = ({mu_true} - {mu0})/{se_z:.4f} = {shift:.4f}")
print(f"    power = {power:.4f}, so beta = {1 - power:.4f}")
print(f"    A {power * 100:.0f}% chance of detecting a real shift of "
      f"{mu_true - mu0} g with n = {n_z}.")

5. Execution and Results

====================================================================
EXPERIMENT 14: Estimation and confidence intervals
====================================================================

  Sample of n = 20 heights (cm)
  mean = 72.100   sample sd = 2.222   standard error = sd/sqrt(n) = 0.4968

  POINT ESTIMATE vs INTERVAL ESTIMATE
    A point estimate is the single number 72.10. It is almost
    certainly not exactly right. An interval estimate admits that.

  CONFIDENCE INTERVALS -- population sd unknown, so use t
    90%: t(0.050, 19) = 1.7291, margin = 0.8591  ->  (71.241, 72.959)
    95%: t(0.025, 19) = 2.0930, margin = 1.0399  ->  (71.060, 73.140)
    99%: t(0.005, 19) = 2.8609, margin = 1.4214  ->  (70.679, 73.521)

  WHAT 95% CONFIDENCE ACTUALLY MEANS
    If we repeated this sampling many times and built an interval each
    time, about 95% of those intervals would contain the true
    population mean. It does NOT mean there is a 95% probability that
    the true mean lies in THIS interval -- the true mean is a fixed
    number, not a random one. Saying otherwise loses marks.

    Notice the interval widens as confidence rises: more certainty
    costs precision. It narrows as n grows, in proportion to 1/sqrt(n)
    -- to halve the width you need four times the data.

====================================================================
EXPERIMENT 15: Hypothesis testing
====================================================================

  THE PROCEDURE, every time:
    1. State H0 (no effect) and H1 (the claim being tested)
    2. Choose the significance level alpha, usually 0.05
    3. Compute the test statistic
    4. Find the p-value, or compare against the critical value
    5. Decide: p < alpha means reject H0
    6. State the conclusion in the words of the original problem

--------------------------------------------------------------------
1. ONE-SAMPLE z-TEST -- population sd KNOWN, large sample
--------------------------------------------------------------------

  A machine should fill packets to a mean of 70 g, with a known population
  sd of 3.0 g. A sample of 40 packets averages 71.2 g.

  H0: mu = 70        H1: mu != 70       alpha = 0.05
  z = (xbar - mu0) / (sigma / sqrt(n))
    = (71.2 - 70) / (3.0 / sqrt(40))
    = 1.20 / 0.4743 = 2.5298
  critical value: +/- 1.96      p-value = 0.0114
  Decision: reject H0
  Conclusion: the mean fill weight differs significantly from 70 g.
  Excel: =2*(1-NORM.S.DIST(ABS(z),TRUE))

--------------------------------------------------------------------
2. TWO-SAMPLE t-TEST -- population sd UNKNOWN
--------------------------------------------------------------------

  Do two teaching methods give different mean scores?

  Group A: n = 10, mean = 81.20, variance = 17.067
  Group B: n = 10, mean = 73.50, variance = 9.167

  H0: mu_A = mu_B      H1: mu_A != mu_B      alpha = 0.05

  pooled variance = [(n_A-1)s_A^2 + (n_B-1)s_B^2] / (n_A + n_B - 2)
                  = [9 x 17.067 + 9 x 9.167] / 18
                  = 13.1167
  standard error  = sqrt(pooled x (1/n_A + 1/n_B)) = 1.6197
  t = (mean_A - mean_B) / se = 7.70 / 1.6197 = 4.7541
  df = 18      two-tailed p = 0.000159
  Decision: reject H0
  Conclusion: the two methods do differ significantly in mean score.
  Excel: =T.TEST(range_A, range_B, 2, 2)

--------------------------------------------------------------------
3. CHI-SQUARE TEST OF INDEPENDENCE
--------------------------------------------------------------------

  Is region independent of purchase type?
  H0: they are independent    H1: they are associated    alpha = 0.05

            Observed                Expected
  Region        Premium   Standard      Premium   Standard
  North              30         70        33.33      66.67
  South              45         55        33.33      66.67
  East               25         75        33.33      66.67

  Expected = (row total x column total) / grand total
  chi-square = sum (O - E)^2 / E = 9.7500
  df = (rows - 1)(columns - 1) = (3-1)(2-1) = 2
  critical value at 0.05 with 2 df = 5.991      p = 0.0076
  Decision: reject H0
  Conclusion: region and purchase type are associated.

  ASSUMPTION: every expected frequency should be at least 5.
  Smallest expected here = 33.33 -- satisfied.
  Excel: =CHISQ.TEST(observed_range, expected_range)

--------------------------------------------------------------------
4. F-TEST FOR EQUALITY OF TWO VARIANCES
--------------------------------------------------------------------

  Do the two groups have equal variances? (This is the assumption the pooled
  t-test above relies on, so it is worth checking.)

  H0: sigma_A^2 = sigma_B^2      H1: they differ      alpha = 0.05

  F = larger variance / smaller variance = s_A^2 / s_B^2
    = 17.0667 / 9.1667 = 1.8618
  df = (9, 9)      two-tailed p = 0.3682
  Decision: fail to reject H0
  Conclusion: the variances are not significantly different, so the pooled t-test above was appropriate.
  Excel: =F.TEST(range_A, range_B)

====================================================================
TYPE I and TYPE II ERRORS, and POWER
====================================================================

                        | H0 is actually TRUE    | H0 is actually FALSE
  ----------------------|------------------------|----------------------
  We REJECT H0          | Type I error (alpha)   | Correct  (power)
  We FAIL TO REJECT H0  | Correct                | Type II error (beta)

  alpha = P(Type I)  -- convicting an innocent person. We choose this, 0.05.
  beta  = P(Type II) -- acquitting a guilty one. Follows from the design.
  Power = 1 - beta   -- the chance of detecting a real effect.

  Lowering alpha to 0.01 makes a Type I error rarer but a Type II error more
  likely. The only way to reduce both at once is a larger sample.

  Worked example -- power of test 1 if the true mean were 71.5:
    shift = (mu_true - mu0)/se = (71.5 - 70)/0.4743 = 3.1623
    power = 0.8854, so beta = 0.1146
    A 89% chance of detecting a real shift of 1.5 g with n = 40.

PSPP

PSPP is the free SPSS alternative named in the syllabus, and may be what your lab has installed.

Enter variable definitions in Variable View first, then data in Data View — the opposite order to a spreadsheet, and the usual source of confusion.

Lab exam tips

1. Label everything. Chart titles, axis labels, legends. Marks are given for a readable output, not just a correct number.

2. Show the formula, not only the result. Examiners often ask you to widen a column or press Ctrl+` to reveal formulas.

3. Interpret in a text box next to each result. "r = 0.87, a strong positive linear relationship; this does not establish causation" earns more than the number alone.

4. State H₀ and H₁ in the sheet for every test, then the statistic, then the p-value, then the decision, then a conclusion in words.

5. Check the assumptions and say you did — expected frequencies ≥ 5 for chi-square, approximate normality for t-tests.

6. Round sensibly. Two or three decimals. Copying fifteen digits from the cell looks careless.

7. Save frequently under a filename with your roll number.

What the practical record should contain

For each experiment, the five parts set out above: 1. Question, the task as set; 2. Aim, in one line; 3. Steps, the method, with the Excel functions; 4. Programme, the program (or the sheet's formulas); 5. Execution and Results, what it gave, and the result in words.

Written-out instructions

This part of the lab is a written procedure rather than a program.

Each program, on its own page

The same experiments, one page each, so a program can be reached by what it does rather than by its number.