ML Uncertainty Quantification in Batch and Sequential Settings: What Changes When Integrating e-values into Conformal Prediction

import numpy as np
from sklearn.ensemble import RandomForestClassifier
 
def split_conformal_predict(clf, cal_X, cal_y, test_X, alpha=0.05):
    """
    Split Conformal Prediction (p-value 방식) 기본 구조.
    보정 집합으로 비적합 점수를 계산하고
    임계값 이상의 레이블을 예측 집합으로 반환합니다.
    """
    # 소프트맥스 확률 기반 비적합 점수: 1 - 정답 클래스의 예측 확률
    cal_probs = clf.predict_proba(cal_X)
    cal_scores = 1 - cal_probs[np.arange(len(cal_y)), cal_y]
 
    # (1-α) 분위수를 임계값으로 사용
    n = len(cal_scores)
    threshold = np.quantile(cal_scores, np.ceil((n + 1) * (1 - alpha)) / n)
 
    # 테스트 데이터에 대해 임계값 이하인 레이블을 예측 집합으로 포함
    test_probs = clf.predict_proba(test_X)
    prediction_sets = []
    for probs in test_probs:
        scores = 1 - probs
        pred_set = np.where(scores <= threshold)[0].tolist()
        prediction_sets.append(pred_set)
 
    return prediction_sets

E[E | H₀] ≤ 1

import numpy as np
from typing import List
 
def compute_conformal_evalue(cal_scores: np.ndarray, test_score: float) -> float:
    """
    보정 집합의 비적합 점수와 테스트 점수로부터 e-value를 계산합니다.
 
    직관: 테스트 점수가 보정 집합보다 얼마나 극단적인가를 '배율'로 표현합니다.
    E = 1/p 변환은 마르코프 부등식으로 유효성이 보장되지만,
    더 정교한 베팅 전략을 사용하면 검출력(power)이 향상됩니다.
    """
    n = len(cal_scores)
    # 테스트 점수보다 크거나 같은 보정 점수의 비율 → p-value
    p_value = (np.sum(cal_scores >= test_score) + 1) / (n + 1)
    return 1.0 / p_value
 
 
def conformal_evalue_predict(clf, cal_X: np.ndarray, cal_y: np.ndarray,
                              test_X: np.ndarray,
                              threshold: float = 1.0) -> List[List[int]]:
    """
    e-value 기반 예측 집합 구성.
    e-value > threshold인 레이블을 집합에 포함합니다.
    """
    cal_probs = clf.predict_proba(cal_X)
    test_probs_all = clf.predict_proba(test_X)  # 배치 추론으로 사전 계산
    n_classes = cal_probs.shape[1]
 
    prediction_sets = []
    for test_probs in test_probs_all:
        pred_set = []
        for y_candidate in range(n_classes):
            cal_scores_y = 1 - cal_probs[cal_y == y_candidate, y_candidate]
            test_score_y = 1 - test_probs[y_candidate]
 
            if len(cal_scores_y) > 0:
                e_val = compute_conformal_evalue(cal_scores_y, test_score_y)
                if e_val > threshold:
                    pred_set.append(y_candidate)
 
        prediction_sets.append(pred_set)
 
    return prediction_sets

import numpy as np
from dataclasses import dataclass, field
from typing import List, Tuple
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
 
# ── 데이터 준비 ──────────────────────────────────────────────────
X, y = make_classification(n_samples=2000, n_features=20, random_state=42)
X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, test_size=0.6, random_state=42
)
X_cal, X_test_all, y_cal, y_test_all = train_test_split(
    X_temp, y_temp, test_size=0.5, random_state=42
)
 
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)
 
# 순차 배치 시뮬레이션 — 5개 병원에서 데이터가 순차 도착
n_batches = 5
batch_size = len(X_test_all) // n_batches
sequential_batches: List[Tuple[np.ndarray, np.ndarray]] = [
    (
        X_test_all[i * batch_size:(i + 1) * batch_size],
        y_test_all[i * batch_size:(i + 1) * batch_size],
    )
    for i in range(n_batches)
]
 
 
# ── 배치 e-value 계산 ────────────────────────────────────────────
 
