精度評価実験の分析を複数被験者対応

config/experiment.yaml

@@ -1,5 +1,5 @@
 experiment:
-  experiment_id: phase1_test_001
+  experiment_id: test
   location: laboratory
   operator: student
   notes: "Phase 2 pedestrian detection and tracking test"

src/accuracy_report.py

@@ -14,6 +14,26 @@ Threshold: ``theta_th = |bias| + k * SD`` of the frame errors for conditions wit
 ``|true_yaw| <= --pool-max-angle`` (pitch 0), per distance. ``--basis worst`` (default) takes the
 worst single condition, ``pooled`` pools all of them. k=3 and the pooling range are provisional
 choices to be agreed on with the supervisor.
+
+Multiple subjects
+-----------------
+Give every subject's marker file(s); the subject name is taken from, in this order, a
+``LABEL=path`` prefix on ``--trials``, the ``subject`` column written by ``accuracy_session plan
+--subject``, or the marker file name::
+
+    python -m src.accuracy_report --trials A=data/analysis/accuracy/trials_a.csv ^
+        B=data/analysis/accuracy/trials_b.csv C=data/analysis/accuracy/trials_c.csv ^
+        --face-csv "data/logs/face/*_face_frames.csv"
+
+The thresholds, regression and ``accuracy_conditions.csv`` pool the frames of all subjects (so the
+SD contains the between-subject bias differences, which is how individual bias is averaged out).
+With two or more subjects the report additionally contains per-subject thresholds
+(``accuracy_subject_thresholds.csv`` / ``accuracy_subject_conditions.csv``), a leave-one-subject-out
+check (``accuracy_loso.csv``), the +/- direction asymmetry of the gain (``accuracy_asymmetry*.csv``),
+a between-subject comparison of the gain (``accuracy_subject_differences.csv``) and per-subject plots.
+The face logs of different subjects must not overlap in time (they are joined by timestamp).
+``--exclude-subjects``, ``--subjects``, ``--alias OLD=NEW``, ``--flip-subject NAME`` and
+``--exclude-seq 12 S01:7`` (seq of one subject only) adjust what is analysed.
 """
 
 from __future__ import annotations
@@ -23,23 +43,33 @@ import bisect
 import csv
 import glob
 import math
+import random
 import statistics
 import sys
 from collections import Counter
 from dataclasses import dataclass
 from datetime import datetime
 from pathlib import Path
-from typing import Sequence
+from typing import Mapping, Sequence
 
 NAN = float("nan")
 TRIAL_COLUMNS = ["session_id", "seq", "condition_id", "distance_m", "true_yaw_deg", "pitch_deg", "repeat", "subject",
                  "status", "track_id", "baseline_n", "baseline_mean", "baseline_sd", "hold_n", "hold_mean", "hold_sd",
                  "delta_deg", "error_deg", "detection_rate"]
-CONDITION_COLUMNS = ["distance_m", "pitch_deg", "true_yaw_deg", "n_trials", "n_frames", "bias_deg", "sd_trials_deg",
+CONDITION_BASE_COLUMNS = ["distance_m", "pitch_deg", "true_yaw_deg", "n_trials", "n_frames", "bias_deg",
-                     "sd_within_deg", "sd_total_deg", "rmse_frame_deg", "p95_abs_error_deg", "mean_detection_rate",
+                          "sd_trials_deg", "sd_within_deg", "sd_total_deg", "rmse_frame_deg", "p95_abs_error_deg",
-                     "exceed_rate"]
+                          "mean_detection_rate", "exceed_rate"]
-THRESHOLD_COLUMNS = ["distance_m", "n_frames", "bias_deg", "sd_deg", "theta_pooled_deg", "theta_worst_deg",
+CONDITION_COLUMNS = CONDITION_BASE_COLUMNS + ["n_subjects", "sd_subjects_deg"]
+THRESHOLD_BASE_COLUMNS = ["distance_m", "n_frames", "bias_deg", "sd_deg", "theta_pooled_deg", "theta_worst_deg",
                           "theta_th_deg", "slope", "intercept_deg", "r2", "min_detectable_deg"]
+THRESHOLD_COLUMNS = THRESHOLD_BASE_COLUMNS + ["n_subjects"]
+SUBJECT_CONDITION_COLUMNS = ["subject"] + CONDITION_BASE_COLUMNS
+SUBJECT_THRESHOLD_COLUMNS = ["subject"] + THRESHOLD_BASE_COLUMNS
+LOSO_COLUMNS = ["held_out", "distance_m", "n_train_subjects", "theta_train_deg", "false_alarm_rate",
+                "min_detectable_deg"]
+ASYMMETRY_CELL_COLUMNS = ["subject", "distance_m", "abs_angle_deg", "gain_pos", "gain_neg", "gain_diff"]
+ASYMMETRY_COLUMNS = ["subject", "n_cells", "gain_pos", "gain_neg", "gain_diff", "p_signflip"]
+DIFFERENCE_COLUMNS = ["distance_m", "n_subjects", "kruskal_h", "p_permutation"]
 
 
 @dataclass(frozen=True)
@@ -57,6 +87,7 @@ class Options:
     basis: str = "worst"
     detect_ratio: float = 0.95
     pitch_ref: float = 0.0
+    flip_subjects: tuple[str, ...] = ()
 
 
 @dataclass
@@ -113,6 +144,14 @@ def expand_paths(patterns: Sequence[str]) -> list[Path]:
     return paths
 
 
+def split_label(token: str) -> tuple[str, str]:
+    """``LABEL=path`` -> (LABEL, path); anything else -> ("", token)."""
+    label, sep, rest = token.partition("=")
+    if sep and label and rest and not any(c in label for c in "\\/:*?[]"):
+        return label, rest
+    return "", token
+
+
 def load_face_series(paths: Sequence[Path], default_fps: float = 30.0) -> FaceSeries:
     rows: list[tuple[float, int, bool, float | None, float | None, float]] = []
     for path in paths:
@@ -139,14 +178,38 @@ def load_face_series(paths: Sequence[Path], default_fps: float = 30.0) -> FaceSe
                       [r[4] for r in rows], [r[5] for r in rows])
 
 
-def load_trials(paths: Sequence[Path], exclude_seqs: Sequence[int] = ()) -> list[dict[str, str]]:
+def _parse_exclusions(tokens: Sequence[int | str]) -> tuple[set[int], set[tuple[str, int]]]:
-    """Last log row per (session, seq); only status ``ok`` survives."""
+    """``12`` excludes seq 12 of every subject, ``S01:12`` only that of subject S01."""
-    last: dict[tuple[str, int], dict[str, str]] = {}
+    general: set[int] = set()
+    specific: set[tuple[str, int]] = set()
+    for token in tokens:
+        text = str(token)
+        if ":" in text:
+            name, _, number = text.rpartition(":")
+            specific.add((name, int(number)))
+        else:
+            general.add(int(text))
+    return general, specific
+
+
+def load_trials(paths: Sequence[Path], exclude_seqs: Sequence[int | str] = (),
+                labels: Mapping[Path, str] | None = None,
+                aliases: Mapping[str, str] | None = None) -> list[dict[str, str]]:
+    """Last log row per (subject, session, seq); only status ``ok`` survives.
+
+    The subject is the ``LABEL`` given for the file, else the ``subject`` column, else the file name.
+    """
+    labels, aliases = labels or {}, aliases or {}
+    last: dict[tuple[str, str, int], dict[str, str]] = {}
     for path in paths:
         with path.open(encoding="utf-8-sig", newline="") as stream:
             for row in csv.DictReader(stream):
-                last[(row.get("session_id", ""), int(row["seq"]))] = row
+                subject = labels.get(path) or (row.get("subject") or "").strip() or path.stem
-    return [row for (_, seq), row in sorted(last.items()) if row["status"] == "ok" and seq not in set(exclude_seqs)]
+                row["subject"] = aliases.get(subject, subject)
+                last[(row["subject"], row.get("session_id", ""), int(row["seq"]))] = row
+    general, specific = _parse_exclusions(exclude_seqs)
+    return [row for (subject, _, seq), row in sorted(last.items())
+            if row["status"] == "ok" and seq not in general and (subject, seq) not in specific]
 
 
 def evaluate_trial(trial: dict[str, str], series: FaceSeries, opt: Options) -> tuple[dict[str, object], list[float]]:
@@ -191,7 +254,8 @@ def evaluate_trial(trial: dict[str, str], series: FaceSeries, opt: Options) -> t
     if len(hold) < opt.min_frames:
         record["status"] = "hold_insufficient"
         return record, []
-    sign = -1.0 if opt.flip_sign else 1.0
+    flipped = opt.flip_sign or str(record["subject"]) in opt.flip_subjects
+    sign = -1.0 if flipped else 1.0
     base_mean = float(record["baseline_mean"])  # type: ignore[arg-type]
     deltas = [_wrap(sign * (y - base_mean)) for y in hold_yaw]
     fps = opt.fps or series.fps[hold[0]]
@@ -209,6 +273,17 @@ def _errors(record: dict[str, object], deltas: Sequence[float]) -> list[float]:
     return [d - float(record["true_yaw_deg"]) for d in deltas]  # type: ignore[arg-type]
 
 
+def _subject(record: dict[str, object]) -> str:
+    return str(record.get("subject") or "")
+
+
+def split_by_subject(usable: Usable) -> dict[str, Usable]:
+    groups: dict[str, Usable] = {}
+    for item in usable:
+        groups.setdefault(_subject(item[0]), []).append(item)
+    return groups
+
+
 def summarize_conditions(usable: Usable) -> list[dict[str, object]]:
     groups: dict[tuple[float, float, float], Usable] = {}
     for record, deltas in usable:
@@ -218,6 +293,9 @@ def summarize_conditions(usable: Usable) -> list[dict[str, object]]:
     for (distance, pitch, true), items in sorted(groups.items()):
         frame_errors = [e for record, deltas in items for e in _errors(record, deltas)]
         within = [float(r["hold_sd"]) ** 2 for r, _ in items if not math.isnan(float(r["hold_sd"]))]  # type: ignore[arg-type]
+        by_subject: dict[str, list[float]] = {}
+        for r, _ in items:
+            by_subject.setdefault(_subject(r), []).append(float(r["error_deg"]))  # type: ignore[arg-type]
         rows.append({
             "distance_m": distance, "pitch_deg": pitch, "true_yaw_deg": true, "n_trials": len(items),
             "n_frames": len(frame_errors), "bias_deg": _mean([float(r["error_deg"]) for r, _ in items]),  # type: ignore[arg-type]
@@ -226,7 +304,8 @@ def summarize_conditions(usable: Usable) -> list[dict[str, object]]:
             "rmse_frame_deg": math.sqrt(_mean([e * e for e in frame_errors])),
             "p95_abs_error_deg": percentile(sorted(abs(e) for e in frame_errors), 95),
             "mean_detection_rate": _mean([float(r["detection_rate"]) for r, _ in items]),  # type: ignore[arg-type]
-            "exceed_rate": NAN,
+            "exceed_rate": NAN, "n_subjects": len(by_subject),
+            "sd_subjects_deg": _sd([_mean(v) for v in by_subject.values()]),
         })
     return rows
 
@@ -264,6 +343,7 @@ def compute_thresholds(usable: Usable, opt: Options) -> dict[float, dict[str, fl
             "theta_pooled_deg": pooled, "theta_worst_deg": worst,
             "theta_th_deg": worst if opt.basis == "worst" else pooled,
             "slope": slope, "intercept_deg": intercept, "r2": r2, "min_detectable_deg": NAN,
+            "n_subjects": len({_subject(r) for r, _ in pooled_sel}),
         }
     return result
 
@@ -288,6 +368,162 @@ def add_exceed_rates(conditions: list[dict[str, object]], usable: Usable, thresh
                 break
 
 
+# ---------------------------------------------------------------------------------------------
+# Multi-subject analysis
+# ---------------------------------------------------------------------------------------------
+def subject_analysis(usable: Usable, opt: Options) -> dict[str, dict[str, object]]:
+    """Conditions and thresholds computed from each subject alone."""
+    result: dict[str, dict[str, object]] = {}
+    for subject, items in sorted(split_by_subject(usable).items()):
+        conditions = summarize_conditions(items)
+        thresholds = compute_thresholds(items, opt)
+        add_exceed_rates(conditions, items, thresholds, opt)
+        result[subject] = {"conditions": conditions, "thresholds": thresholds}
+    return result
+
+
+def leave_one_subject_out(usable: Usable, opt: Options) -> list[dict[str, object]]:
+    """Threshold made from the other subjects, applied to the held-out one (false alarm at 0 deg, min angle)."""
+    groups = split_by_subject(usable)
+    if len(groups) < 2:
+        return []
+    rows: list[dict[str, object]] = []
+    for held, test in sorted(groups.items()):
+        train = [item for subject, items in groups.items() if subject != held for item in items]
+        trained = compute_thresholds(train, opt)
+        carried = {d: {"theta_th_deg": t["theta_th_deg"], "min_detectable_deg": NAN} for d, t in trained.items()}
+        conditions = summarize_conditions(test)
+        add_exceed_rates(conditions, test, carried, opt)
+        for distance, info in sorted(carried.items()):
+            zero = [c for c in conditions if c["distance_m"] == distance and c["pitch_deg"] == opt.pitch_ref
+                    and float(c["true_yaw_deg"]) == 0.0]  # type: ignore[arg-type]
+            if not any(c["distance_m"] == distance for c in conditions):
+                continue
+            rows.append({"held_out": held, "distance_m": distance,
+                         "n_train_subjects": len({_subject(r) for r, _ in train}),
+                         "theta_train_deg": info["theta_th_deg"],
+                         "false_alarm_rate": zero[0]["exceed_rate"] if zero else NAN,
+                         "min_detectable_deg": info["min_detectable_deg"]})
+    return rows
+
+
+def gain_records(usable: Usable, opt: Options) -> list[dict[str, object]]:
+    """Per-trial gain = (estimated change toward the set direction) / |set angle|; 1 is ideal."""
+    rows: list[dict[str, object]] = []
+    for record, _ in usable:
+        true = float(record["true_yaw_deg"])  # type: ignore[arg-type]
+        if float(record["pitch_deg"]) != opt.pitch_ref or true == 0.0:  # type: ignore[arg-type]
+            continue
+        rows.append({"subject": _subject(record), "distance_m": float(record["distance_m"]),  # type: ignore[arg-type]
+                     "abs_angle_deg": abs(true), "side": "pos" if true > 0 else "neg",
+                     "gain": float(record["delta_deg"]) * math.copysign(1.0, true) / abs(true)})  # type: ignore[arg-type]
+    return rows
+
+
+def asymmetry_cells(gains: Sequence[dict[str, object]]) -> list[dict[str, object]]:
+    cells: dict[tuple[str, float, float], dict[str, list[float]]] = {}
+    for g in gains:
+        key = (str(g["subject"]), float(g["distance_m"]), float(g["abs_angle_deg"]))  # type: ignore[arg-type]
+        cells.setdefault(key, {"pos": [], "neg": []})[str(g["side"])].append(float(g["gain"]))  # type: ignore[arg-type]
+    rows: list[dict[str, object]] = []
+    for (subject, distance, angle), sides in sorted(cells.items()):
+        if sides["pos"] and sides["neg"]:
+            pos, neg = _mean(sides["pos"]), _mean(sides["neg"])
+            rows.append({"subject": subject, "distance_m": distance, "abs_angle_deg": angle,
+                         "gain_pos": pos, "gain_neg": neg, "gain_diff": pos - neg})
+    return rows
+
+
+def _signflip_p(diffs: Sequence[float], *, seed: int = 0, exact_max: int = 16, samples: int = 20000) -> float:
+    """Two-sided paired sign-flip permutation test of mean(diffs) = 0 (exact for small n)."""
+    n = len(diffs)
+    if n < 2 or all(d == 0 for d in diffs):
+        return NAN
+    observed, eps = abs(sum(diffs)), 1e-12
+    if n <= exact_max:
+        hits = sum(abs(sum(d if (mask >> i) & 1 else -d for i, d in enumerate(diffs))) >= observed - eps
+                   for mask in range(1 << n))
+        return hits / (1 << n)
+    rng = random.Random(seed)
+    hits = sum(abs(sum(d if rng.random() < 0.5 else -d for d in diffs)) >= observed - eps for _ in range(samples))
+    return (hits + 1) / (samples + 1)
+
+
+def asymmetry_summary(cells: Sequence[dict[str, object]]) -> list[dict[str, object]]:
+    """Per-subject and overall +/- gain difference. The p value treats cells as independent (reference only)."""
+    rows: list[dict[str, object]] = []
+    subjects = sorted({str(c["subject"]) for c in cells})
+    for name, selected in [(s, [c for c in cells if c["subject"] == s]) for s in subjects] + [("ALL", list(cells))]:
+        if not selected:
+            continue
+        diffs = [float(c["gain_diff"]) for c in selected]  # type: ignore[arg-type]
+        rows.append({"subject": name, "n_cells": len(selected),
+                     "gain_pos": _mean([float(c["gain_pos"]) for c in selected]),  # type: ignore[arg-type]
+                     "gain_neg": _mean([float(c["gain_neg"]) for c in selected]),  # type: ignore[arg-type]
+                     "gain_diff": _mean(diffs), "p_signflip": _signflip_p(diffs)})
+    return rows
+
+
+def _ranks(values: Sequence[float]) -> list[float]:
+    order = sorted(range(len(values)), key=lambda i: values[i])
+    ranks = [0.0] * len(values)
+    i = 0
+    while i < len(order):
+        j = i
+        while j + 1 < len(order) and values[order[j + 1]] == values[order[i]]:
+            j += 1
+        for k in range(i, j + 1):
+            ranks[order[k]] = (i + j) / 2.0 + 1.0
+        i = j + 1
+    return ranks
+
+
+def _kruskal_h(rank_groups: Sequence[Sequence[float]]) -> float:
+    total = sum(len(g) for g in rank_groups)
+    return 12.0 / (total * (total + 1)) * sum(sum(g) ** 2 / len(g) for g in rank_groups) - 3.0 * (total + 1)
+
+
+def kruskal_permutation(groups: Sequence[Sequence[float]], *, seed: int = 0, samples: int = 5000) -> tuple[float, float]:
+    """Kruskal-Wallis H and its permutation p value (no scipy needed, valid with ties)."""
+    flat = [v for g in groups for v in g]
+    ranks = _ranks(flat)
+    sizes = [len(g) for g in groups]
+
+    def split(values: Sequence[float]) -> list[list[float]]:
+        out, start = [], 0
+        for size in sizes:
+            out.append(list(values[start:start + size]))
+            start += size
+        return out
+
+    observed = _kruskal_h(split(ranks))
+    rng, hits = random.Random(seed), 0
+    shuffled = list(ranks)
+    for _ in range(samples):
+        rng.shuffle(shuffled)
+        hits += _kruskal_h(split(shuffled)) >= observed - 1e-12
+    return observed, (hits + 1) / (samples + 1)
+
+
+def subject_differences(gains: Sequence[dict[str, object]], min_abs_angle: float = 10.0) -> list[dict[str, object]]:
+    """Per distance: does the gain (|set angle| >= 10 deg) differ between subjects?"""
+    rows: list[dict[str, object]] = []
+    for distance in sorted({float(g["distance_m"]) for g in gains}):  # type: ignore[arg-type]
+        by_subject: dict[str, list[float]] = {}
+        for g in gains:
+            if float(g["distance_m"]) == distance and float(g["abs_angle_deg"]) >= min_abs_angle:  # type: ignore[arg-type]
+                by_subject.setdefault(str(g["subject"]), []).append(float(g["gain"]))  # type: ignore[arg-type]
+        groups = [v for v in by_subject.values() if len(v) >= 2]
+        if len(groups) < 2:
+            continue
+        h, p = kruskal_permutation(groups)
+        rows.append({"distance_m": distance, "n_subjects": len(groups), "kruskal_h": h, "p_permutation": p})
+    return rows
+
+
+# ---------------------------------------------------------------------------------------------
+# Output
+# ---------------------------------------------------------------------------------------------
 def _fmt(value: object, digits: int = 2) -> str:
     if isinstance(value, float):
         return "-" if math.isnan(value) else f"{value:.{digits}f}"
@@ -303,18 +539,24 @@ def _write_csv(path: Path, columns: Sequence[str], rows: Sequence[dict[str, obje
 
 
 def build_markdown(records: list[dict[str, object]], conditions: list[dict[str, object]],
-                   thresholds: dict[float, dict[str, float]], opt: Options) -> str:
+                   thresholds: dict[float, dict[str, float]], opt: Options,
+                   multi: dict[str, object] | None = None) -> str:
     lines = ["# 振り向き角度推定の精度実験 結果", ""]
     status = Counter(str(r["status"]) for r in records)
     lines += [f"- 記録された有効試行: {len(records)}、解析に使用: {status['ok']}", "- 除外: " + (
         ", ".join(f"{k}={v}" for k, v in status.items() if k != "ok") or "なし")]
-    bad = [str(r["seq"]) for r in records if r["status"] != "ok"]
+    bad = [str(r["seq"]) if not multi else f"{_subject(r)}:{r['seq']}" for r in records if r["status"] != "ok"]
     if bad:
         lines.append(f"- 除外された seq: {', '.join(bad)}")
-    lines += [f"- 解析設定: signal={opt.signal}, flip_sign={opt.flip_sign}, settle={opt.settle:g}s, "
+    if multi:
+        counts: dict[str, Counter[str]] = multi["counts"]  # type: ignore[assignment]
+        lines.append("- 被験者: " + ", ".join(f"{s} ({c['ok']}試行)" for s, c in sorted(counts.items())))
+    flip = f", flip_subjects={','.join(opt.flip_subjects)}" if opt.flip_subjects else ""
+    lines += [f"- 解析設定: signal={opt.signal}, flip_sign={opt.flip_sign}{flip}, settle={opt.settle:g}s, "
               f"baseline_trim={opt.baseline_trim:g}s, k={opt.k:g}, pool_max_angle={opt.pool_max_angle:g}°, "
-              f"basis={opt.basis}", "", "誤差の定義: error = (保持区間の平均yaw − 基準yaw) − 真値角(度)。", ""]
+              f"basis={opt.basis}", "", "誤差の定義: error = (保持区間の平均yaw − 基準yaw) − 真値角(度)。"
-    lines += ["## しきい値 (距離別・pitch 基準条件)", "",
+              + ("閾値・回帰は全被験者のフレームをプールして算出 (SDには被験者間のバイアス差が含まれる)。" if multi else ""), ""]
+    lines += ["## しきい値 (距離別・pitch 基準条件" + (", 全被験者プール)" if multi else ")"), "",
               "| 距離[m] | frames | bias | SD | θ_pooled | θ_worst | **θ_th** | 傾き | 切片 | R² | 判定可能な最小角 |",
               "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"]
     warn = []
@@ -328,12 +570,17 @@ def build_markdown(records: list[dict[str, object]], conditions: list[dict[str,
     if thresholds:
         lines += ["", f"全距離で共通に使う保守的な値: θ_th = {max(t['theta_th_deg'] for t in thresholds.values()):.2f}°"]
     lines += [f"- {w}" for w in warn]
-    lines += ["", "## 条件別の結果", "",
+    if multi:
-              "| 距離[m] | pitch | 真値[°] | 試行 | frames | bias | SD(試行間) | SD(試行内) | SD(全体) | RMSE | 95%誤差 | 検出率 | θ超過率 |",
+        lines += _multi_markdown(multi)
-              "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"]
+    extra_head = " 被験者 | SD(被験者間) |" if multi else ""
+    extra_rule = "---:|---:|" if multi else ""
+    lines += ["", "## 条件別の結果" + (" (全被験者プール)" if multi else ""), "",
+              "| 距離[m] | pitch | 真値[°] | 試行 |" + extra_head + " frames | bias | SD(試行間) | SD(試行内) | SD(全体) | RMSE | 95%誤差 | 検出率 | θ超過率 |",
+              "|---:|---:|---:|---:|" + extra_rule + "---:|---:|---:|---:|---:|---:|---:|---:|---:|"]
     for c in conditions:
+        extra = f" {c['n_subjects']} | {_fmt(c['sd_subjects_deg'])} |" if multi else ""
         lines.append(f"| {float(c['distance_m']):g} | {float(c['pitch_deg']):+g} | {float(c['true_yaw_deg']):+g} | "  # type: ignore[arg-type]
-                     f"{c['n_trials']} | {c['n_frames']} | {_fmt(c['bias_deg'])} | {_fmt(c['sd_trials_deg'])} | "
+                     f"{c['n_trials']} |{extra} {c['n_frames']} | {_fmt(c['bias_deg'])} | {_fmt(c['sd_trials_deg'])} | "
                      f"{_fmt(c['sd_within_deg'])} | {_fmt(c['sd_total_deg'])} | {_fmt(c['rmse_frame_deg'])} | "
                      f"{_fmt(c['p95_abs_error_deg'])} | {_fmt(c['mean_detection_rate'], 3)} | {_fmt(c['exceed_rate'], 3)} |")
     lines += ["", "注: θ超過率は保持区間のフレームのうち |Δyaw| ≥ θ_th の割合 (真値0°では誤検出率、それ以外では検出率)。",
@@ -341,6 +588,49 @@ def build_markdown(records: list[dict[str, object]], conditions: list[dict[str,
     return "\n".join(lines) + "\n"
 
 
+def _multi_markdown(multi: dict[str, object]) -> list[str]:
+    lines: list[str] = ["", "## 被験者別しきい値 (各被験者単独で算出)", "",
+                        "| 被験者 | 距離[m] | frames | bias | SD | θ_th | 傾き | R² | 判定可能な最小角 |",
+                        "|---|---:|---:|---:|---:|---:|---:|---:|---:|"]
+    per_subject: dict[str, dict[str, object]] = multi["per_subject"]  # type: ignore[assignment]
+    by_distance: dict[float, list[float]] = {}
+    warn: list[str] = []
+    for subject, info in per_subject.items():
+        for distance, t in sorted(info["thresholds"].items()):  # type: ignore[attr-defined]
+            lines.append(f"| {subject} | {distance:g} | {int(t['n_frames'])} | {_fmt(t['bias_deg'])} | {_fmt(t['sd_deg'])} | "
+                         f"**{_fmt(t['theta_th_deg'])}** | {_fmt(t['slope'], 3)} | {_fmt(t['r2'], 3)} | "
+                         f"{_fmt(t['min_detectable_deg'])} |")
+            by_distance.setdefault(distance, []).append(t["theta_th_deg"])
+            if not math.isnan(t["slope"]) and t["slope"] < 0:
+                warn.append(f"{subject} / {distance:g} m: 傾きが負です。符号定義が逆の可能性があります (--flip-subject {subject})。")
+    if by_distance:
+        lines += ["", "距離別の被験者間レンジ(θ_th): " + ", ".join(
+            f"{d:g} m: {min(v):.1f}〜{max(v):.1f}°" for d, v in sorted(by_distance.items())),
+            f"全被験者・全距離で最も保守的な値: θ_th = {max(max(v) for v in by_distance.values()):.2f}°"]
+    lines += [f"- {w}" for w in warn]
+    loso: list[dict[str, object]] = multi["loso"]  # type: ignore[assignment]
+    if loso:
+        lines += ["", "## 被験者除外 (leave-one-subject-out) 検証", "",
+                  "他の被験者だけで求めたθ_thを除外した被験者に適用した場合の、真値0°の誤検出率と判定可能な最小角。", "",
+                  "| 除外した被験者 | 距離[m] | 学習θ_th | 誤検出率(0°) | 判定可能な最小角 |", "|---|---:|---:|---:|---:|"]
+        lines += [f"| {r['held_out']} | {float(r['distance_m']):g} | {_fmt(r['theta_train_deg'])} | "  # type: ignore[arg-type]
+                  f"{_fmt(r['false_alarm_rate'], 3)} | {_fmt(r['min_detectable_deg'])} |" for r in loso]
+    asym: list[dict[str, object]] = multi["asymmetry"]  # type: ignore[assignment]
+    if asym:
+        lines += ["", "## 正負方向の非対称 (ゲイン = 設定方向への推定変化量 / 設定角, 1が理想)", "",
+                  "| 被験者 | セル数 | ゲイン(正) | ゲイン(負) | 差(正−負) | p(符号反転検定) |", "|---|---:|---:|---:|---:|---:|"]
+        lines += [f"| {r['subject']} | {r['n_cells']} | {_fmt(r['gain_pos'], 3)} | {_fmt(r['gain_neg'], 3)} | "
+                  f"{_fmt(r['gain_diff'], 3)} | {_fmt(r['p_signflip'], 3)} |" for r in asym]
+        lines += ["", "セル = (被験者×距離×|設定角|)。pはセルを独立とみなした参考値で、被験者数が少ない場合は記述統計として扱う。"]
+    diffs: list[dict[str, object]] = multi["differences"]  # type: ignore[assignment]
+    if diffs:
+        lines += ["", "## 被験者間差 (|設定角|≥10°のゲイン, Kruskal-Wallis, 並べ替え検定)", "",
+                  "| 距離[m] | 被験者数 | H | p |", "|---:|---:|---:|---:|"]
+        lines += [f"| {float(r['distance_m']):g} | {r['n_subjects']} | {_fmt(r['kruskal_h'])} | {_fmt(r['p_permutation'], 3)} |"  # type: ignore[arg-type]
+                  for r in diffs]
+    return lines
+
+
 def make_plots(usable: Usable, thresholds: dict[float, dict[str, float]], out_dir: Path, opt: Options) -> list[Path]:
     try:
         import matplotlib
@@ -413,9 +703,120 @@ def make_plots(usable: Usable, thresholds: dict[float, dict[str, float]], out_di
     return paths
 
 
+def make_subject_plots(usable: Usable, thresholds: dict[float, dict[str, float]], cells: Sequence[dict[str, object]],
+                       out_dir: Path, opt: Options) -> list[Path]:
+    """Per-subject versions of the figures (only meaningful with two or more subjects)."""
+    try:
+        import matplotlib
+        matplotlib.use("Agg", force=True)
+        import matplotlib.pyplot as plt
+    except ImportError:
+        return []
+    ref = [(r, d) for r, d in usable if float(r["pitch_deg"]) == opt.pitch_ref]  # type: ignore[arg-type]
+    subjects = sorted({_subject(r) for r, _ in ref})
+    distances = sorted({float(r["distance_m"]) for r, _ in ref})  # type: ignore[arg-type]
+    if len(subjects) < 2 or not distances:
+        return []
+    out_dir.mkdir(parents=True, exist_ok=True)
+    axis_label_size, tick_label_size = 14, 13
+    cmap = plt.get_cmap("tab10")
+    color = {s: cmap(i % 10) for i, s in enumerate(subjects)}
+    paths: list[Path] = []
+    n = len(distances)
+
+    fig, axes = plt.subplots(1, n, figsize=(5.2 * n, 5), sharex=True, sharey=True, squeeze=False)
+    lim = max(abs(float(r["true_yaw_deg"])) for r, _ in ref) * 1.15 + 1  # type: ignore[arg-type]
+    for ax, distance in zip(axes[0], distances):
+        ax.plot([-lim, lim], [-lim, lim], "k--", linewidth=1)
+        theta = thresholds.get(distance, {}).get("theta_th_deg")
+        if theta is not None:
+            ax.axhspan(-theta, theta, color="red", alpha=0.08)
+        for i, subject in enumerate(subjects):
+            jitter = (i - (len(subjects) - 1) / 2) * 0.8
+            pts = [(float(r["true_yaw_deg"]) + jitter, float(r["delta_deg"])) for r, _ in ref  # type: ignore[arg-type]
+                   if _subject(r) == subject and float(r["distance_m"]) == distance]  # type: ignore[arg-type]
+            ax.scatter([p[0] for p in pts], [p[1] for p in pts], s=18, alpha=0.8, color=color[subject], label=subject)
+        ax.set_title(f"{distance:g} m")
+        ax.set_xlabel("true yaw [deg]", fontsize=axis_label_size)
+        ax.tick_params(axis="both", labelsize=tick_label_size)
+        ax.grid(alpha=0.3)
+    axes[0][0].set_ylabel("estimated yaw change [deg]", fontsize=axis_label_size)
+    axes[0][-1].legend()
+    fig.tight_layout()
+    paths.append(out_dir / "accuracy_scatter_by_subject.png")
+    fig.savefig(paths[-1])
+    plt.close(fig)
+
+    fig, axes = plt.subplots(n, 1, figsize=(10, 3.4 * n), sharex=True, squeeze=False)
+    angles = sorted({float(r["true_yaw_deg"]) for r, _ in ref})  # type: ignore[arg-type]
+    step = min((b - a for a, b in zip(angles, angles[1:])), default=5.0)
+    width = 0.8 * step / len(subjects)
+    for ax, distance in zip(axes[:, 0], distances):
+        for i, subject in enumerate(subjects):
+            offset = (i - (len(subjects) - 1) / 2) * width
+            data, positions = [], []
+            for angle in angles:
+                values = [float(r["error_deg"]) for r, _ in ref if _subject(r) == subject  # type: ignore[arg-type]
+                          and float(r["distance_m"]) == distance and float(r["true_yaw_deg"]) == angle]  # type: ignore[arg-type]
+                if values:
+                    data.append(values)
+                    positions.append(angle + offset)
+            if data:
+                box = ax.boxplot(data, positions=positions, widths=width * 0.9, patch_artist=True, showfliers=False,
+                                 manage_ticks=False)
+                for patch in box["boxes"]:
+                    patch.set_facecolor(color[subject])
+                    patch.set_alpha(0.6)
+        theta = thresholds.get(distance, {}).get("theta_th_deg")
+        ax.axhline(0, color="k", linewidth=0.8)
+        if theta is not None:
+            ax.axhline(theta, color="r", linestyle=":")
+            ax.axhline(-theta, color="r", linestyle=":")
+        ax.set_ylabel(f"error [deg] ({distance:g} m)", fontsize=axis_label_size)
+        ax.tick_params(axis="both", labelsize=tick_label_size)
+        ax.grid(alpha=0.3)
+    axes[-1, 0].set_xticks(angles)
+    axes[-1, 0].set_xlabel("true yaw [deg]", fontsize=axis_label_size)
+    handles = [plt.Line2D([0], [0], color=color[s], lw=6, alpha=0.6) for s in subjects]
+    axes[0, 0].legend(handles, subjects, ncol=len(subjects))
+    fig.tight_layout()
+    paths.append(out_dir / "accuracy_error_by_subject.png")
+    fig.savefig(paths[-1])
+    plt.close(fig)
+
+    if cells:
+        fig, axes = plt.subplots(1, n, figsize=(5.2 * n, 4.4), sharey=True, squeeze=False)
+        for ax, distance in zip(axes[0], distances):
+            for subject in subjects:
+                rows = sorted((c for c in cells if c["subject"] == subject and c["distance_m"] == distance),
+                              key=lambda c: float(c["abs_angle_deg"]))  # type: ignore[arg-type]
+                ax.plot([c["abs_angle_deg"] for c in rows], [c["gain_pos"] for c in rows], "-o", color=color[subject],
+                        markersize=4, label=f"{subject} (+)")
+                ax.plot([c["abs_angle_deg"] for c in rows], [c["gain_neg"] for c in rows], "--s", color=color[subject],
+                        markersize=4, label=f"{subject} (-)")
+            ax.axhline(1.0, color="k", linewidth=0.8)
+            ax.set_title(f"{distance:g} m")
+            ax.set_xlabel("|true yaw| [deg]", fontsize=axis_label_size)
+            ax.tick_params(axis="both", labelsize=tick_label_size)
+            ax.grid(alpha=0.3)
+        axes[0][0].set_ylabel("gain (estimated / true)", fontsize=axis_label_size)
+        axes[0][-1].legend(fontsize=10)
+        fig.tight_layout()
+        paths.append(out_dir / "accuracy_gain_by_side.png")
+        fig.savefig(paths[-1])
+        plt.close(fig)
+    return paths
+
+
 def run_report(trial_paths: Sequence[Path], face_paths: Sequence[Path], out_dir: Path, opt: Options,
-               exclude_seqs: Sequence[int] = (), plots: bool = True) -> dict[str, object]:
+               exclude_seqs: Sequence[int | str] = (), plots: bool = True,
-    trials = load_trials(trial_paths, exclude_seqs)
+               labels: Mapping[Path, str] | None = None, include_subjects: Sequence[str] = (),
+               exclude_subjects: Sequence[str] = (), aliases: Mapping[str, str] | None = None) -> dict[str, object]:
+    trials = load_trials(trial_paths, exclude_seqs, labels, aliases)
+    if include_subjects:
+        trials = [t for t in trials if t["subject"] in set(include_subjects)]
+    if exclude_subjects:
+        trials = [t for t in trials if t["subject"] not in set(exclude_subjects)]
     series = load_face_series(face_paths)
     evaluated = [evaluate_trial(t, series, opt) for t in trials]
     records = [rec for rec, _ in evaluated]
@@ -423,25 +824,54 @@ def run_report(trial_paths: Sequence[Path], face_paths: Sequence[Path], out_dir:
     conditions = summarize_conditions(usable)
     thresholds = compute_thresholds(usable, opt)
     add_exceed_rates(conditions, usable, thresholds, opt)
+    subjects = sorted({_subject(r) for r in records})
+    multi: dict[str, object] | None = None
+    cells: list[dict[str, object]] = []
+    if len(split_by_subject(usable)) >= 2:
+        counts: dict[str, Counter[str]] = {}
+        for r in records:
+            counts.setdefault(_subject(r), Counter())[str(r["status"])] += 1
+        gains = gain_records(usable, opt)
+        cells = asymmetry_cells(gains)
+        multi = {"counts": counts, "per_subject": subject_analysis(usable, opt),
+                 "loso": leave_one_subject_out(usable, opt), "cells": cells, "asymmetry": asymmetry_summary(cells),
+                 "differences": subject_differences(gains)}
     out_dir.mkdir(parents=True, exist_ok=True)
     _write_csv(out_dir / "accuracy_trials.csv", TRIAL_COLUMNS, records)
     _write_csv(out_dir / "accuracy_conditions.csv", CONDITION_COLUMNS, conditions)
     _write_csv(out_dir / "accuracy_thresholds.csv", THRESHOLD_COLUMNS,
                [{"distance_m": d, **t} for d, t in sorted(thresholds.items())])
-    markdown = build_markdown(records, conditions, thresholds, opt)
+    if multi:
+        per: dict[str, dict[str, object]] = multi["per_subject"]  # type: ignore[assignment]
+        _write_csv(out_dir / "accuracy_subject_conditions.csv", SUBJECT_CONDITION_COLUMNS,
+                   [{"subject": s, **c} for s, info in per.items() for c in info["conditions"]])  # type: ignore[attr-defined]
+        _write_csv(out_dir / "accuracy_subject_thresholds.csv", SUBJECT_THRESHOLD_COLUMNS,
+                   [{"subject": s, "distance_m": d, **t} for s, info in per.items()
+                    for d, t in sorted(info["thresholds"].items())])  # type: ignore[attr-defined]
+        _write_csv(out_dir / "accuracy_loso.csv", LOSO_COLUMNS, multi["loso"])  # type: ignore[arg-type]
+        _write_csv(out_dir / "accuracy_asymmetry_cells.csv", ASYMMETRY_CELL_COLUMNS, cells)
+        _write_csv(out_dir / "accuracy_asymmetry.csv", ASYMMETRY_COLUMNS, multi["asymmetry"])  # type: ignore[arg-type]
+        _write_csv(out_dir / "accuracy_subject_differences.csv", DIFFERENCE_COLUMNS, multi["differences"])  # type: ignore[arg-type]
+    markdown = build_markdown(records, conditions, thresholds, opt, multi)
     (out_dir / "accuracy_report.md").write_text(markdown, encoding="utf-8")
-    figures = make_plots(usable, thresholds, out_dir / "plots", opt) if plots else []
+    figures: list[Path] = []
+    if plots:
+        figures = make_plots(usable, thresholds, out_dir / "plots", opt)
+        if multi:
+            figures += make_subject_plots(usable, thresholds, cells, out_dir / "plots", opt)
     return {"records": records, "conditions": conditions, "thresholds": thresholds, "markdown": markdown,
-            "figures": figures}
+            "figures": figures, "subjects": subjects, "multi": multi}
 
 
 def main(argv: Sequence[str] | None = None, out=print) -> int:  # type: ignore[no-untyped-def]
     parser = argparse.ArgumentParser(description="Accuracy report from trial markers and face_frames logs")
-    parser.add_argument("--trials", nargs="+", required=True, help="trials CSV file(s) from accuracy_session")
+    parser.add_argument("--trials", nargs="+", required=True,
+                        help="trials CSV file(s) from accuracy_session; prefix LABEL= to name the subject (A=trials_a.csv)")
     parser.add_argument("--face-csv", nargs="+", required=True, help="*_face_frames.csv file(s); globs allowed")
     parser.add_argument("--out", type=Path, default=Path("data/analysis/accuracy/report"))
     parser.add_argument("--signal", choices=("yaw", "raw_yaw"), default="yaw")
-    parser.add_argument("--flip-sign", action="store_true")
+    parser.add_argument("--flip-sign", action="store_true", help="flip the yaw sign of every subject")
+    parser.add_argument("--flip-subject", nargs="*", default=[], help="flip the yaw sign of these subjects only")
     parser.add_argument("--settle", type=float, default=3.0, help="seconds after the cue ignored (turn in progress)")
     parser.add_argument("--baseline-trim", type=float, default=0.5)
     parser.add_argument("--min-frames", type=int, default=10)
@@ -452,17 +882,37 @@ def main(argv: Sequence[str] | None = None, out=print) -> int:  # type: ignore[n
     parser.add_argument("--pool-max-angle", type=float, default=15.0)
     parser.add_argument("--basis", choices=("worst", "pooled"), default="worst")
     parser.add_argument("--detect-ratio", type=float, default=0.95)
-    parser.add_argument("--exclude-seq", type=int, nargs="*", default=[])
+    parser.add_argument("--exclude-seq", nargs="*", default=[],
+                        help="seq numbers to drop; SUBJECT:seq drops it for that subject only")
+    parser.add_argument("--subjects", nargs="*", default=[], help="analyse only these subjects")
+    parser.add_argument("--exclude-subjects", nargs="*", default=[], help="drop these subjects")
+    parser.add_argument("--alias", nargs="*", default=[], help="merge subject names, e.g. S01b=S01")
     parser.add_argument("--no-plots", action="store_true")
     args = parser.parse_args(argv)
-    trial_paths, face_paths = expand_paths(args.trials), expand_paths(args.face_csv)
+    labels: dict[Path, str] = {}
+    trial_paths: list[Path] = []
+    for token in args.trials:
+        label, pattern = split_label(token)
+        found = expand_paths([pattern])
+        trial_paths.extend(found)
+        if label:
+            labels.update({path: label for path in found})
+    face_paths = expand_paths(args.face_csv)
     missing = [str(p) for p in trial_paths + face_paths if not p.is_file()]
     if missing or not trial_paths or not face_paths:
         print(f"file not found: {', '.join(missing) or 'no input files'}", file=sys.stderr)
         return 2
+    try:
+        aliases = dict(item.split("=", 1) for item in args.alias)
+        _parse_exclusions(args.exclude_seq)
+    except ValueError:
+        print("invalid --alias (use OLD=NEW) or --exclude-seq (use N or SUBJECT:N)", file=sys.stderr)
+        return 2
     opt = Options(args.signal, args.flip_sign, args.settle, args.baseline_trim, args.min_frames, args.face_every,
-                  args.track_id, args.fps, args.k, args.pool_max_angle, args.basis, args.detect_ratio)
+                  args.track_id, args.fps, args.k, args.pool_max_angle, args.basis, args.detect_ratio,
-    result = run_report(trial_paths, face_paths, args.out, opt, args.exclude_seq, not args.no_plots)
+                  flip_subjects=tuple(args.flip_subject))
+    result = run_report(trial_paths, face_paths, args.out, opt, args.exclude_seq, not args.no_plots, labels,
+                        args.subjects, args.exclude_subjects, aliases)
     if not result["records"]:
         print("no usable trials (status ok) in the trials file", file=sys.stderr)
         return 1

src/accuracy_report.py.orig

@@ -0,0 +1,475 @@
+"""Accuracy-experiment report: join trial markers with ``*_face_frames.csv`` by timestamp.
+
+Usage (after ``src.accuracy_session`` recorded the markers)::
+
+    python -m src.accuracy_report --trials data/analysis/accuracy/trials_<session>.csv ^
+        --face-csv "data/logs/face/*_face_frames.csv"
+
+Per trial the baseline yaw (subject facing the camera) is averaged, the hold-phase yaw is
+compared against it, and ``error = (hold - baseline) - true_yaw``. Outputs (``--out``):
+``accuracy_trials.csv``, ``accuracy_conditions.csv``, ``accuracy_thresholds.csv``,
+``accuracy_report.md`` and ``plots/*.png``.
+
+Threshold: ``theta_th = |bias| + k * SD`` of the frame errors for conditions with
+``|true_yaw| <= --pool-max-angle`` (pitch 0), per distance. ``--basis worst`` (default) takes the
+worst single condition, ``pooled`` pools all of them. k=3 and the pooling range are provisional
+choices to be agreed on with the supervisor.
+"""
+
+from __future__ import annotations
+
+import argparse
+import bisect
+import csv
+import glob
+import math
+import statistics
+import sys
+from collections import Counter
+from dataclasses import dataclass
+from datetime import datetime
+from pathlib import Path
+from typing import Sequence
+
+NAN = float("nan")
+TRIAL_COLUMNS = ["session_id", "seq", "condition_id", "distance_m", "true_yaw_deg", "pitch_deg", "repeat", "subject",
+                 "status", "track_id", "baseline_n", "baseline_mean", "baseline_sd", "hold_n", "hold_mean", "hold_sd",
+                 "delta_deg", "error_deg", "detection_rate"]
+CONDITION_COLUMNS = ["distance_m", "pitch_deg", "true_yaw_deg", "n_trials", "n_frames", "bias_deg", "sd_trials_deg",
+                     "sd_within_deg", "sd_total_deg", "rmse_frame_deg", "p95_abs_error_deg", "mean_detection_rate",
+                     "exceed_rate"]
+THRESHOLD_COLUMNS = ["distance_m", "n_frames", "bias_deg", "sd_deg", "theta_pooled_deg", "theta_worst_deg",
+                     "theta_th_deg", "slope", "intercept_deg", "r2", "min_detectable_deg"]
+
+
+@dataclass(frozen=True)
+class Options:
+    signal: str = "yaw"
+    flip_sign: bool = False
+    settle: float = 1.5
+    baseline_trim: float = 0.5
+    min_frames: int = 10
+    face_every: int = 1
+    track_id: int | None = None
+    fps: float | None = None
+    k: float = 3.0
+    pool_max_angle: float = 15.0
+    basis: str = "worst"
+    detect_ratio: float = 0.95
+    pitch_ref: float = 0.0
+
+
+@dataclass
+class FaceSeries:
+    times: list[float]
+    track_ids: list[int]
+    pose: list[bool]
+    yaw: list[float | None]
+    raw_yaw: list[float | None]
+    fps: list[float]
+
+
+def _ts(value: object) -> float | None:
+    try:
+        return datetime.fromisoformat(str(value).strip()).timestamp()
+    except ValueError:
+        return None
+
+
+def _float(value: object) -> float | None:
+    try:
+        number = float(str(value))
+    except ValueError:
+        return None
+    return number if math.isfinite(number) else None
+
+
+def _mean(values: Sequence[float]) -> float:
+    return statistics.fmean(values) if values else NAN
+
+
+def _sd(values: Sequence[float]) -> float:
+    return statistics.stdev(values) if len(values) >= 2 else NAN
+
+
+def _wrap(angle: float) -> float:
+    return (angle + 180.0) % 360.0 - 180.0
+
+
+def percentile(sorted_values: Sequence[float], q: float) -> float:
+    if not sorted_values:
+        return NAN
+    position = (len(sorted_values) - 1) * q / 100.0
+    lower, upper = math.floor(position), math.ceil(position)
+    return float(sorted_values[lower] + (sorted_values[upper] - sorted_values[lower]) * (position - lower))
+
+
+def expand_paths(patterns: Sequence[str]) -> list[Path]:
+    """Expand globs ourselves because Windows shells do not."""
+    paths: list[Path] = []
+    for pattern in patterns:
+        matches = sorted(glob.glob(pattern)) if any(c in pattern for c in "*?[") else [pattern]
+        paths.extend(Path(m) for m in matches)
+    return paths
+
+
+def load_face_series(paths: Sequence[Path], default_fps: float = 30.0) -> FaceSeries:
+    rows: list[tuple[float, int, bool, float | None, float | None, float]] = []
+    for path in paths:
+        file_rows: list[tuple[float, int, bool, float | None, float | None]] = []
+        frames: list[tuple[float, float]] = []
+        with path.open(encoding="utf-8-sig", newline="") as stream:
+            for row in csv.DictReader(stream):
+                t, track = _ts(row.get("timestamp")), _float(row.get("track_id"))
+                if t is None or track is None:
+                    continue
+                file_rows.append((t, int(track), str(row.get("pose_estimated", "")).lower() == "true",
+                                  _float(row.get("yaw")), _float(row.get("raw_yaw"))))
+                frame = _float(row.get("frame_number"))
+                if frame is not None:
+                    frames.append((t, frame))
+        fps = default_fps
+        if len(frames) >= 2:
+            dt, df = frames[-1][0] - frames[0][0], frames[-1][1] - frames[0][1]
+            if dt > 0 and df > 0:
+                fps = df / dt
+        rows.extend((*r, fps) for r in file_rows)
+    rows.sort(key=lambda r: r[0])
+    return FaceSeries([r[0] for r in rows], [r[1] for r in rows], [r[2] for r in rows], [r[3] for r in rows],
+                      [r[4] for r in rows], [r[5] for r in rows])
+
+
+def load_trials(paths: Sequence[Path], exclude_seqs: Sequence[int] = ()) -> list[dict[str, str]]:
+    """Last log row per (session, seq); only status ``ok`` survives."""
+    last: dict[tuple[str, int], dict[str, str]] = {}
+    for path in paths:
+        with path.open(encoding="utf-8-sig", newline="") as stream:
+            for row in csv.DictReader(stream):
+                last[(row.get("session_id", ""), int(row["seq"]))] = row
+    return [row for (_, seq), row in sorted(last.items()) if row["status"] == "ok" and seq not in set(exclude_seqs)]
+
+
+def evaluate_trial(trial: dict[str, str], series: FaceSeries, opt: Options) -> tuple[dict[str, object], list[float]]:
+    """Return the trial record and the frame-level deltas (degrees, signed, hold phase)."""
+    record: dict[str, object] = {column: NAN for column in TRIAL_COLUMNS}
+    record.update({
+        "session_id": trial.get("session_id", ""), "seq": int(trial["seq"]), "condition_id": trial["condition_id"],
+        "distance_m": float(trial["distance_m"]), "true_yaw_deg": float(trial["true_yaw_deg"]),
+        "pitch_deg": float(trial["pitch_deg"]), "repeat": trial.get("repeat", ""), "subject": trial.get("subject", ""),
+        "status": "ok", "track_id": "",
+    })
+    t_base, t_cue, t_end = (_ts(trial.get(key)) for key in ("t_baseline_start", "t_turn_cue", "t_hold_end"))
+    if t_base is None or t_cue is None or t_end is None:
+        record["status"] = "bad_timestamps"
+        return record, []
+    b0, b1, h0, h1 = t_base + opt.baseline_trim, t_cue, t_cue + opt.settle, t_end
+    values = series.raw_yaw if opt.signal == "raw_yaw" else series.yaw
+
+    def span(a: float, b: float) -> tuple[int, int]:
+        return bisect.bisect_left(series.times, a), bisect.bisect_right(series.times, b)
+
+    lo, hi = span(b0, h1)
+    counts = Counter(series.track_ids[i] for i in range(lo, hi) if series.pose[i] and values[i] is not None)
+    track = opt.track_id if opt.track_id is not None else (counts.most_common(1)[0][0] if counts else None)
+    if track is None or counts.get(track, 0) == 0:
+        record["status"] = "no_track"
+        return record, []
+    record["track_id"] = track
+
+    def pick(a: float, b: float) -> list[int]:
+        s, e = span(a, b)
+        return [i for i in range(s, e) if series.track_ids[i] == track and series.pose[i] and values[i] is not None]
+
+    base, hold = pick(b0, b1), pick(h0, h1)
+    base_yaw = [float(values[i]) for i in base]  # type: ignore[arg-type]
+    hold_yaw = [float(values[i]) for i in hold]  # type: ignore[arg-type]
+    record.update({"baseline_n": len(base), "hold_n": len(hold), "baseline_mean": _mean(base_yaw),
+                   "baseline_sd": _sd(base_yaw), "hold_mean": _mean(hold_yaw), "hold_sd": _sd(hold_yaw)})
+    if len(base) < opt.min_frames:
+        record["status"] = "baseline_insufficient"
+        return record, []
+    if len(hold) < opt.min_frames:
+        record["status"] = "hold_insufficient"
+        return record, []
+    sign = -1.0 if opt.flip_sign else 1.0
+    base_mean = float(record["baseline_mean"])  # type: ignore[arg-type]
+    deltas = [_wrap(sign * (y - base_mean)) for y in hold_yaw]
+    fps = opt.fps or series.fps[hold[0]]
+    expected = max(1.0, (h1 - h0) * fps / max(1, opt.face_every))
+    record["detection_rate"] = min(1.0, len(hold) / expected)
+    record["delta_deg"] = _mean(deltas)
+    record["error_deg"] = float(record["delta_deg"]) - float(record["true_yaw_deg"])  # type: ignore[arg-type]
+    return record, deltas
+
+
+Usable = list[tuple[dict[str, object], list[float]]]
+
+
+def _errors(record: dict[str, object], deltas: Sequence[float]) -> list[float]:
+    return [d - float(record["true_yaw_deg"]) for d in deltas]  # type: ignore[arg-type]
+
+
+def summarize_conditions(usable: Usable) -> list[dict[str, object]]:
+    groups: dict[tuple[float, float, float], Usable] = {}
+    for record, deltas in usable:
+        key = (float(record["distance_m"]), float(record["pitch_deg"]), float(record["true_yaw_deg"]))  # type: ignore[arg-type]
+        groups.setdefault(key, []).append((record, deltas))
+    rows: list[dict[str, object]] = []
+    for (distance, pitch, true), items in sorted(groups.items()):
+        frame_errors = [e for record, deltas in items for e in _errors(record, deltas)]
+        within = [float(r["hold_sd"]) ** 2 for r, _ in items if not math.isnan(float(r["hold_sd"]))]  # type: ignore[arg-type]
+        rows.append({
+            "distance_m": distance, "pitch_deg": pitch, "true_yaw_deg": true, "n_trials": len(items),
+            "n_frames": len(frame_errors), "bias_deg": _mean([float(r["error_deg"]) for r, _ in items]),  # type: ignore[arg-type]
+            "sd_trials_deg": _sd([float(r["error_deg"]) for r, _ in items]),  # type: ignore[arg-type]
+            "sd_within_deg": math.sqrt(_mean(within)) if within else NAN, "sd_total_deg": _sd(frame_errors),
+            "rmse_frame_deg": math.sqrt(_mean([e * e for e in frame_errors])),
+            "p95_abs_error_deg": percentile(sorted(abs(e) for e in frame_errors), 95),
+            "mean_detection_rate": _mean([float(r["detection_rate"]) for r, _ in items]),  # type: ignore[arg-type]
+            "exceed_rate": NAN,
+        })
+    return rows
+
+
+def _line_fit(points: Sequence[tuple[float, float]]) -> tuple[float, float, float]:
+    xs, ys = [p[0] for p in points], [p[1] for p in points]
+    if len(points) < 2 or len(set(xs)) < 2:
+        return NAN, NAN, NAN
+    mx, my = statistics.fmean(xs), statistics.fmean(ys)
+    sxx = sum((x - mx) ** 2 for x in xs)
+    sxy = sum((x - mx) * (y - my) for x, y in points)
+    syy = sum((y - my) ** 2 for y in ys)
+    slope = sxy / sxx
+    return slope, my - slope * mx, (sxy * sxy / (sxx * syy) if syy > 0 else NAN)
+
+
+def compute_thresholds(usable: Usable, opt: Options) -> dict[float, dict[str, float]]:
+    result: dict[float, dict[str, float]] = {}
+    for distance in sorted({float(r["distance_m"]) for r, _ in usable}):  # type: ignore[arg-type]
+        ref = [(r, d) for r, d in usable if float(r["distance_m"]) == distance  # type: ignore[arg-type]
+               and float(r["pitch_deg"]) == opt.pitch_ref]  # type: ignore[arg-type]
+        pooled_sel = [(r, d) for r, d in ref if abs(float(r["true_yaw_deg"])) <= opt.pool_max_angle]  # type: ignore[arg-type]
+        all_errors = [e for r, d in pooled_sel for e in _errors(r, d)]
+        if len(all_errors) < 2:
+            continue
+        pooled = abs(_mean(all_errors)) + opt.k * _sd(all_errors)
+        by_condition: dict[float, list[float]] = {}
+        for r, d in pooled_sel:
+            by_condition.setdefault(float(r["true_yaw_deg"]), []).extend(_errors(r, d))  # type: ignore[arg-type]
+        worst_values = [abs(_mean(e)) + opt.k * _sd(e) for e in by_condition.values() if len(e) >= 2]
+        worst = max(worst_values) if worst_values else pooled
+        slope, intercept, r2 = _line_fit([(float(r["true_yaw_deg"]), float(r["delta_deg"])) for r, _ in ref])  # type: ignore[arg-type]
+        result[distance] = {
+            "n_frames": float(len(all_errors)), "bias_deg": _mean(all_errors), "sd_deg": _sd(all_errors),
+            "theta_pooled_deg": pooled, "theta_worst_deg": worst,
+            "theta_th_deg": worst if opt.basis == "worst" else pooled,
+            "slope": slope, "intercept_deg": intercept, "r2": r2, "min_detectable_deg": NAN,
+        }
+    return result
+
+
+def add_exceed_rates(conditions: list[dict[str, object]], usable: Usable, thresholds: dict[float, dict[str, float]],
+                     opt: Options) -> None:
+    """exceed_rate = share of hold frames with |delta| >= theta_th (false-positive rate at 0 deg)."""
+    frames: dict[tuple[float, float, float], list[float]] = {}
+    for record, deltas in usable:
+        key = (float(record["distance_m"]), float(record["pitch_deg"]), float(record["true_yaw_deg"]))  # type: ignore[arg-type]
+        frames.setdefault(key, []).extend(deltas)
+    for row in conditions:
+        theta = thresholds.get(float(row["distance_m"]), {}).get("theta_th_deg")  # type: ignore[arg-type]
+        data = frames[(float(row["distance_m"]), float(row["pitch_deg"]), float(row["true_yaw_deg"]))]  # type: ignore[arg-type]
+        row["exceed_rate"] = sum(abs(d) >= theta for d in data) / len(data) if theta is not None and data else NAN
+    for distance, info in thresholds.items():
+        ref = [c for c in conditions if c["distance_m"] == distance and c["pitch_deg"] == opt.pitch_ref
+               and float(c["true_yaw_deg"]) != 0.0 and not math.isnan(float(c["exceed_rate"]))]  # type: ignore[arg-type]
+        for angle in sorted({abs(float(c["true_yaw_deg"])) for c in ref}):  # type: ignore[arg-type]
+            if all(float(c["exceed_rate"]) >= opt.detect_ratio for c in ref if abs(float(c["true_yaw_deg"])) >= angle):  # type: ignore[arg-type]
+                info["min_detectable_deg"] = angle
+                break
+
+
+def _fmt(value: object, digits: int = 2) -> str:
+    if isinstance(value, float):
+        return "-" if math.isnan(value) else f"{value:.{digits}f}"
+    return str(value)
+
+
+def _write_csv(path: Path, columns: Sequence[str], rows: Sequence[dict[str, object]]) -> None:
+    with path.open("w", encoding="utf-8-sig", newline="") as stream:
+        writer = csv.DictWriter(stream, fieldnames=list(columns), extrasaction="ignore")
+        writer.writeheader()
+        for row in rows:
+            writer.writerow({k: (_fmt(v, 4) if isinstance(v, float) else v) for k, v in row.items()})
+
+
+def build_markdown(records: list[dict[str, object]], conditions: list[dict[str, object]],
+                   thresholds: dict[float, dict[str, float]], opt: Options) -> str:
+    lines = ["# 振り向き角度推定の精度実験 結果", ""]
+    status = Counter(str(r["status"]) for r in records)
+    lines += [f"- 記録された有効試行: {len(records)}、解析に使用: {status['ok']}", "- 除外: " + (
+        ", ".join(f"{k}={v}" for k, v in status.items() if k != "ok") or "なし")]
+    bad = [str(r["seq"]) for r in records if r["status"] != "ok"]
+    if bad:
+        lines.append(f"- 除外された seq: {', '.join(bad)}")
+    lines += [f"- 解析設定: signal={opt.signal}, flip_sign={opt.flip_sign}, settle={opt.settle:g}s, "
+              f"baseline_trim={opt.baseline_trim:g}s, k={opt.k:g}, pool_max_angle={opt.pool_max_angle:g}°, "
+              f"basis={opt.basis}", "", "誤差の定義: error = (保持区間の平均yaw − 基準yaw) − 真値角(度)。", ""]
+    lines += ["## しきい値 (距離別・pitch 基準条件)", "",
+              "| 距離[m] | frames | bias | SD | θ_pooled | θ_worst | **θ_th** | 傾き | 切片 | R² | 判定可能な最小角 |",
+              "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"]
+    warn = []
+    for distance, t in sorted(thresholds.items()):
+        lines.append(f"| {distance:g} | {int(t['n_frames'])} | {_fmt(t['bias_deg'])} | {_fmt(t['sd_deg'])} | "
+                     f"{_fmt(t['theta_pooled_deg'])} | {_fmt(t['theta_worst_deg'])} | **{_fmt(t['theta_th_deg'])}** | "
+                     f"{_fmt(t['slope'], 3)} | {_fmt(t['intercept_deg'])} | {_fmt(t['r2'], 3)} | "
+                     f"{_fmt(t['min_detectable_deg'])} |")
+        if not math.isnan(t["slope"]) and t["slope"] < 0:
+            warn.append(f"{distance:g} m: 傾きが負です。角度の符号定義が逆の可能性があります (--flip-sign)。")
+    if thresholds:
+        lines += ["", f"全距離で共通に使う保守的な値: θ_th = {max(t['theta_th_deg'] for t in thresholds.values()):.2f}°"]
+    lines += [f"- {w}" for w in warn]
+    lines += ["", "## 条件別の結果", "",
+              "| 距離[m] | pitch | 真値[°] | 試行 | frames | bias | SD(試行間) | SD(試行内) | SD(全体) | RMSE | 95%誤差 | 検出率 | θ超過率 |",
+              "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"]
+    for c in conditions:
+        lines.append(f"| {float(c['distance_m']):g} | {float(c['pitch_deg']):+g} | {float(c['true_yaw_deg']):+g} | "  # type: ignore[arg-type]
+                     f"{c['n_trials']} | {c['n_frames']} | {_fmt(c['bias_deg'])} | {_fmt(c['sd_trials_deg'])} | "
+                     f"{_fmt(c['sd_within_deg'])} | {_fmt(c['sd_total_deg'])} | {_fmt(c['rmse_frame_deg'])} | "
+                     f"{_fmt(c['p95_abs_error_deg'])} | {_fmt(c['mean_detection_rate'], 3)} | {_fmt(c['exceed_rate'], 3)} |")
+    lines += ["", "注: θ超過率は保持区間のフレームのうち |Δyaw| ≥ θ_th の割合 (真値0°では誤検出率、それ以外では検出率)。",
+              "検出率は姿勢推定に成功したフレーム数 / 期待フレーム数 (フレームレートはログから推定)。"]
+    return "\n".join(lines) + "\n"
+
+
+def make_plots(usable: Usable, thresholds: dict[float, dict[str, float]], out_dir: Path, opt: Options) -> list[Path]:
+    try:
+        import matplotlib
+        matplotlib.use("Agg", force=True)
+        import matplotlib.pyplot as plt
+    except ImportError:
+        return []
+    out_dir.mkdir(parents=True, exist_ok=True)
+    axis_label_size = 14
+    tick_label_size = 13
+    ref = [(r, d) for r, d in usable if float(r["pitch_deg"]) == opt.pitch_ref]  # type: ignore[arg-type]
+    distances = sorted({float(r["distance_m"]) for r, _ in ref})  # type: ignore[arg-type]
+    if not distances:
+        return []
+    paths: list[Path] = []
+    fig, ax = plt.subplots(figsize=(6, 6))
+    for distance in distances:
+        pts = [(float(r["true_yaw_deg"]), float(r["delta_deg"])) for r, _ in ref if float(r["distance_m"]) == distance]  # type: ignore[arg-type]
+        ax.scatter([p[0] for p in pts], [p[1] for p in pts], s=18, alpha=0.7, label=f"{distance:g} m")
+    lim = max(abs(float(r["true_yaw_deg"])) for r, _ in ref) * 1.15 + 1  # type: ignore[arg-type]
+    ax.plot([-lim, lim], [-lim, lim], "k--", linewidth=1, label="y = x")
+    ax.set(title="Estimated vs true")
+    ax.set_xlabel("true yaw [deg]", fontsize=axis_label_size)
+    ax.set_ylabel("estimated yaw change [deg]", fontsize=axis_label_size)
+    ax.tick_params(axis="both", labelsize=tick_label_size)
+    ax.legend()
+    ax.grid(alpha=0.3)
+    fig.tight_layout()
+    paths.append(out_dir / "accuracy_scatter.png")
+    fig.savefig(paths[-1])
+    plt.close(fig)
+    fig, axes = plt.subplots(len(distances), 1, figsize=(8, 3.2 * len(distances)), squeeze=False)
+    for ax, distance in zip(axes[:, 0], distances):
+        by_angle: dict[float, list[float]] = {}
+        for r, d in ref:
+            if float(r["distance_m"]) == distance:  # type: ignore[arg-type]
+                by_angle.setdefault(float(r["true_yaw_deg"]), []).extend(_errors(r, d))  # type: ignore[arg-type]
+        angles = sorted(by_angle)
+        ax.boxplot([by_angle[a] for a in angles], positions=list(range(len(angles))), showfliers=False)
+        ax.set_xticks(list(range(len(angles))))
+        ax.set_xticklabels([f"{a:+g}" for a in angles])
+        theta = thresholds.get(distance, {}).get("theta_th_deg")
+        ax.axhline(0, color="k", linewidth=0.8)
+        if theta is not None:
+            ax.axhline(theta, color="r", linestyle=":")
+            ax.axhline(-theta, color="r", linestyle=":")
+        ax.set(title=f"Frame error by true yaw ({distance:g} m)")
+        ax.set_xlabel("true yaw [deg]", fontsize=axis_label_size)
+        ax.set_ylabel("error [deg]", fontsize=axis_label_size)
+        ax.tick_params(axis="both", labelsize=tick_label_size)
+    fig.tight_layout()
+    paths.append(out_dir / "accuracy_error_by_condition.png")
+    fig.savefig(paths[-1])
+    plt.close(fig)
+    fig, ax = plt.subplots(figsize=(7, 4))
+    for distance in distances:
+        errors = [e for r, d in ref if float(r["distance_m"]) == distance  # type: ignore[arg-type]
+                  and abs(float(r["true_yaw_deg"])) <= opt.pool_max_angle for e in _errors(r, d)]  # type: ignore[arg-type]
+        if errors:
+            ax.hist(errors, bins=40, histtype="step", density=True, label=f"{distance:g} m")
+        ax.set(title=f"Frame error distribution (|true| <= {opt.pool_max_angle:g} deg)")
+        ax.set_xlabel("error [deg]", fontsize=axis_label_size)
+        ax.set_ylabel("density", fontsize=axis_label_size)
+        ax.tick_params(axis="both", labelsize=tick_label_size)
+    ax.legend()
+    fig.tight_layout()
+    paths.append(out_dir / "accuracy_error_hist.png")
+    fig.savefig(paths[-1])
+    plt.close(fig)
+    return paths
+
+
+def run_report(trial_paths: Sequence[Path], face_paths: Sequence[Path], out_dir: Path, opt: Options,
+               exclude_seqs: Sequence[int] = (), plots: bool = True) -> dict[str, object]:
+    trials = load_trials(trial_paths, exclude_seqs)
+    series = load_face_series(face_paths)
+    evaluated = [evaluate_trial(t, series, opt) for t in trials]
+    records = [rec for rec, _ in evaluated]
+    usable: Usable = [(rec, deltas) for rec, deltas in evaluated if rec["status"] == "ok"]
+    conditions = summarize_conditions(usable)
+    thresholds = compute_thresholds(usable, opt)
+    add_exceed_rates(conditions, usable, thresholds, opt)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    _write_csv(out_dir / "accuracy_trials.csv", TRIAL_COLUMNS, records)
+    _write_csv(out_dir / "accuracy_conditions.csv", CONDITION_COLUMNS, conditions)
+    _write_csv(out_dir / "accuracy_thresholds.csv", THRESHOLD_COLUMNS,
+               [{"distance_m": d, **t} for d, t in sorted(thresholds.items())])
+    markdown = build_markdown(records, conditions, thresholds, opt)
+    (out_dir / "accuracy_report.md").write_text(markdown, encoding="utf-8")
+    figures = make_plots(usable, thresholds, out_dir / "plots", opt) if plots else []
+    return {"records": records, "conditions": conditions, "thresholds": thresholds, "markdown": markdown,
+            "figures": figures}
+
+
+def main(argv: Sequence[str] | None = None, out=print) -> int:  # type: ignore[no-untyped-def]
+    parser = argparse.ArgumentParser(description="Accuracy report from trial markers and face_frames logs")
+    parser.add_argument("--trials", nargs="+", required=True, help="trials CSV file(s) from accuracy_session")
+    parser.add_argument("--face-csv", nargs="+", required=True, help="*_face_frames.csv file(s); globs allowed")
+    parser.add_argument("--out", type=Path, default=Path("data/analysis/accuracy/report"))
+    parser.add_argument("--signal", choices=("yaw", "raw_yaw"), default="yaw")
+    parser.add_argument("--flip-sign", action="store_true")
+    parser.add_argument("--settle", type=float, default=3.0, help="seconds after the cue ignored (turn in progress)")
+    parser.add_argument("--baseline-trim", type=float, default=0.5)
+    parser.add_argument("--min-frames", type=int, default=10)
+    parser.add_argument("--face-every", type=int, default=1, help="face.process_every_n_frames of the run")
+    parser.add_argument("--track-id", type=int, default=None)
+    parser.add_argument("--fps", type=float, default=None, help="override the frame rate used for detection rate")
+    parser.add_argument("--k", type=float, default=3.0)
+    parser.add_argument("--pool-max-angle", type=float, default=15.0)
+    parser.add_argument("--basis", choices=("worst", "pooled"), default="worst")
+    parser.add_argument("--detect-ratio", type=float, default=0.95)
+    parser.add_argument("--exclude-seq", type=int, nargs="*", default=[])
+    parser.add_argument("--no-plots", action="store_true")
+    args = parser.parse_args(argv)
+    trial_paths, face_paths = expand_paths(args.trials), expand_paths(args.face_csv)
+    missing = [str(p) for p in trial_paths + face_paths if not p.is_file()]
+    if missing or not trial_paths or not face_paths:
+        print(f"file not found: {', '.join(missing) or 'no input files'}", file=sys.stderr)
+        return 2
+    opt = Options(args.signal, args.flip_sign, args.settle, args.baseline_trim, args.min_frames, args.face_every,
+                  args.track_id, args.fps, args.k, args.pool_max_angle, args.basis, args.detect_ratio)
+    result = run_report(trial_paths, face_paths, args.out, opt, args.exclude_seq, not args.no_plots)
+    if not result["records"]:
+        print("no usable trials (status ok) in the trials file", file=sys.stderr)
+        return 1
+    out(str(result["markdown"]))
+    out(f"-> {args.out}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

src/accuracy_session.py

@@ -44,6 +44,11 @@ INTERVAL_AUDIO = {
     "block": Path("assets/audio/block.wav"),
     "auto": Path("assets/audio/auto_interval.wav"),
 }
+TRIAL_AUDIO = {
+    "trial_start": Path("assets/audio/1.wav"),
+    "turn": Path("assets/audio/2.wav"),
+    "trial_end": Path("assets/audio/3.wav"),
+}
 # Angle announcement is composed from separate sign and number WAV files.
 # Prepare these files under assets/audio/ using exactly these names.
 ANGLE_SIGN_AUDIO = {
@@ -184,6 +189,14 @@ class IntervalAudioPlayer:
                     self._sounds[kind] = mixer.Sound(str(path))
                 except Exception as exc:
                     LOGGER.error("Failed to load interval audio %s: %s", path, exc)
+            for kind, path in TRIAL_AUDIO.items():
+                if not path.is_file():
+                    LOGGER.warning("Trial audio file does not exist: %s", path)
+                    continue
+                try:
+                    self._sounds[kind] = mixer.Sound(str(path))
+                except Exception as exc:
+                    LOGGER.error("Failed to load trial audio %s: %s", path, exc)
             for sign, path in ANGLE_SIGN_AUDIO.items():
                 if not path.is_file():
                     LOGGER.warning("Angle sign audio file does not exist (%s): %s", sign, path)
@@ -347,7 +360,7 @@ class SessionRunner:
                 continue
             if first_trial_move_pending:
                 self.interval_audio.play("block")
-                self._countdown(20.0, label="initial move")
+                self._countdown(10.0, label="initial move")
                 first_trial_move_pending = False
             elif block_change and self.block_interval > 0:
                 self.interval_audio.play("block")
@@ -374,17 +387,17 @@ class SessionRunner:
 
     def _execute(self, row: PlanRow) -> tuple[str, str, str]:
         t_base = self.clock()
-        self.beep("start")
+        self.interval_audio.play("trial_start")
         self.out(f"  baseline {self.baseline_sec:g} s: face the camera and stay still")
         self.sleep(self.baseline_sec)
         t_cue = self.clock()
-        self.beep("turn")
+        self.interval_audio.play("turn")
         self.out(f"  TURN to {_sign(row.true_yaw_deg)} deg (head only) — you have {self.turn_grace_sec:g} s to turn")
         self.sleep(self.turn_grace_sec)
         self.out(f"  HOLD {self.hold_sec:g} s")
         self.sleep(self.hold_sec)
         t_end = self.clock()
-        self.beep("end")
+        #self.interval_audio.play("trial_end")
         self.out("  return to the front")
         return t_base, t_cue, t_end
 

tests/test_accuracy_report.py

@@ -7,7 +7,8 @@ from pathlib import Path
 import pytest
 
 from src.accuracy_report import (
-    Options, evaluate_trial, load_face_series, load_trials, main, percentile, run_report,
+    Options, evaluate_trial, kruskal_permutation, load_face_series, load_trials, main, percentile, run_report,
+    split_label, _signflip_p,
 )
 from src.accuracy_session import SessionRunner, build_plan
 
@@ -16,8 +17,8 @@ FPS = 30.0
 
 
 class FakeTime:
-    def __init__(self) -> None:
+    def __init__(self, hour: int = 12) -> None:
-        self.now = datetime(2026, 10, 4, 12, 0, 0, tzinfo=TZ)
+        self.now = datetime(2026, 10, 4, hour, 0, 0, tzinfo=TZ)
 
     def sleep(self, seconds: float) -> None:
         self.now += timedelta(seconds=seconds)
@@ -26,10 +27,15 @@ class FakeTime:
         return self.now.isoformat(timespec="milliseconds")
 
 
-def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: float = 0.5, face_dropout: float = 0.0):
+def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: float = 0.5, face_dropout: float = 0.0,
-    """Run a fake session and write a matching face_frames CSV (30 fps, yaw ramps after each cue)."""
+              hour: int = 12, seed: int = 0):
+    """Run a fake session and write a matching face_frames CSV (30 fps, yaw ramps after each cue).
+
+    ``hour`` shifts the wall clock so that several simulated subjects never overlap in time.
+    """
+    tmp_path.mkdir(parents=True, exist_ok=True)
     plan = build_plan([1, 3], [0, 5, -5, 10, -10], 3, seed=3)
-    fake, trials = FakeTime(), tmp_path / "trials.csv"
+    fake, trials = FakeTime(hour), tmp_path / "trials.csv"
     start = fake.now
 
     def operator_pause(_prompt: str) -> str:
@@ -39,10 +45,10 @@ def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise:
     SessionRunner(plan, trials, "S", input_fn=operator_pause, out=lambda _m: None, sleep=fake.sleep, clock=fake.clock,
                   beep=lambda _k: None).run()
     end = fake.now + timedelta(seconds=2)
-    rng = random.Random(0)
+    rng = random.Random(seed)
     marks = [(datetime.fromisoformat(r["t_turn_cue"]), float(r["true_yaw_deg"]), datetime.fromisoformat(r["t_hold_end"]))
              for r in csv.DictReader(trials.open(encoding="utf-8-sig"))]
-    path = tmp_path / "20261004_120000_exp_face_frames.csv"
+    path = tmp_path / f"20261004_{hour:02d}0000_exp_face_frames.csv"
     with path.open("w", encoding="utf-8-sig", newline="") as stream:
         writer = csv.writer(stream)
         writer.writerow(["frame_number", "timestamp", "track_id", "pose_estimated", "yaw", "raw_yaw"])
@@ -83,6 +89,9 @@ def test_full_pipeline_recovers_scale_bias_and_threshold(tmp_path: Path) -> None
     for name in ("accuracy_trials.csv", "accuracy_conditions.csv", "accuracy_thresholds.csv", "accuracy_report.md"):
         assert (tmp_path / "out" / name).is_file()
     assert "θ_th" in result["markdown"]
+    # one subject: no multi-subject sections or files
+    assert result["multi"] is None and "被験者別" not in result["markdown"]
+    assert not (tmp_path / "out" / "accuracy_loso.csv").exists()
 
 
 def test_flip_sign_is_flagged_by_negative_slope(tmp_path: Path) -> None:
@@ -137,3 +146,97 @@ def test_cli_end_to_end_with_glob_and_plots(tmp_path: Path, capsys: pytest.Captu
     assert "しきい値" in capsys.readouterr().out
     assert (out / "plots" / "accuracy_scatter.png").is_file()  # matplotlib is in requirements.txt
     assert main(["--trials", str(trials), "--face-csv", str(tmp_path / "missing.csv")]) == 2
+
+
+# ---------------------------------------------------------------------------------------------
+# multiple subjects
+# ---------------------------------------------------------------------------------------------
+def _two_subjects(tmp_path: Path, scale_b: float = 0.80):
+    ta, fa = _simulate(tmp_path / "a", scale=0.95, offset=-19.0, hour=12)
+    tb, fb = _simulate(tmp_path / "b", scale=scale_b, offset=-5.0, hour=14, seed=1)
+    return ta, fa, tb, fb
+
+
+def test_split_label() -> None:
+    assert split_label("A=trials_a.csv") == ("A", "trials_a.csv")
+    assert split_label(r"C:\data\x=1.csv") == ("", r"C:\data\x=1.csv")
+    assert split_label("trials.csv") == ("", "trials.csv")
+
+
+def test_subject_labels_aliases_and_exclusions(tmp_path: Path) -> None:
+    header = "session_id,seq,block,condition_id,distance_m,true_yaw_deg,pitch_deg,repeat,subject,status,t_baseline_start,t_turn_cue,t_hold_end\n"
+    a, b = tmp_path / "a.csv", tmp_path / "b.csv"
+    a.write_text(header + "S,1,1,c,1,10,0,1,S01,ok,a,b,c\nS,2,1,c,1,10,0,2,S01,ok,a,b,c\n", encoding="utf-8-sig")
+    b.write_text(header + "T,1,1,c,1,10,0,1,,ok,a,b,c\nT,2,1,c,1,10,0,2,,ok,a,b,c\n", encoding="utf-8-sig")
+    assert {r["subject"] for r in load_trials([a, b])} == {"S01", "b"}      # column, then file name
+    rows = load_trials([a, b], labels={b: "S02"}, aliases={"S01": "A"})
+    assert {r["subject"] for r in rows} == {"A", "S02"} and len(rows) == 4
+    kept = load_trials([a, b], exclude_seqs=["S02:2", 1], labels={b: "S02"})
+    assert [(r["subject"], r["seq"]) for r in kept] == [("S01", "2")]
+    only = load_trials([a, b], exclude_seqs=["S02:2"], labels={b: "S02"})
+    assert [(r["subject"], r["seq"]) for r in only] == [("S01", "1"), ("S01", "2"), ("S02", "1")]
+
+
+def test_multi_subject_pools_frames_and_reports_each_subject(tmp_path: Path) -> None:
+    ta, fa, tb, fb = _two_subjects(tmp_path)
+    result = run_report([ta, tb], [fa, fb], tmp_path / "out", Options(), plots=False, labels={ta: "A", tb: "B"})
+    assert result["subjects"] == ["A", "B"] and len(result["records"]) == 60
+    assert all(r["status"] == "ok" for r in result["records"])
+    multi = result["multi"]
+    per = multi["per_subject"]
+    assert abs(per["A"]["thresholds"][1.0]["slope"] - 0.95) < 0.05
+    assert abs(per["B"]["thresholds"][1.0]["slope"] - 0.80) < 0.05
+    pooled = result["thresholds"][1.0]
+    assert pooled["n_subjects"] == 2 and 0.80 < pooled["slope"] < 0.95
+    cond = {(c["distance_m"], c["true_yaw_deg"]): c for c in result["conditions"]}
+    assert cond[(1.0, 10.0)]["n_subjects"] == 2 and cond[(1.0, 10.0)]["n_trials"] == 6
+    assert cond[(1.0, 10.0)]["sd_subjects_deg"] > 0.5        # +10 deg reads 9.5 (A) vs 8.0 (B)
+    # leave-one-subject-out: one row per held-out subject and distance
+    loso = multi["loso"]
+    assert {(r["held_out"], r["distance_m"]) for r in loso} == {("A", 1.0), ("A", 3.0), ("B", 1.0), ("B", 3.0)}
+    assert all(r["n_train_subjects"] == 1 for r in loso)
+    # the +/- gain is symmetric in this simulation; the per-subject gains differ
+    asym = {r["subject"]: r for r in multi["asymmetry"]}
+    assert set(asym) == {"A", "B", "ALL"} and abs(asym["ALL"]["gain_diff"]) < 0.05
+    assert asym["A"]["gain_pos"] > asym["B"]["gain_pos"] + 0.05
+    assert multi["differences"] and all(r["p_permutation"] < 0.05 for r in multi["differences"])
+    for name in ("accuracy_subject_conditions.csv", "accuracy_subject_thresholds.csv", "accuracy_loso.csv",
+                 "accuracy_asymmetry.csv", "accuracy_asymmetry_cells.csv", "accuracy_subject_differences.csv"):
+        assert (tmp_path / "out" / name).is_file(), name
+    header = (tmp_path / "out" / "accuracy_conditions.csv").read_text(encoding="utf-8-sig").splitlines()[0]
+    assert header.endswith("exceed_rate,n_subjects,sd_subjects_deg")
+    for text in ("被験者別しきい値", "leave-one-subject-out", "正負方向の非対称", "A (30試行)", "B (30試行)"):
+        assert text in result["markdown"]
+
+
+def test_flip_only_the_named_subject(tmp_path: Path) -> None:
+    ta, fa, tb, fb = _two_subjects(tmp_path, scale_b=-0.95)
+    labels = {ta: "A", tb: "B"}
+    raw = run_report([ta, tb], [fa, fb], tmp_path / "o1", Options(), plots=False, labels=labels)
+    assert raw["multi"]["per_subject"]["B"]["thresholds"][1.0]["slope"] < 0
+    assert "--flip-subject B" in raw["markdown"]
+    fixed = run_report([ta, tb], [fa, fb], tmp_path / "o2", Options(flip_subjects=("B",)), plots=False, labels=labels)
+    per = fixed["multi"]["per_subject"]
+    assert per["A"]["thresholds"][1.0]["slope"] > 0.9 and per["B"]["thresholds"][1.0]["slope"] > 0.9
+    assert "--flip-subject B" not in fixed["markdown"]
+    only_a = run_report([ta, tb], [fa, fb], tmp_path / "o3", Options(), plots=False, labels=labels,
+                        exclude_subjects=["B"])
+    assert only_a["subjects"] == ["A"] and only_a["multi"] is None
+
+
+def test_cli_multi_subject_with_labels_and_subject_plots(tmp_path: Path, capsys: pytest.CaptureFixture[str]) -> None:
+    ta, fa, tb, fb = _two_subjects(tmp_path)
+    out = tmp_path / "report"
+    assert main(["--trials", f"A={ta}", f"B={tb}", "--face-csv", str(fa), str(fb), "--out", str(out)]) == 0
+    assert "被験者別しきい値" in capsys.readouterr().out
+    for name in ("accuracy_scatter.png", "accuracy_scatter_by_subject.png", "accuracy_error_by_subject.png",
+                 "accuracy_gain_by_side.png"):
+        assert (out / "plots" / name).is_file(), name
+    assert main(["--trials", f"A={ta}", "--face-csv", str(fa), "--alias", "broken"]) == 2
+
+
+def test_permutation_helpers() -> None:
+    assert _signflip_p([1.0, 1.0, 1.0, 1.0]) == pytest.approx(2 / 16)
+    assert math.isnan(_signflip_p([0.0, 0.0, 0.0]))
+    h, p = kruskal_permutation([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
+    assert h == pytest.approx(3.857, abs=0.01) and 0.05 < p < 0.2

tests/test_accuracy_report.py.orig

@@ -0,0 +1,139 @@
+import csv
+import math
+import random
+from datetime import datetime, timedelta, timezone
+from pathlib import Path
+
+import pytest
+
+from src.accuracy_report import (
+    Options, evaluate_trial, load_face_series, load_trials, main, percentile, run_report,
+)
+from src.accuracy_session import SessionRunner, build_plan
+
+TZ = timezone(timedelta(hours=9))
+FPS = 30.0
+
+
+class FakeTime:
+    def __init__(self) -> None:
+        self.now = datetime(2026, 10, 4, 12, 0, 0, tzinfo=TZ)
+
+    def sleep(self, seconds: float) -> None:
+        self.now += timedelta(seconds=seconds)
+
+    def clock(self) -> str:
+        return self.now.isoformat(timespec="milliseconds")
+
+
+def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: float = 0.5, face_dropout: float = 0.0):
+    """Run a fake session and write a matching face_frames CSV (30 fps, yaw ramps after each cue)."""
+    plan = build_plan([1, 3], [0, 5, -5, 10, -10], 3, seed=3)
+    fake, trials = FakeTime(), tmp_path / "trials.csv"
+    start = fake.now
+
+    def operator_pause(_prompt: str) -> str:
+        fake.sleep(4.0)  # time the subject needs to return and the operator needs to press Enter
+        return ""
+
+    SessionRunner(plan, trials, "S", input_fn=operator_pause, out=lambda _m: None, sleep=fake.sleep, clock=fake.clock,
+                  beep=lambda _k: None).run()
+    end = fake.now + timedelta(seconds=2)
+    rng = random.Random(0)
+    marks = [(datetime.fromisoformat(r["t_turn_cue"]), float(r["true_yaw_deg"]), datetime.fromisoformat(r["t_hold_end"]))
+             for r in csv.DictReader(trials.open(encoding="utf-8-sig"))]
+    path = tmp_path / "20261004_120000_exp_face_frames.csv"
+    with path.open("w", encoding="utf-8-sig", newline="") as stream:
+        writer = csv.writer(stream)
+        writer.writerow(["frame_number", "timestamp", "track_id", "pose_estimated", "yaw", "raw_yaw"])
+        t, frame = start - timedelta(seconds=1), 0
+        while t < end:
+            frame += 1
+            true = 0.0
+            for cue, angle, stop in marks:
+                if cue <= t <= stop:
+                    true = angle * min(1.0, (t - cue).total_seconds() / 0.8)
+            yaw = offset + scale * true + rng.gauss(0, noise)
+            ok = rng.random() >= face_dropout
+            writer.writerow([frame, t.isoformat(timespec="milliseconds"), 1, ok, f"{yaw:.4f}" if ok else "", f"{yaw:.4f}" if ok else ""])
+            t += timedelta(seconds=1 / FPS)
+    return trials, path
+
+
+def test_percentile() -> None:
+    assert percentile([0.0, 10.0], 50) == 5.0 and math.isnan(percentile([], 50))
+
+
+def test_full_pipeline_recovers_scale_bias_and_threshold(tmp_path: Path) -> None:
+    trials, face = _simulate(tmp_path)
+    result = run_report([trials], [face], tmp_path / "out", Options(), plots=False)
+    records = result["records"]
+    assert len(records) == 30 and all(r["status"] == "ok" for r in records)
+    # a constant -19 deg offset must cancel; 0.95 scale makes +10 read about 9.5 -> bias about -0.5
+    cond = {(c["distance_m"], c["true_yaw_deg"]): c for c in result["conditions"]}
+    assert abs(cond[(1.0, 0.0)]["bias_deg"]) < 0.3
+    assert abs(cond[(3.0, 10.0)]["bias_deg"] + 0.5) < 0.3
+    th = result["thresholds"][1.0]
+    assert 0.95 - 0.05 < th["slope"] < 0.95 + 0.05 and th["r2"] > 0.99
+    assert 1.0 < th["theta_th_deg"] < 3.5
+    assert th["theta_worst_deg"] >= th["theta_pooled_deg"] - 1e-9
+    assert th["min_detectable_deg"] == 5.0           # 5 deg turns exceed theta_th in >= 95 % of frames
+    assert cond[(1.0, 0.0)]["exceed_rate"] < 0.05    # false-positive rate at 0 deg
+    assert cond[(1.0, 10.0)]["mean_detection_rate"] > 0.95
+    for name in ("accuracy_trials.csv", "accuracy_conditions.csv", "accuracy_thresholds.csv", "accuracy_report.md"):
+        assert (tmp_path / "out" / name).is_file()
+    assert "θ_th" in result["markdown"]
+
+
+def test_flip_sign_is_flagged_by_negative_slope(tmp_path: Path) -> None:
+    trials, face = _simulate(tmp_path, scale=-1.0)
+    result = run_report([trials], [face], tmp_path / "out", Options(), plots=False)
+    assert result["thresholds"][1.0]["slope"] < 0 and "--flip-sign" in result["markdown"]
+    fixed = run_report([trials], [face], tmp_path / "out2", Options(flip_sign=True), plots=False)
+    assert fixed["thresholds"][1.0]["slope"] > 0.9
+
+
+def test_dropout_lowers_detection_rate_and_missing_data_is_excluded(tmp_path: Path) -> None:
+    trials, face = _simulate(tmp_path, face_dropout=0.3)
+    result = run_report([trials], [face], tmp_path / "out", Options(), plots=False)
+    rates = [r["detection_rate"] for r in result["records"] if r["status"] == "ok"]
+    assert 0.55 < sum(rates) / len(rates) < 0.85
+    # a face log that ends early leaves later trials without any track
+    series = load_face_series([face])
+    cut = series.times[len(series.times) // 2]
+    short = tmp_path / "short_face_frames.csv"
+    with face.open(encoding="utf-8-sig", newline="") as src, short.open("w", encoding="utf-8-sig", newline="") as dst:
+        reader, writer = csv.reader(src), csv.writer(dst)
+        for index, row in enumerate(reader):
+            if index == 0 or datetime.fromisoformat(row[1]).timestamp() <= cut:
+                writer.writerow(row)
+    partial = run_report([trials], [short], tmp_path / "out3", Options(), plots=False)
+    statuses = {r["status"] for r in partial["records"]}
+    assert "no_track" in statuses and "ok" in statuses
+
+
+def test_last_row_wins_and_exclude_seq(tmp_path: Path) -> None:
+    path = tmp_path / "t.csv"
+    header = "session_id,seq,block,condition_id,distance_m,true_yaw_deg,pitch_deg,repeat,subject,status,t_baseline_start,t_turn_cue,t_hold_end\n"
+    path.write_text(header + "S,1,1,c,1,10,0,1,,ok,a,b,c\nS,1,1,c,1,10,0,1,,discarded,,,\nS,2,1,c,1,10,0,2,,ok,a,b,c\n"
+                    "S,3,1,c,1,10,0,3,,ok,a,b,c\n", encoding="utf-8-sig")
+    assert [r["seq"] for r in load_trials([path])] == ["2", "3"]
+    assert [r["seq"] for r in load_trials([path], exclude_seqs=[3])] == ["2"]
+
+
+def test_bad_timestamps_do_not_crash(tmp_path: Path) -> None:
+    trials, face = _simulate(tmp_path)
+    series = load_face_series([face])
+    record, deltas = evaluate_trial({"seq": "1", "condition_id": "x", "distance_m": "1", "true_yaw_deg": "5",
+                                     "pitch_deg": "0", "t_baseline_start": "", "t_turn_cue": "", "t_hold_end": ""},
+                                    series, Options())
+    assert record["status"] == "bad_timestamps" and deltas == []
+
+
+def test_cli_end_to_end_with_glob_and_plots(tmp_path: Path, capsys: pytest.CaptureFixture[str]) -> None:
+    trials, face = _simulate(tmp_path)
+    out = tmp_path / "report"
+    assert main(["--trials", str(trials), "--face-csv", str(tmp_path / "*_face_frames.csv"), "--out", str(out)]) == 0
+    assert "しきい値" in capsys.readouterr().out
+    assert (out / "plots" / "accuracy_scatter.png").is_file()  # matplotlib is in requirements.txt
+    assert main(["--trials", str(trials), "--face-csv", str(tmp_path / "missing.csv")]) == 2