モデルの適用範囲（AD）とは？材料データでの設定・検証方法をPythonで解説

"""Reproducible numeric-feature AD tutorial; all demo observations are synthetic.
Run: python ad_practical.py --demo --out ad_output
Real CSV: python ad_practical.py --data my_data --out my_results
Required: train.csv, calibration.csv, test.csv, candidates.csv.
"""
from pathlib import Path
import argparse
import json
import numpy as np
import pandas as pd
import sklearn
import joblib
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import NearestNeighbors
from sklearn.ensemble import RandomForestRegressor
from sklearn.svm import OneClassSVM
from sklearn.metrics import mean_squared_error, mean_absolute_error

FEATURES = ["temperature_c", "time_min", "solids_pct"]
TARGET = "strength_mpa"
K = 10
TOLERANCE = 2.5  # illustrative engineering requirement, not a universal value
MIN_CALIBRATION = 25
ABS_ERROR_LIMIT = 2.5  # separate per-observation error limit, illustrative only

def demo(folder):
    rng = np.random.default_rng(42)
    folder.mkdir(parents=True, exist_ok=False)
    def observations(n_inside, n_outside, prefix):
        n = n_inside + n_outside
        temp = np.r_[rng.uniform(70, 110, n_inside), rng.uniform(115, 145, n_outside)]
        time = rng.uniform(15, 45, n)
        solids = rng.uniform(20, 40, n)
        # Fictional polymer-film relationship and iid observation noise.
        strength = (35 + .18*(temp-90) + .12*(time-30) + .30*(solids-30)
                    - .018*np.maximum(temp-110, 0)**2
                    + .002*(temp-90)*(time-30) + rng.normal(0, 1.2, n))
        return pd.DataFrame({"sample_id": [f"{prefix}{i:03d}" for i in range(n)],
            FEATURES[0]: temp, FEATURES[1]: time, FEATURES[2]: solids, TARGET: strength})
    observations(180, 0, "TR").to_csv(folder/"train.csv", index=False)
    observations(80, 40, "CA").to_csv(folder/"calibration.csv", index=False)
    observations(100, 50, "TE").to_csv(folder/"test.csv", index=False)
    pd.DataFrame({"sample_id": ["A", "B", "C", "D"],
        "temperature_c": [90, 140, 75, 105], "time_min": [30, 30, 44, 18],
        "solids_pct": [30, 30, 39, 23]}).to_csv(folder/"candidates.csv", index=False)

def validate_frame(frame, features, target=None):
    required = ["sample_id", *features] + ([target] if target else [])
    missing = set(required) - set(frame.columns)
    if missing:
        raise ValueError(f"Missing columns: {missing}")
    if frame.empty or frame.sample_id.isna().any():
        raise ValueError("Empty data or missing sample_id")
    frame = frame.copy()
    frame["sample_id"] = frame.sample_id.astype(str).str.strip()
    if frame.sample_id.eq("").any() or frame.sample_id.duplicated().any():
        raise ValueError("Blank or duplicate sample_id")
    numeric = features + ([target] if target else [])
    frame[numeric] = frame[numeric].apply(pd.to_numeric, errors="raise")
    if not np.isfinite(frame[numeric].to_numpy(dtype=float)).all():
        raise ValueError("Numeric columns contain missing/infinite values")
    return frame


def add_decision_columns(result, metadata):
    # Aggregate validation evidence is separate from each candidate's distance.
    confirmed = (metadata["calibration_criterion_met"] and
                 metadata["test_inside_criterion_met"])
    result = result.copy()
    result["validation_status"] = "criterion_met" if confirmed else "provisional"
    result["review_reason"] = np.select(
        [np.full(len(result), not confirmed), ~result.ad_inside,
         result.outside_any_feature_range],
        ["validation_not_established", "outside_distance_boundary",
         "outside_feature_range_check"],
        default="check_material_and_measurement_context")
    return result