def compute_batch_evalue(clf, cal_X: np.ndarray, cal_y: np.ndarray,
                          batch_X: np.ndarray, batch_y: np.ndarray) -> float:
    """
    배치 데이터로부터 배치 e-value를 계산합니다.
 
    구현: 각 샘플의 e-value(= 1/p-value)의 기하평균을 사용합니다.
    이론적으로 독립 e-value의 곱도 유효한 e-value이며,
    기하평균은 배치 크기 차이를 보정한 샘플당 증거 강도입니다.
 
    ⚠️ 주의: arXiv:2503.13050의 배치 알고리즘은 베팅 함수(betting function)를
    명시적으로 최적화합니다. 이 단순 구현은 이론적 보장은 유지하지만
    검출력이 낮을 수 있습니다. 프로덕션 적용 전 원논문 확인을 권장합니다.
    """
    cal_probs = clf.predict_proba(cal_X)
    cal_scores = 1 - cal_probs[np.arange(len(cal_y)), cal_y]
 
    batch_probs = clf.predict_proba(batch_X)
    batch_scores = 1 - batch_probs[np.arange(len(batch_y)), batch_y]
 
    n_cal = len(cal_scores)
    p_values = np.array([
        (np.sum(cal_scores >= s) + 1) / (n_cal + 1)
        for s in batch_scores
    ])
 
    # 각 샘플의 e-value 기하평균 (배치 크기 정규화)
    individual_evalues = 1.0 / p_values
    batch_evalue = float(np.exp(np.mean(np.log(individual_evalues))))
    return batch_evalue
 
 
# ── 애니타임 모니터 ──────────────────────────────────────────────
 
@dataclass
class BatchEValueMonitor:
    """
    배치 단위로 순차 도착하는 데이터에 대해
    e-value를 누적 곱으로 관리하는 모니터.
    언제든 검정을 수행할 수 있는 애니타임 유효성을 제공합니다.
    """
    alpha: float = 0.05
    cumulative_evalue: float = 1.0
    batch_history: List[float] = field(default_factory=list)
 
    def update(self, batch_evalue: float) -> dict:
        """
        새 배치의 e-value로 누적값을 업데이트합니다.
        귀무가설 기각 여부와 현재 유효 오류율을 반환합니다.
        """
        self.cumulative_evalue *= batch_evalue
        self.batch_history.append(batch_evalue)
 
        # Markov 부등식: P(E_cumul ≥ 1/α) ≤ α
        reject = self.cumulative_evalue >= (1.0 / self.alpha)
        current_alpha = 1.0 / self.cumulative_evalue
 
        return {
            "batch_num": len(self.batch_history),
            "batch_evalue": batch_evalue,
            "cumulative_evalue": self.cumulative_evalue,
            "reject_null": reject,
            "current_effective_alpha": min(current_alpha, 1.0),
        }
 
 
# ── 순차 실행 예시 ───────────────────────────────────────────────
 
monitor = BatchEValueMonitor(alpha=0.05)
 
for batch_X, batch_y in sequential_batches:
    e_val = compute_batch_evalue(clf, X_cal, y_cal, batch_X, batch_y)
    result = monitor.update(e_val)
 
    print(f"배치 {result['batch_num']:2d} | "
          f"e-value: {result['batch_evalue']:.3f} | "
          f"누적: {result['cumulative_evalue']:.3f} | "
          f"기각: {result['reject_null']}")
 
    if result["reject_null"]:
        print("→ 귀무가설 기각: 통계적으로 유의미한 효과 확인됨")
        break  # 이 시점 이후 추가 배치가 와도 이미 얻은 보장은 유효합니다

import numpy as np
from typing import List, Optional, Dict
 
def compute_fuzzy_membership(e_value: float, threshold: float = 1.0) -> float:
    """
    e-value를 [0, 1] 범위의 퍼지 멤버십 값으로 변환합니다.
 
    변환 공식: log(1 + (e_value - threshold)) / log(10)
    - threshold 미만 → 0.0 (예측 집합 미포함)
    - e_value = threshold + 9일 때 → 1.0 (포화 기준)
    - log 스케일로 e-value가 커질수록 포화되는 형태 → 극단값 안정화
 
    포화 기준(+9)은 "e-value가 threshold보다 9 높을 때 완전 포함"을 의미하는
    임의 상수입니다. 엄격한 의료 진단처럼 빠른 포화가 필요하다면
    np.log1p(4.0)으로 낮추고, 넓은 분포의 데이터라면 더 크게 설정할 수 있습니다.
    """
    if e_value < threshold:
        return 0.0
    return min(1.0, np.log1p(e_value - threshold) / np.log1p(9.0))
 
 
