#!/usr/bin/env python3
"""Create instructor/reference answers from the shipped CSVs.
This is a reference analysis, not a grader of a learner's own work.
Run: python tools/analyze.py
"""
from __future__ import annotations
from collections import Counter
import csv
import json
import math
from pathlib import Path
import platform
import statistics
import numpy as np
import scipy
from scipy import stats
try:
from .generate_data import SHIPMENT_SEED, TRIAL_SEED, clean_shipment_rows
except ImportError:
from generate_data import SHIPMENT_SEED, TRIAL_SEED, clean_shipment_rows
ROOT = Path(__file__).resolve().parents[1]
def finite_values(values) -> list[float]:
result = [float(value) for value in values]
if not result or not all(math.isfinite(value) for value in result):
raise ValueError("Need nonempty finite numeric values; handle missingness explicitly.")
return result
def describe(values) -> dict:
x = finite_values(values)
q1, q3 = np.quantile(x, [0.25, 0.75], method="linear")
return {"n": len(x), "sum": sum(x), "mean": statistics.mean(x),
"median": statistics.median(x),
"sample_sd": statistics.stdev(x) if len(x) >= 2 else None,
"minimum": min(x), "q1_linear": float(q1), "q3_linear": float(q3),
"maximum": max(x), "iqr_linear": float(q3 - q1)}
def mean_t_interval(values, confidence=0.95) -> dict:
x = finite_values(values)
if len(x) < 2 or not 0 < confidence < 1:
raise ValueError("Need at least two values and confidence strictly between 0 and 1.")
mean = statistics.mean(x)
se = statistics.stdev(x) / math.sqrt(len(x))
df = len(x) - 1
critical = float(stats.t.ppf((1 + confidence) / 2, df))
margin = critical * se
return {"confidence_level": confidence, "mean": mean, "n": len(x),
"df": df, "standard_error": se, "critical_t": critical,
"margin": margin, "low": mean - margin, "high": mean + margin}
def welch_difference(first, second, confidence=0.95) -> dict:
"""Two-sided Welch test and matched CI for first minus second."""
a, b = finite_values(first), finite_values(second)
if min(len(a), len(b)) < 2 or not 0 < confidence < 1:
raise ValueError("Need two or more observations per group and a valid confidence level.")
first_term = statistics.variance(a) / len(a)
second_term = statistics.variance(b) / len(b)
variance = first_term + second_term
if variance <= 0:
raise ValueError("Welch inference requires nonzero estimated variance.")
se = math.sqrt(variance)
df = variance ** 2 / (first_term ** 2 / (len(a) - 1)
+ second_term ** 2 / (len(b) - 1))
difference = statistics.mean(a) - statistics.mean(b)
t = difference / se
pvalue = float(2 * stats.t.sf(abs(t), df))
critical = float(stats.t.ppf((1 + confidence) / 2, df))
return {"first_n": len(a), "second_n": len(b), "difference": difference,
"standard_error": se, "t_statistic": t, "df": df,
"two_sided_p": pvalue, "confidence_level": confidence,
"ci_low": difference - critical * se,
"ci_high": difference + critical * se,
"null_difference": 0.0, "alternative": "two-sided"}
def linear_description(x_values, y_values) -> dict:
"""OLS line and Pearson r: descriptive summaries, no inference or causality."""
x, y = finite_values(x_values), finite_values(y_values)
if len(x) != len(y) or len(x) < 3:
raise ValueError("Need equally long sequences with at least three observations.")
xbar, ybar = statistics.mean(x), statistics.mean(y)
ssx = sum((v - xbar) ** 2 for v in x)
ssy = sum((v - ybar) ** 2 for v in y)
if ssx == 0 or ssy == 0:
raise ValueError("Both variables must vary to calculate Pearson r and R-squared.")
cross = sum((a - xbar) * (b - ybar) for a, b in zip(x, y))
slope = cross / ssx
intercept = ybar - slope * xbar
r = cross / math.sqrt(ssx * ssy)
residuals = [b - (intercept + slope * a) for a, b in zip(x, y)]
sse = sum(value ** 2 for value in residuals)
return {"n": len(x), "pearson_r": r, "slope_minutes_per_km": slope,
"intercept_minutes": intercept, "r_squared": 1 - sse / ssy,
"rmse_minutes_denominator_n": math.sqrt(sse / len(x)),
"observed_distance_min_km": min(x), "observed_distance_max_km": max(x),
"fitted_minutes_at_20_km": intercept + slope * 20,
"fitted_minutes_at_40_km": intercept + slope * 40,
"interpretation": "Association within this fictional observational snapshot, not a causal effect.",
"checks_before_inference": ["Inspect scatterplot for approximate linearity and influential observations.",
"Inspect residuals for changing spread and systematic patterns.",
"Shipments share week conditions; independent errors are not established.",
"Service mix and unmeasured operational factors may confound relationships.",
"Do not extrapolate beyond observed distances or infer real-world performance."]}
def group_summary(rows) -> dict:
total = len(rows)
late = sum(row["late"] for row in rows)
ratings = [row["satisfaction"] for row in rows if row["satisfaction"] is not None]
return {"n": total, "actual_minutes": describe(row["actual_minutes"] for row in rows),
"distance_km": describe(row["distance_km"] for row in rows),
"late_count": late, "late_denominator": total, "late_fraction": late / total,
"priority_count": sum(row["service"] == "priority" for row in rows),
"satisfaction_respondents": len(ratings), "satisfaction_missing": total - len(ratings),
"satisfaction_missing_fraction": (total - len(ratings)) / total,
"respondent_mean_rating": statistics.mean(ratings) if ratings else None,
"rating_counts_among_respondents": {str(k): ratings.count(k) for k in range(1, 6)}}
def read_csv(path: Path) -> list[dict[str, str]]:
with path.open(newline="", encoding="utf-8") as handle:
return list(csv.DictReader(handle))
def build_results(root: Path = ROOT) -> dict:
raw = read_csv(root / "data/shipments_raw.csv")
shipments, cleaning = clean_shipment_rows(raw)
clean_csv, clean_audit = clean_shipment_rows(read_csv(root / "data/shipments_clean.csv"))
if shipments != clean_csv:
raise ValueError("Shipped clean file does not match the documented raw cleaning process.")
trial = read_csv(root / "data/checklist_trial.csv")
if len({row["person_id"] for row in trial}) != len(trial):
raise ValueError("Trial person IDs must be unique.")
if {row["group"] for row in trial} != {"control", "checklist"}:
raise ValueError("Unexpected trial group.")
control = [float(row["duration_minutes"]) for row in trial if row["group"] == "control"]
checklist = [float(row["duration_minutes"]) for row in trial if row["group"] == "checklist"]
if len(control) != 20 or len(checklist) != 20 or min(control + checklist) <= 0:
raise ValueError("Expected 20 positive durations per randomized group.")
welch = welch_difference(checklist, control)
welch["direction"] = "checklist minus control, minutes; negative means shorter"
welch["absolute_difference_minutes"] = welch["difference"]
welch["absolute_magnitude_minutes"] = abs(welch["difference"])
welch["relative_difference_fraction_of_control"] = welch["difference"] / statistics.mean(control)
welch["relative_difference_percent_of_control"] = 100 * welch["relative_difference_fraction_of_control"]
welch["relative_ci_note"] = "No confidence interval for the relative change is computed; the Welch interval is in minutes."
control_ci = mean_t_interval(control)
# These independently implemented library calls verify our explicit formulas.
scipy_welch = stats.ttest_ind(checklist, control, equal_var=False, alternative="two-sided")
library_ci = scipy_welch.confidence_interval(confidence_level=0.95)
if not all(math.isclose(a, b, abs_tol=1e-10) for a, b in [
(welch["t_statistic"], float(scipy_welch.statistic)),
(welch["df"], float(scipy_welch.df)), (welch["two_sided_p"], float(scipy_welch.pvalue)),
(welch["ci_low"], float(library_ci.low)), (welch["ci_high"], float(library_ci.high))]):
raise ArithmeticError("Manual Welch formulas differ from SciPy.")
worked = describe([8, 9, 10, 11, 12, 30])
worked["values"] = [8, 9, 10, 11, 12, 30]
worked["sample_variance"] = statistics.variance(worked["values"])
worked["population_sd_descriptive_only"] = statistics.pstdev(worked["values"])
worked["tukey_lower_fence"] = worked["q1_linear"] - 1.5 * worked["iqr_linear"]
worked["tukey_upper_fence"] = worked["q3_linear"] + 1.5 * worked["iqr_linear"]
worked["beyond_fences"] = [x for x in worked["values"] if x < worked["tukey_lower_fence"] or x > worked["tukey_upper_fence"]]
return {"metadata": {"fixture_name": "Northstar Dispatch: entirely fictional teaching data",
"shipment_seed": SHIPMENT_SEED, "trial_seed": TRIAL_SEED,
"python": platform.python_version(), "numpy": np.__version__, "scipy": scipy.__version__,
"sample_sd_denominator": "n-1", "quartile_method": "NumPy linear: index (n-1)*p, linearly interpolated",
"fractions": "0–1 unless explicitly labelled percent", "results_unrounded": True},
"audits": {"raw_cleaning": cleaning,
"clean_file": {"rows": len(shipments), "unique_shipment_ids": len({r["shipment_id"] for r in shipments}),
"weeks": sorted({r["week"] for r in shipments}),
"rows_by_week": dict(sorted(Counter(str(r["week"]) for r in shipments).items())),
"additional_cleaning_required": any(clean_audit[k] for k in ["formatting_cells_normalized", "late_flags_recomputed", "exact_duplicate_rows_removed", "invalid_rows_quarantined"])},
"trial_file": {"rows": len(trial), "unique_person_ids": len({r["person_id"] for r in trial}),
"group_counts": dict(sorted(Counter(r["group"] for r in trial).items())),
"missing_duration_count": 0}},
"worked_list": worked,
"shipments": {"overall": group_summary(shipments),
"by_hub": {key: group_summary([r for r in shipments if r["hub"] == key]) for key in ["East", "West"]},
"by_service": {key: group_summary([r for r in shipments if r["service"] == key]) for key in ["priority", "standard"]},
"by_hub_service": {f"{hub}_{service}": group_summary([r for r in shipments if r["hub"] == hub and r["service"] == service]) for hub in ["East", "West"] for service in ["priority", "standard"]},
"by_week": {str(week): group_summary([r for r in shipments if r["week"] == week]) for week in range(1, 9)},
"by_late_status": {label: group_summary([r for r in shipments if r["late"] == late]) for late, label in [(0, "on_time"), (1, "late")]},
"observational_distance_time": linear_description([r["distance_km"] for r in shipments], [r["actual_minutes"] for r in shipments]),
"inference_warning": "Descriptive snapshot, not a probability sample. Shared weeks and different route/service mixes preclude naive independent-row causal hub tests.",
"satisfaction_warning": "Ratings are ordinal and response is selectively missing. The respondent mean is a coding-based summary of respondents, not all customers. Missing ratings are never zero or imputed."},
"trial": {"control": describe(control), "checklist": describe(checklist),
"control_mean_t_interval_95": control_ci,
"welch_checklist_minus_control": welch,
"design": "40 independent simulated people from one fictional teaching process, one task each, balanced random assignment. Not a representative sample of any real workforce.",
"generating_model": "Duration = 24 + independent Normal(0,4.2) person error - 3 if assigned checklist; rounded to 0.01 minute. Each simulated person's two potential outcomes differ by -3 before rounding; only the assigned outcome is observed. No rerolls.",
"inference_scope": "Randomization illustrates causal identification within this simulated experiment. t methods assume independent errors and approximately normal within-group outcomes; equal population variances are not required by Welch. No claim about real people or shipment delivery performance.",
"p_value_warning": "The p-value is not the probability that the null is true, that chance caused the result, or that a future rollout will succeed.",
"ci_warning": "95% describes the long-run coverage of the method under its assumptions, not a 95% probability for a fixed unknown parameter after this interval is observed. Non-significance would not establish equivalence."},
"probability_toy": {"note": "Separate stipulated 1,000-shipment toy, not estimates from the CSV.",
"total": 1000, "late": 100, "on_time": 900,
"flagged_late": 80, "unflagged_late": 20,
"flagged_on_time": 90, "unflagged_on_time": 810,
"base_late_probability": 0.10,
"p_flag_given_late": 0.80, "p_flag_given_on_time": 0.10,
"p_flag": 170 / 1000, "p_late_given_flag": 80 / 170,
"p_late_or_flag": 190 / 1000}}
def main() -> None:
result = build_results()
output = ROOT / "outputs/reference_results.json"
output.parent.mkdir(exist_ok=True)
output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n", encoding="utf-8")
print(f"Wrote {output.relative_to(ROOT)} (instructor/reference answers).")
print("Trial difference (checklist - control):", result["trial"]["welch_checklist_minus_control"])
if __name__ == "__main__":
main()