def predict_saved(bundle_path, candidate_path, out):
    # Load only trusted bundles produced by your own workflow.
    bundle = joblib.load(bundle_path)
    metadata = bundle["metadata"]
    if metadata["sklearn_version"] != sklearn.__version__:
        raise ValueError("Use the recorded scikit-learn version for this saved bundle")
    features = metadata["features"]
    frame = validate_frame(pd.read_csv(candidate_path, dtype={"sample_id": str}), features)
    z = bundle["scaler"].transform(frame[features])
    distances, indices = bundle["neighbors"].kneighbors(z)
    result = frame.copy()
    result["prediction_mpa"] = bundle["regression_model"].predict(z)
    result["knn_distance"] = distances.mean(axis=1)
    result["threshold"] = metadata["selected_threshold"]
    result["ad_inside"] = result.knn_distance <= result.threshold
    ids = np.asarray(bundle["training_sample_ids"])
    result["nearest_sample_id"] = ids[indices[:, 0]]
    result["neighbor_ids"] = [";".join(ids[ix]) for ix in indices]
    result["outside_any_feature_range"] = (
        (frame[features] < pd.Series(metadata["training_min"])) |
        (frame[features] > pd.Series(metadata["training_max"]))).any(axis=1)
    result = add_decision_columns(result, metadata)
    out.mkdir(parents=True, exist_ok=False)
    result.to_csv(out/"candidates_results.csv", index=False)
    print(result.to_string(index=False))
    return result


def load(folder):
    frames = {name: pd.read_csv(folder/f"{name}.csv", dtype={"sample_id": str, "batch_id": str})
              for name in ["train", "calibration", "test", "candidates"]}
    for name, frame in frames.items():
        frames[name] = validate_frame(frame, FEATURES, None if name == "candidates" else TARGET)
    # If batch IDs are supplied, all labeled files must include them.
    labeled = [frames[n] for n in ["train", "calibration", "test"]]
    if any("batch_id" in frame for frame in labeled):
        if not all("batch_id" in frame for frame in labeled):
            raise ValueError("batch_id must appear in all three labeled files")
        for frame in labeled:
            if frame.batch_id.isna().any():
                raise ValueError("Missing batch_id")
            frame["batch_id"] = frame.batch_id.astype(str).str.strip()
            if frame.batch_id.eq("").any():
                raise ValueError("Blank batch_id")
        for i, j in [(0, 1), (0, 2), (1, 2)]:
            if set(labeled[i].batch_id) & set(labeled[j].batch_id):
                raise ValueError("Batch IDs overlap across labeled files")
    for first, second in [("train", "calibration"), ("train", "test"), ("calibration", "test")]:
        if set(frames[first].sample_id) & set(frames[second].sample_id):
            raise ValueError(f"{first}/{second}: sample IDs overlap")
    if len(frames["train"]) <= K:
        raise ValueError("Training count must exceed K for self-excluded distances")
    if (frames["train"][FEATURES].nunique() < 2).any():
        raise ValueError("Constant training features: remove with a separate change check")
    if frames["train"].duplicated(FEATURES).any():
        raise ValueError("Duplicate training feature rows: review independent batches/replicates")
    return frames