def fuzzy_conformal_predict(clf, cal_X: np.ndarray, cal_y: np.ndarray,
                             test_X: np.ndarray,
                             ambiguous_labels: Optional[np.ndarray] = None
                             ) -> List[Dict[int, float]]:
    """
    퍼지 컨포멀 예측 집합을 구성합니다.
 
    ambiguous_labels: (n_cal_samples, n_classes) 형태의 소프트 레이블 행렬.
                      각 셀은 해당 클래스에 대한 전문가 동의율(0~1)입니다.
                      None이면 하드 레이블을 원-핫으로 처리합니다.
    반환: 각 테스트 샘플에 대한 {레이블: 멤버십 값} 딕셔너리 리스트
    """
    cal_probs = clf.predict_proba(cal_X)
    test_probs_all = clf.predict_proba(test_X)  # 배치 추론으로 사전 계산
    n_classes = cal_probs.shape[1]
 
    if ambiguous_labels is None:
        ambiguous_labels = np.eye(n_classes)[cal_y]
 
    results = []
    for test_probs in test_probs_all:
        memberships: Dict[int, float] = {}
 
        for y_candidate in range(n_classes):
            soft_weights = ambiguous_labels[:, y_candidate]
            cal_scores = 1 - cal_probs[:, y_candidate]
            weighted_scores = cal_scores * soft_weights
 
            valid_mask = soft_weights > 0
            if valid_mask.sum() < 5:  # 최소 샘플 수 보장
                memberships[y_candidate] = 0.0
                continue
 
            test_score = 1 - test_probs[y_candidate]
            e_val = compute_conformal_evalue(weighted_scores[valid_mask], test_score)
            memberships[y_candidate] = compute_fuzzy_membership(e_val)
 
        results.append(memberships)
 
    return results
 
 
# ── 사용 예시: 3명의 전문가 레이블이 있는 의료 진단 ──────────────
 
n_classes = clf.n_classes_
n_experts = 3
 
# 전문가 레이블 시뮬레이션 — 실제로는 어노테이션 툴에서 가져옵니다
np.random.seed(42)
expert_labels = np.random.randint(0, n_classes, size=(len(y_cal), n_experts))
 
# 클래스별 소프트 레이블 행렬 구성
# ambiguous[i, c] = 샘플 i에 대해 클래스 c로 분류한 전문가 비율
ambiguous = np.zeros((len(y_cal), n_classes))
for c in range(n_classes):
    ambiguous[:, c] = np.mean(expert_labels == c, axis=1)
# 각 행의 합 = 1.0 (클래스별 동의율의 합)
 
fuzzy_preds = fuzzy_conformal_predict(
    clf, X_cal, y_cal, X_test_all[:10],
    ambiguous_labels=ambiguous,
)
 
for i, memberships in enumerate(fuzzy_preds[:3]):
    print(f"\n테스트 샘플 {i}:")
    for label, membership in sorted(memberships.items(), key=lambda x: -x[1]):
        bar = "█" * int(membership * 20)
        print(f"  클래스 {label}: {membership:.3f} |{bar:<20}|")

import matplotlib.pyplot as plt
 
cal_sizes = [100, 500, 1000]
coverages = []
 
for cal_size in cal_sizes:
    cal_X_sub = X_cal[:cal_size]
    cal_y_sub = y_cal[:cal_size]
    pred_sets = conformal_evalue_predict(clf, cal_X_sub, cal_y_sub, X_test_all)
    coverage = np.mean([
        y_true in pred_set
        for y_true, pred_set in zip(y_test_all, pred_sets)
    ])
    coverages.append(coverage)
 
plt.figure(figsize=(8, 4))
plt.plot(cal_sizes, coverages, marker="o", label="실제 커버리지")
plt.axhline(y=0.95, color="red", linestyle="--", label="목표 커버리지 (95%)")
plt.xlabel("보정 집합 크기")
plt.ylabel("실제 커버리지")
plt.title("보정 집합 크기와 커버리지 관계")
plt.legend()
plt.tight_layout()
plt.show()