def metrics(y, prediction, mask, label):
    count = int(mask.sum())
    return {"subset": label, "n": count,
        "rmse_mpa": float(np.sqrt(mean_squared_error(y[mask], prediction[mask]))) if count else None,
        "mae_mpa": float(mean_absolute_error(y[mask], prediction[mask])) if count else None,
        "error_limit_mpa": ABS_ERROR_LIMIT,
        "error_exceedance_fraction": float((np.abs(y[mask]-prediction[mask]) > ABS_ERROR_LIMIT).mean()) if count else None}

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--demo", action="store_true")
    parser.add_argument("--data", default="ad_demo_data")
    parser.add_argument("--out", default="ad_output")
    parser.add_argument("--bundle", help="Trusted saved bundle; use with --candidates")
    parser.add_argument("--candidates", help="New candidate CSV; no strength labels required")
    args = parser.parse_args()
    if args.bundle:
        if not args.candidates or args.demo:
            parser.error("--bundle requires --candidates and cannot combine with --demo")
        predict_saved(args.bundle, args.candidates, Path(args.out))
        return
    if args.candidates:
        parser.error("--candidates requires --bundle")
    folder, out = Path(args.data), Path(args.out)
    if args.demo:
        demo(folder)  # generate only in a new folder
    frames = load(folder)
    out.mkdir(parents=True, exist_ok=False)
    train = frames["train"]
    scaler = StandardScaler().fit(train[FEATURES])  # train only
    z = {name: scaler.transform(frame[FEATURES]) for name, frame in frames.items()}
    nn = NearestNeighbors(n_neighbors=K).fit(z["train"])
    # X=None queries training samples and excludes each row itself.
    train_distances, _ = nn.kneighbors(n_neighbors=K)
    train_scores = train_distances.mean(axis=1)
    preliminary = float(np.quantile(train_scores, .95))
    model = RandomForestRegressor(n_estimators=200, min_samples_leaf=3,
                                  random_state=42, n_jobs=1).fit(z["train"], train[TARGET])
    predictions = {name: model.predict(z[name]) for name in ["calibration", "test", "candidates"]}
    distances, indices, scores = {}, {}, {}
    for name in ["calibration", "test", "candidates"]:
        distances[name], indices[name] = nn.kneighbors(z[name])
        scores[name] = distances[name].mean(axis=1)
    # Calibrate using only calibration labels, with fixed K/model/tolerance.
    cal_y = frames["calibration"][TARGET].to_numpy()
    grid = np.unique(np.r_[np.quantile(train_scores, [.5, .75, .9, .95, .99]),
                           np.quantile(scores["calibration"], [.5, .6, .7, .8, .9, 1])])
    rows = []
    for threshold in grid:
        mask = scores["calibration"] <= threshold
        row = metrics(cal_y, predictions["calibration"], mask, "calibration")
        row.update(threshold=float(threshold), accepted_fraction=float(mask.mean()))
        rows.append(row)
    tuning = pd.DataFrame(rows)
    eligible = tuning[(tuning.n >= MIN_CALIBRATION) & (tuning.rmse_mpa <= TOLERANCE)]
    if eligible.empty:
        chosen, established = preliminary, False
    else:
        chosen = float(eligible.sort_values(["accepted_fraction", "threshold"]).iloc[-1].threshold)
        established = True
    tuning.to_csv(out/"calibration_thresholds.csv", index=False)
    result_tables = {}
    for name in ["calibration", "test", "candidates"]:
        result = frames[name].copy()
        result["prediction_mpa"] = predictions[name]
        result["knn_distance"] = scores[name]
        result["threshold"] = chosen
        result["ad_inside"] = scores[name] <= chosen
        result["nearest_sample_id"] = train.sample_id.iloc[indices[name][:, 0]].to_numpy()
        result["neighbor_ids"] = [";".join(train.sample_id.iloc[ix]) for ix in indices[name]]
        result["outside_any_feature_range"] = (
            (result[FEATURES] < train[FEATURES].min()) |
            (result[FEATURES] > train[FEATURES].max())).any(axis=1)
        if name != "candidates":
            result["absolute_error_mpa"] = np.abs(result[TARGET] - result.prediction_mpa)
        result.to_csv(out/f"{name}_results.csv", index=False)
        result_tables[name] = result
    test_y = frames["test"][TARGET].to_numpy()
    inside = scores["test"] <= chosen
    summary = [metrics(test_y, predictions["test"], np.ones(len(test_y), dtype=bool), "all"),
               metrics(test_y, predictions["test"], inside, "AD_inside"),
               metrics(test_y, predictions["test"], ~inside, "AD_outside")]
    pd.DataFrame(summary).to_csv(out/"test_summary.csv", index=False)
    # Final-test diagnostic curve only. Never use it to tune/reselect the threshold.
    order = np.argsort(scores["test"])
    curve = []
    for n in range(min(25, len(order)), len(order)+1):
        ids = order[:n]
        curve.append({"n": n, "accepted_fraction": n/len(order),
                      "rmse_mpa": float(np.sqrt(mean_squared_error(test_y[ids], predictions["test"][ids])))})
    pd.DataFrame(curve).to_csv(out/"test_coverage_rmse.csv", index=False)
    # Illustrative alternative, not optimized and not used for threshold selection.
    oc = OneClassSVM(kernel="rbf", gamma="scale", nu=.05).fit(z["train"])
    comparison = frames["candidates"][["sample_id"]].copy()
    comparison["knn_inside"] = scores["candidates"] <= chosen
    comparison["ocsvm_inside"] = oc.predict(z["candidates"]) == 1
    comparison["ocsvm_score"] = oc.decision_function(z["candidates"])
    comparison.to_csv(out/"ocsvm_comparison.csv", index=False)
    metadata = {"features": FEATURES, "k": K, "preliminary_threshold": preliminary,
        "training_min": train[FEATURES].min().to_dict(),
        "training_max": train[FEATURES].max().to_dict(),
        "absolute_error_limit_mpa": ABS_ERROR_LIMIT,
        "selected_threshold": chosen, "calibration_criterion_met": established,
        "tolerance_rmse_mpa": TOLERANCE, "min_calibration_n": MIN_CALIBRATION,
        "scaler_mean": scaler.mean_.tolist(), "scaler_scale": scaler.scale_.tolist(),
        "sklearn_version": sklearn.__version__, "demo": args.demo,
        "test_inside_fraction": float(inside.mean()),
        "test_inside_criterion_met": bool(inside.sum() >= MIN_CALIBRATION and
            summary[1]["rmse_mpa"] <= TOLERANCE)}
    for name, result in result_tables.items():
        result_tables[name] = add_decision_columns(result, metadata)
        result_tables[name].to_csv(out/f"{name}_results.csv", index=False)
    group_rows = []
    for column in ["batch_id", "material_family", "measurement_date"]:
        if column in frames["test"]:
            groups = frames["test"][column].fillna("<missing>").astype(str)
            for group in sorted(groups.unique()):
                mask = (groups == group).to_numpy() & inside
                row = metrics(test_y, predictions["test"], mask, "AD_inside")
                row.update(group_column=column, group_value=group)
                group_rows.append(row)
    if group_rows:
        pd.DataFrame(group_rows).to_csv(out/"test_group_summary.csv", index=False)
    (out/"metadata.json").write_text(json.dumps(metadata, ensure_ascii=False, indent=2))
    joblib.dump({"scaler": scaler, "regression_model": model, "neighbors": nn,
                 "ocsvm": oc, "training_sample_ids": train.sample_id.tolist(),
                 "metadata": metadata}, out/"ad_bundle.joblib")
    fig, ax = plt.subplots(figsize=(6, 3.5))
    ax.scatter(scores["test"], np.abs(test_y-predictions["test"]), s=14, alpha=.7)
    ax.axvline(chosen, color="red", linestyle="--", label="Frozen threshold")
    ax.set(xlabel="Mean k-neighbor distance", ylabel="Absolute error (MPa)")
    ax.legend(); fig.tight_layout(); fig.savefig(out/"distance_error.svg")
    fig.savefig(out/"distance_error.png", dpi=180)
    plt.close(fig)
    fig, ax = plt.subplots(figsize=(6, 3.5))
    ax.plot([r["accepted_fraction"] for r in curve], [r["rmse_mpa"] for r in curve])
    ax.axhline(TOLERANCE, color="red", linestyle="--", label="Illustrative tolerance")
    ax.set(xlabel="Accepted fraction (final test)", ylabel="RMSE (MPa)")
    ax.legend(); fig.tight_layout(); fig.savefig(out/"coverage_rmse.svg")
    fig.savefig(out/"coverage_rmse.png", dpi=180)
    plt.close(fig)
    print(json.dumps(metadata, ensure_ascii=False, indent=2))
    print(pd.DataFrame(summary).round(3).to_string(index=False))
    print(result_tables["candidates"][["sample_id", "prediction_mpa", "knn_distance", "ad_inside",
                                     "nearest_sample_id"]].round(3).to_string(index=False))

if __name__ == "__main__":
    main()