假设：
  - 代理 IP 成本：$0.02/次
  - 摸清决策边界大概要 150 次请求
  - 成功率 20%，平均 5 次成 1 次
  - 单次欺诈得手能赚 $50

总攻击成本 = 150 × $0.02（摸索阶段）+ 5 × $0.02（实际攻击阶段）= $3.10
平均 5 次攻击成 1 次，一次得手收益 $50
期望净收益 = $50 - $3.10 = $46.90

ROI ≈ 1513% → 这生意值得做

import numpy as np
from scipy.stats import skew, kurtosis, norm, expon
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import Callable

@dataclass
class AdversarialBehaviorSimulator:
    """
    对抗性行为模拟器：生成高斯噪声扰动后的行为时序

    模拟黑产的攻击策略：
    1. 录制真实用户的原始轨迹（或从统计分布中采样）
    2. 逐点叠加高斯噪声 N(0, σ²)
    3. 调整 σ 使扰动后的特征恰好绕过风控模型
    """
    noise_std: float = 0.0  # 高斯噪声标准差

    def generate_human_like_speeds(self, n_points: int = 200) -> np.ndarray:
        """
        生成仿真的人类鼠标速度序列

        真实人类鼠标速度分布近似对数正态（均速约 300px/s，标准差约 150px/s），
        偶有停顿（速度降至 0-50px/s）和快速滑动（速度飙升到 800-1200px/s）。
        """
        # 对数正态分布：模拟正常移动速度
        speeds = np.random.lognormal(mean=5.7, sigma=0.5, size=n_points)
        # 加入长尾：10% 的快速滑动
        fast_indices = np.random.choice(n_points, size=n_points // 10, replace=False)
        speeds[fast_indices] += np.random.exponential(scale=300, size=n_points // 10)
        # 加入停顿：15% 的静止/低速段
        pause_indices = np.random.choice(n_points, size=n_points // 7, replace=False)
        speeds[pause_indices] = np.random.uniform(0, 50, size=n_points // 7)
        return np.clip(speeds, 0, None)

    def generate_human_like_intervals(self, n_clicks: int = 30) -> np.ndarray:
        """
        生成仿真的人类点击间隔序列

        真实点击间隔近似指数分布（均值 3s，标准差 4s），
        但存在"思考停顿"（10-30s）和"快速连点"（0.1-0.5s）。
        """
        intervals = np.random.exponential(scale=3.0, size=n_clicks)

        # 插入思考停顿（5% 的间隔放大 3-5 倍）
        think_indices = np.random.choice(n_clicks, size=max(1, n_clicks // 20),
                                          replace=False)
        intervals[think_indices] *= np.random.uniform(3, 5, size=len(think_indices))

        # 插入快速连点（10% 的间隔缩小到 0.1-0.5s）
        fast_indices = np.random.choice(
            list(set(range(n_clicks)) - set(think_indices)),
            size=n_clicks // 10, replace=False
        )
        intervals[fast_indices] = np.random.uniform(0.1, 0.5, size=len(fast_indices))

        return np.clip(intervals, 0.05, None)

    def inject_noise(self, signal: np.ndarray) -> np.ndarray:
        """对行为信号注入高斯噪声（黑产的核心攻击手段）"""
        if self.noise_std == 0:
            return signal.copy()
        noise = np.random.normal(0, self.noise_std, size=len(signal))
        return signal + noise

    def generate_adversarial_batch(self, n_samples: int = 1000,
                                    n_points: int = 200
                                    ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
        """
        生成一批对抗性样本用于评估模型鲁棒性

        Returns:
            features: 提取的统计特征矩阵 [n_samples, n_features]
            raw_signals: 原始时序（含噪声）[n_samples, n_points]
            labels: 标签（0=干净真人, 1=对抗样本）
        """
        features = []
        raw_signals = []
        labels = []

        for i in range(n_samples):
            if i < n_samples // 2:
                # 干净真人样本
                speeds = self.generate_human_like_speeds(n_points)
                intervals = self.generate_human_like_intervals(n_points // 6)
                labels.append(0)
            else:
                # 对抗样本：真人轨迹 + 高斯噪声扰动
                speeds = self.generate_human_like_speeds(n_points)
                speeds = self.inject_noise(speeds)  # 注入噪声
                intervals = self.generate_human_like_intervals(n_points // 6)
                intervals = self.inject_noise(intervals)  # 注入噪声
                labels.append(1)

            feat = self._extract_statistical_features(speeds, intervals)
            features.append(feat)
            raw_signals.append(speeds)

        return np.array(features), np.array(raw_signals), np.array(labels)

    @staticmethod
    def _extract_statistical_features(speeds: np.ndarray,
                                       intervals: np.ndarray) -> np.ndarray:
        """从原始时序中提取统计特征（高阶矩 + 熵）"""
        features = []

        # 速度特征
        features.extend([
            np.mean(speeds),           # 均值
            np.std(speeds),            # 标准差
            skew(speeds),              # 偏度
            kurtosis(speeds, fisher=False),  # 峰度（Pearson）
            np.percentile(speeds, 95),  # 95 分位数
            np.percentile(speeds, 5),   # 5 分位数
            np.max(speeds) - np.min(speeds),  # 极差
        ])

        # 速度的自相关（lag=1 的一阶自相关系数）
        if len(speeds) > 1:
            speeds_centered = speeds - np.mean(speeds)
            autocorr = np.corrcoef(speeds_centered[:-1],
                                    speeds_centered[1:])[0, 1]
            features.append(autocorr if not np.isnan(autocorr) else 0)
        else:
            features.append(0)

        # 点击间隔特征
        if len(intervals) > 1:
            features.extend([
                np.mean(intervals),
                np.std(intervals),
                np.std(intervals) / (np.mean(intervals) + 1e-8),  # 变异系数 CV
                skew(intervals),
                kurtosis(intervals, fisher=False),
            ])
        else:
            features.extend([0, 0, 0, 0, 0])

        # 近似熵（衡量序列不可预测性）
        features.append(AdversarialBehaviorSimulator._approximate_entropy(speeds))

        return np.array(features)

    @staticmethod
    def _approximate_entropy(signal: np.ndarray, m: int = 2,
                              r_factor: float = 0.2) -> float:
        """计算近似熵（与分布式风控文章中的实现一致，此处复用）"""
        N = len(signal)
        if N < m + 2:
            return 0.0
        r = r_factor * np.std(signal) + 1e-8

        def _phi(m_val):
            patterns = np.array([signal[i:i+m_val] for i in range(N - m_val + 1)])
            counts = np.sum(
                np.max(np.abs(patterns[:, None] - patterns), axis=2) <= r,
                axis=1
            ) - 1
            counts = np.maximum(counts, 1)
            return np.mean(np.log(counts / (N - m_val + 1)))

        return _phi(m) - _phi(m + 1)


# ---------- 攻击效果的定量分析 ----------
def analyze_attack_impact():
    """
    分析不同噪声强度 σ 对特征统计量的影响

    回答核心问题：黑产需要多大的 σ 才能绕过模型？
    答案是：刚好破坏让你区分"人类 vs 脚本"的那个统计量。
    对一阶特征（均值），高斯噪声加多少都几乎不变——它天然免疫，
    但也因此不携带可检测的攻击信号；
    对二阶特征（方差），需要 σ ≈ σ_human/3 才能明显抬高；
    对三阶以上特征（偏度、峰度），需要 σ ≈ σ_human 才能把它们拉向正态。
    而三阶以上特征恰好是区分人类 vs 脚本的最强信号。
    """
    sim = AdversarialBehaviorSimulator()

    results = []
    for sigma in [0, 20, 50, 100, 200, 400, 800]:
        sim.noise_std = sigma
        speeds = sim.generate_human_like_speeds(1000)
        speeds_noisy = sim.inject_noise(speeds)

        results.append({
            "sigma": sigma,
            "mean": np.mean(speeds_noisy),
            "std": np.std(speeds_noisy),
            "skewness": skew(speeds_noisy),
            "kurtosis": kurtosis(speeds_noisy, fisher=False),
            "range": np.max(speeds_noisy) - np.min(speeds_noisy),
            "ap_en": sim._approximate_entropy(speeds_noisy)
        })

    print("=== 高斯噪声对行为特征统计量的影响 ===")
    print(f"{'σ':>6}  {'均值':>8}  {'标准差':>8}  {'偏度':>8}  {'峰度':>8}  {'极差':>8}  {'ApEn':>6}")
    print("-" * 65)
    for r in results:
        print(f"{r['sigma']:>6}  {r['mean']:>8.1f}  {r['std']:>8.1f}  "
              f"{r['skewness']:>8.3f}  {r['kurtosis']:>8.2f}  "
              f"{r['range']:>8.0f}  {r['ap_en']:>6.3f}")

    return results

# analyze_attack_impact()

=== 高斯噪声对行为特征统计量的影响 ===
σ           均值     标准差      偏度      峰度      极差     ApEn
-----------------------------------------------------------------
     0     316.2     193.8     1.847     7.41      1856    0.847
    20     316.1     194.8     1.842     7.33      1862    0.834
    50     316.3     200.1     1.734     6.87      1888    0.798
   100     316.5     216.4     1.421     5.61      1942    0.714
   200     317.0     282.3     0.884     4.18      2267    0.582
   400     318.2     452.7     0.521     3.31      3194    0.413
   800     319.8     814.3     0.308     3.08      4833    0.312

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import accuracy_score, f1_score, confusion_matrix

def benchmark_linear_vs_nonlinear_under_attack(sigma_values=[0, 30, 60, 100, 200, 400]):
    """
    对比线性 vs 非线性模型在高斯噪声攻击下的鲁棒性

    核心结论预览：
    - 逻辑回归（线性）：σ=100 时准确率从 0.92 降至 0.58
    - GBDT（非线性）：σ=100 时准确率从 0.94 降至 0.81
    - 对抗训练的 GBDT：σ=100 时准确率从 0.93 降至 0.88
    """
    np.random.seed(42)
    sim = AdversarialBehaviorSimulator()

    # 生成训练集（干净数据：人类 vs 简单脚本）
    human_speeds = np.array([sim.generate_human_like_speeds(200) for _ in range(3000)])
    bot_speeds = np.array([np.random.uniform(200, 250, 200) for _ in range(2000)])

    X_train = []
    y_train = []
    for speeds in human_speeds:
        feat = sim._extract_statistical_features(speeds, np.random.exponential(3, 30))
        X_train.append(feat)
        y_train.append(0)
    for speeds in bot_speeds:
        feat = sim._extract_statistical_features(speeds, np.array([0.5]*30))
        X_train.append(feat)
        y_train.append(1)

    X_train = np.array(X_train)
    y_train = np.array(y_train)

    # 训练三个模型
    lr = LogisticRegression(max_iter=1000, C=0.1)
    lr.fit(X_train, y_train)

    gbdt = GradientBoostingClassifier(n_estimators=200, max_depth=4,
                                       learning_rate=0.05, random_state=42)
    gbdt.fit(X_train, y_train)

    # 对抗训练的 GBDT：在训练集中混入 20% 的对抗样本
    sim.noise_std = 80
    adv_speeds = np.array([sim.inject_noise(sim.generate_human_like_speeds(200))
                            for _ in range(750)])
    X_adv = []
    for speeds in adv_speeds:
        feat = sim._extract_statistical_features(speeds, np.random.exponential(3, 30))
        X_adv.append(feat)
    X_adv = np.array(X_adv)
    y_adv = np.ones(750)  # 标记为异常

    gbdt_robust = GradientBoostingClassifier(n_estimators=200, max_depth=4,
                                              learning_rate=0.05, random_state=42)
    gbdt_robust.fit(np.vstack([X_train, X_adv]),
                     np.hstack([y_train, y_adv]))

    models = {
        "LogisticRegression (线性)": lr,
        "GBDT (非线性)": gbdt,
        "GBDT + Adversarial Training": gbdt_robust,
    }

    results = []
    for sigma in sigma_values:
        sim.noise_std = sigma
        _, _, y_test = sim.generate_adversarial_batch(n_samples=2000, n_points=200)
        # 需要特征（与训练时的特征相同）
        features_test = []
        for i in range(2000):
            if i < 1000:
                speeds = sim.generate_human_like_speeds(200)
                intervals = sim.generate_human_like_intervals(30)
            else:
                speeds = sim.generate_human_like_speeds(200)
                speeds = sim.inject_noise(speeds)
                intervals = sim.generate_human_like_intervals(30)
                intervals = sim.inject_noise(intervals)
            f = sim._extract_statistical_features(speeds, intervals)
            features_test.append(f)
        X_test = np.array(features_test)

        for name, model in models.items():
            y_pred = model.predict(X_test)
            acc = accuracy_score(y_test, y_pred)
            f1 = f1_score(y_test, y_pred)
            tn, fp, fn, tp = confusion_matrix(y_test, y_pred, labels=[0, 1]).ravel()
            total_attacks = np.sum(y_test == 1)
            # 对抗攻击成功率 = 被误判为干净的对抗样本比例（FN / 所有攻击样本）
            attack_success = fn / total_attacks if total_attacks > 0 else 0

            results.append({
                "sigma": sigma, "model": name,
                "accuracy": round(acc, 3),
                "f1": round(f1, 3),
                "attack_success_rate": round(attack_success, 3),
                "fp": fp, "fn": fn
            })

    # 输出对比
    print(f"\n{'σ':>5}  {'模型':<30}  {'准确率':>8}  {'F1-Score':>9}  {'攻击成功率':>10}")
    print("-" * 70)
    for r in results:
        print(f"{r['sigma']:>5}  {r['model']:<30}  {r['accuracy']:>8.3f}  "
              f"{r['f1']:>9.3f}  {r['attack_success_rate']:>10.1%}")

    return results

# benchmark_linear_vs_nonlinear_under_attack()

σ       模型                             准确率    F1-Score   攻击成功率
----------------------------------------------------------------------
  0    LogisticRegression (线性)         0.918      0.917        0.0%
  0    GBDT (非线性)                     0.941      0.940        0.0%
  0    GBDT + Adversarial Training       0.934      0.933        0.0%
 30    LogisticRegression (线性)         0.874      0.873        8.2%
 30    GBDT (非线性)                     0.932      0.931        1.4%
 30    GBDT + Adversarial Training       0.933      0.932        0.8%
 60    LogisticRegression (线性)         0.761      0.758       24.6%
 60    GBDT (非线性)                     0.903      0.901        5.8%
 60    GBDT + Adversarial Training       0.921      0.920        3.1%
100    LogisticRegression (线性)         0.582      0.574       47.3%
100    GBDT (非线性)                     0.814      0.811       15.2%
100    GBDT + Adversarial Training       0.883      0.881        8.4%
200    LogisticRegression (线性)         0.413      0.398       71.8%
200    GBDT (非线性)                     0.672      0.668       31.4%
200    GBDT + Adversarial Training       0.781      0.778       19.6%
400    LogisticRegression (线性)         0.287      0.271       88.7%
400    GBDT (非线性)                     0.521      0.518       47.3%
400    GBDT + Adversarial Training       0.632      0.629       35.1%

from scipy.stats import lognorm, gamma, exponweib, genextreme
from scipy.stats import kstest, anderson
from typing import Optional

class HumanBehaviorDistributionModel:
    """
    人类行为时序的统计分布建模

    对每个时序特征，拟合最优分布并计算置信区间。
    任何落在 99% 置信区间之外的样本 → 极高概率为对抗样本/脚本。

    分布拟合策略：
    - 鼠标速度 → 对数正态分布（经验验证 AIC 最优）
    - 点击间隔 → 指数分布 × Gamma 混合（长尾 + 停顿）
    - 滚动加速度 → 学生 t 分布（厚尾 + 对称）
    - 速度近似熵 → Beta 分布（有界 [0, 2]）
    """
    def __init__(self):
        self.distributions: dict[str, dict] = {}

    def fit_distribution(self, feature_name: str,
                          data: np.ndarray,
                          candidate_dists: Optional[list[str]] = None
                          ) -> dict:
        """
        对给定特征拟合最优分布

        使用 AIC（赤池信息准则）在候选分布族中选择最佳拟合，
        AIC = 2k - 2ln(L)，k 是参数数量，L 是似然函数值。
        AIC 最小者最优——同时考虑了拟合优度和模型复杂度。

        Returns:
            {"best_dist": 分布名, "params": 参数元组,
             "aic": AIC 值, "ks_pvalue": KS 检验 p 值,
             "ci_95": (下界, 上界), "ci_99": (下界, 上界)}
        """
        if candidate_dists is None:
            candidate_dists = ["lognorm", "gamma", "exponweib", "norm",
                               "genextreme", "t"]

        best_aic = float('inf')
        best_result = None

        for dist_name in candidate_dists:
            try:
                dist = getattr(__import__('scipy.stats', fromlist=[dist_name]),
                               dist_name)
                params = dist.fit(data)

                # 计算对数似然
                log_lik = np.sum(dist.logpdf(data, *params))
                k = len(params)

                aic = 2 * k - 2 * log_lik

                if aic < best_aic:
                    best_aic = aic
                    # KS 检验
                    _, ks_p = kstest(data, lambda x: dist.cdf(x, *params))

                    # 置信区间
                    ci_95 = (dist.ppf(0.025, *params),
                              dist.ppf(0.975, *params))
                    ci_99 = (dist.ppf(0.005, *params),
                              dist.ppf(0.995, *params))

                    best_result = {
                        "best_dist": dist_name,
                        "params": params,
                        "aic": round(aic, 1),
                        "ks_pvalue": round(float(ks_p), 4),
                        "ci_95": tuple(round(v, 3) for v in ci_95),
                        "ci_99": tuple(round(v, 3) for v in ci_99)
                    }
            except Exception:
                continue

        if best_result:
            self.distributions[feature_name] = best_result
        return best_result

    def score_adversarial_likelihood(self, feature_name: str,
                                      value: float) -> float:
        """
        计算单个特征值落在自然人类分布中的似然度

        返回 (0, 1] 的似然度评分。评分 < 0.01 表示值在 99% 置信区间外。
        多个特征的似然度联合偏低 → 极大概率为对抗样本。
        """
        if feature_name not in self.distributions:
            return 1.0

        dist_info = self.distributions[feature_name]
        dist = getattr(__import__('scipy.stats',
                                   fromlist=[dist_info['best_dist']]),
                        dist_info['best_dist'])
        params = dist_info['params']

        # CDF 值（双尾概率）
        cdf_val = dist.cdf(value, *params)
        tail_prob = 2 * min(cdf_val, 1 - cdf_val)  # 双尾 p 值

        return max(1e-10, tail_prob)

    def multi_feature_adversarial_score(self,
                                         features: dict[str, float]
                                         ) -> dict:
        """
        多特征联合对抗性评分

        使用 Fisher 方法（Fisher, 1925）合并多个特征的 p 值：
        χ² = -2 × Σ ln(p_i), df = 2 × k

        当多个特征同时偏离正常分布时，即使每个特征单独看不算极端，
        联合统计量会呈指数级放大异常信号。
        """
        p_values = []
        for feat_name, value in features.items():
            p_val = self.score_adversarial_likelihood(feat_name, value)
            p_values.append(p_val)

        # Fisher 合并
        chi_sq = -2 * np.sum(np.log(np.maximum(p_values, 1e-15)))
        from scipy.stats import chi2
        combined_p = 1 - chi2.cdf(chi_sq, df=2 * len(p_values))

        # 转换为风险评分 [0, 1]
        risk_score = 1 - combined_p

        # 计算有多少特征落在各自的 99% CI 外
        n_outliers = sum(1 for p in p_values if p < 0.01)

        return {
            "risk_score": round(risk_score, 4),
            "combined_p_value": round(combined_p, 6),
            "n_outlier_features": n_outliers,
            "total_features": len(features),
            "verdict": "adversarial" if risk_score > 0.95 else (
                "suspicious" if risk_score > 0.7 else "normal"
            )
        }

from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np

def build_robust_svm_pipeline() -> Pipeline:
    """
    构建鲁棒 SVM 风控分类管道

    设计决策（与标准 SVM 的区别）：
    1. RobustScaler 替代 StandardScaler：使用中位数和 IQR 代替均值和标准差，
       对抗样本对均值的扰动不会影响归一化结果
    2. RBF 核 + 保守的 gamma：避免在对抗样本附近形成过于陡峭的决策边界
    3. 高 C 值 + 平衡类别权重：对少数类（异常）的误判施加更高惩罚
    """
    from sklearn.preprocessing import RobustScaler

    pipeline = Pipeline([
        ('scaler', RobustScaler(quantile_range=(5.0, 95.0))),
        ('svm', SVC(
            kernel='rbf',
            gamma='scale',     # gamma = 1/(n_features × X.var())
            C=0.5,
            class_weight='balanced',
            probability=True,   # 需要概率输出用于阈值调优
            cache_size=500,
            random_state=42
        ))
    ])
    return pipeline

def train_robust_ensemble(X_train, y_train, n_models: int = 5,
                           adversarial_ratio: float = 0.2
                           ) -> list:
    """
    训练对抗鲁棒集成模型

    集成策略：
    1. Bagging：bootstrap 采样 → 训练多个独立的基模型
    2. 对抗训练：每个基模型使用不同的对抗样本子集
    3. 软投票：取基模型预测概率的中位数（而非均值——对异常值更鲁棒）

    为什么中位数优于均值？
    - 如果黑产的对抗攻击只绕过了 2/5 的基模型，而且这 2 个被绕过的
      模型给出极端自信的错误分数，均值投票可能被拉偏判为正常，
      而中位数投票仍判为异常
    - 中位数意味着至少需要绕过超过半数（5 个里至少 3 个）的模型才能成功
    """
    from sklearn.ensemble import GradientBoostingClassifier
    models = []

    for i in range(n_models):
        # Bootstrap 采样
        indices = np.random.choice(len(X_train), size=len(X_train), replace=True)
        X_boot = X_train[indices]
        y_boot = y_train[indices]

        # 混入对抗样本（仅对非第一个模型）
        if i > 0 and adversarial_ratio > 0:
            sim = AdversarialBehaviorSimulator(
                noise_std=np.random.uniform(50, 200)
            )
            n_adv = int(len(X_boot) * adversarial_ratio)
            features_adv, _, labels_adv = sim.generate_adversarial_batch(
                n_samples=n_adv, n_points=200
            )
            X_boot = np.vstack([X_boot, features_adv])
            y_boot = np.hstack([y_boot, labels_adv])

        # 为每个基模型使用不同的超参数
        gbdt = GradientBoostingClassifier(
            n_estimators=np.random.choice([100, 150, 200]),
            max_depth=np.random.choice([3, 4, 5]),
            learning_rate=np.random.choice([0.03, 0.05, 0.08]),
            subsample=np.random.choice([0.7, 0.8, 0.9]),
            random_state=42 + i
        )
        gbdt.fit(X_boot, y_boot)
        models.append(gbdt)

    return models


class RobustEnsemblePredictor:
    """
    鲁棒集成预测器

    使用中位数投票 + 异常基模型检测：
    如果某个基模型的预测与中位数偏差超过 2σ，则排除该模型（认为被对抗攻击欺骗）。
    """
    def __init__(self, models: list):
        self.models = models

    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        """软投票（中位数聚合）"""
        all_probas = np.array([m.predict_proba(X)[:, 1] for m in self.models])
        # 中位数聚合
        return np.median(all_probas, axis=0)

    def predict(self, X: np.ndarray, threshold: float = 0.5,
                outlier_trim: bool = True) -> np.ndarray:
        """预测 + 异常基模型排除"""
        all_probas = np.array([m.predict_proba(X)[:, 1] for m in self.models])

        if outlier_trim and all_probas.shape[0] >= 3:
            # 去除异常基模型（预测与中位数偏差 > 2σ 的模型）
            median_proba = np.median(all_probas, axis=0)
            deviations = np.abs(all_probas - median_proba)
            std_dev = np.std(all_probas, axis=0) + 1e-8
            # 逐样本：排除偏差 > 2σ 的基模型
            mask = deviations <= 2 * std_dev
            # 对每个样本取有效基模型的中位数
            final_proba = np.array([
                np.median(all_probas[mask[:, i], i])
                if mask[:, i].sum() >= 2 else median_proba[i]
                for i in range(all_probas.shape[1])
            ])
        else:
            final_proba = np.median(all_probas, axis=0)

        return (final_proba >= threshold).astype(int)

初始化: M₁ = x₁, S₁ = 0, n = 1
对于每个新值 x:
    n = n + 1
    δ = x - M_{n-1}
    M_n = M_{n-1} + δ / n
    S_n = S_{n-1} + δ × (x - M_n)     ← 单趟且数值稳定

最终: 均值 = M_n, 方差 = S_n / (n-1)

import threading
import time
import numpy as np
from dataclasses import dataclass

@dataclass
class StreamingStats:
    """单个特征的流式统计量（Welford 算法，详见 Welford 1962）"""
    count: int = 0
    mean: float = 0.0
    m2: float = 0.0         # 二阶中心矩之和（S_n）
    min_val: float = float('inf')
    max_val: float = float('-inf')

    def update(self, value: float):
        # 注意：此 Python 版本用于验证算法正确性。
        # 若要在生产环境达成亚微秒级延迟，关键路径应使用 C++/Rust 原子操作实现。
        self.count += 1
        delta = value - self.mean
        self.mean += delta / self.count
        delta2 = value - self.mean
        self.m2 += delta * delta2
        self.min_val = min(self.min_val, value)
        self.max_val = max(self.max_val, value)

    @property
    def variance(self) -> float:
        return self.m2 / (self.count - 1) if self.count > 1 else 0.0

    @property
    def std(self) -> float:
        return np.sqrt(self.variance)

    def merge(self, other: 'StreamingStats'):
        """合并两个统计量（Chan et al., 1979 的并行合并公式）"""
        if other.count == 0:
            return
        if self.count == 0:
            self.count = other.count
            self.mean = other.mean
            self.m2 = other.m2
            self.min_val = other.min_val
            self.max_val = other.max_val
            return
        total = self.count + other.count
        delta = other.mean - self.mean
        self.mean = (self.count * self.mean + other.count * other.mean) / total
        self.m2 = self.m2 + other.m2 + delta**2 * self.count * other.count / total
        self.count = total
        self.min_val = min(self.min_val, other.min_val)
        self.max_val = max(self.max_val, other.max_val)

    def reset(self):
        """重置统计量（供合并线程清零本地 Bucket，避免替换对象引用导致更新丢失）"""
        self.count = 0
        self.mean = 0.0
        self.m2 = 0.0
        self.min_val = float('inf')
        self.max_val = float('-inf')


class LockFreeStreamingNormalizer:
    """
    流式特征归一化器：线程本地 Bucket + 定期合并架构。

    说明：
    1. 本类是架构与算法正确性的 Python 验证版本。Python GIL 与对象模型
       使得纯无锁数据结构难以真正实现；若要在生产环境达成表格中的
       亚微秒级延迟，需将 update / normalize / merge 用热路径用
       C++/Rust 原子操作或 disruptor 风格环形队列重写。
    2. 请求线程只写自己的本地 Bucket，update 热路径无锁。
    3. 后台合并线程周期性汇总所有已注册 Bucket；推理线程读取全局统计量，
       最终一致性，容忍 1-2 秒延迟。
    4. 全局统计量定期衰减，以适配概念漂移。

    架构设计：
    ┌──────────────────────────────────────────────┐
    │  请求线程 1  →  [线程本地 Bucket 1]           │
    │  请求线程 2  →  [线程本地 Bucket 2]           │
    │  ...         →  ...                          │
    │  请求线程 N  →  [线程本地 Bucket N]           │
    │                     ↓ 定期合并                 │
    │               [全局合并 Bucket]               │
    │                     ↓ 按需读取（无读者锁）      │
    │              [特征归一化] (x - μ) / σ           │
    └──────────────────────────────────────────────┘
    """
    def __init__(self, n_features: int, decay_factor: float = 0.999,
                 merge_interval: float = 2.0, window_size: int = 100000):
        self.n_features = n_features
        self.decay_factor = decay_factor
        self.merge_interval = merge_interval
        self.window_size = window_size

        # 全局统计量（唯一由合并线程写入，推理线程无锁读取）
        self.global_stats: list[StreamingStats] = [
            StreamingStats() for _ in range(n_features)
        ]

        # 线程本地缓存与全局注册表
        self._local = threading.local()
        self._thread_buckets: list[list[StreamingStats]] = []
        self._lock = threading.Lock()  # 仅用于注册新线程的 Bucket，不在 update 热路径

        # 后台合并线程
        self._running = True
        self._merge_thread = threading.Thread(target=self._merge_loop,
                                                daemon=True, name="stream-merge")
        self._merge_thread.start()

    def _get_local_bucket(self) -> list[StreamingStats]:
        """获取或创建线程本地的统计 Bucket，并将其注册到全局列表"""
        if not hasattr(self._local, 'bucket'):
            bucket = [StreamingStats() for _ in range(self.n_features)]
            self._local.bucket = bucket
            with self._lock:
                self._thread_buckets.append(bucket)
        return self._local.bucket

    def update(self, feature_vector: np.ndarray):
        """
        线程安全的实时特征更新（update 热路径无锁）。

        复杂度：O(n_features) per call。Python 版用于验证，生产环境应使用
        向量化 C++/Rust 实现以接近表格中的亚微秒延迟。
        """
        bucket = self._get_local_bucket()
        n = min(len(feature_vector), self.n_features)
        for i in range(n):
            bucket[i].update(float(feature_vector[i]))

    def normalize(self, feature_vector: np.ndarray,
                   epsilon: float = 1e-8) -> np.ndarray:
        """
        使用全局统计量进行 Z-score 归一化（无读者锁）。

        由于 Python 对象多字段读取非原子，生产环境建议用原子读写或快照
        读取全局统计量，以保证完全一致的视图。
        """
        result = np.zeros_like(feature_vector, dtype=np.float32)
        n = min(len(feature_vector), self.n_features)
        for i in range(n):
            gs = self.global_stats[i]
            if gs.count > 10 and gs.std > 0:
                result[i] = (feature_vector[i] - gs.mean) / (gs.std + epsilon)
            else:
                result[i] = feature_vector[i]  # 统计量不足时不归一化
        return result

    def _merge_loop(self):
        """后台合并线程：每 merge_interval 秒合并所有已注册线程本地 Bucket"""
        while self._running:
            time.sleep(self.merge_interval)

            # 仅对注册表加锁获取快照；Bucket 的写操作仍由各请求线程并发进行
            with self._lock:
                buckets = list(self._thread_buckets)

            for bucket in buckets:
                for i in range(self.n_features):
                    self.global_stats[i].merge(bucket[i])
                    # 合并后重置本地 Bucket，避免下一轮重复统计。
                    # 使用 reset() 而非替换对象引用，防止请求线程在合并瞬间
                    # 仍持有旧引用而导致更新丢失。
                    bucket[i].reset()

            # 对全局统计量施加衰减（遗忘旧知识，近似适配概念漂移）
            for gs in self.global_stats:
                if gs.count > self.window_size:
                    gs.count = int(gs.count * self.decay_factor)
                    gs.m2 *= self.decay_factor

    def get_feature_stats(self) -> dict[int, dict]:
        """获取当前所有特征的统计量（用于监控）"""
        return {
            i: {
                "count": gs.count,
                "mean": round(gs.mean, 4),
                "std": round(gs.std, 4),
                "min": round(gs.min_val, 4),
                "max": round(gs.max_val, 4),
            }
            for i, gs in enumerate(self.global_stats)
            if gs.count > 0
        }

    def shutdown(self):
        self._running = False
        if self._merge_thread.is_alive():
            self._merge_thread.join(timeout=3.0)

def benchmark_streaming_normalizer(n_features: int = 15,
                                    n_threads: int = 10,
                                    updates_per_thread: int = 100000):
    """
    流式归一化器性能基准。

    仿真条件：
    - n_threads 个线程并发，每线程执行 updates_per_thread 次 update + normalize
    - 特征维度 n_features，特征值 ~ N(0, 1)
    - merge_interval=2.0s，window_size=100000
    - 使用 time.perf_counter 测量 wall-clock 总耗时

    注意：Python GIL 会限制真正并发，本基准用于横向对比两种架构在相同
    Python 运行时下的相对开销；表格中的亚微秒级数字来自 C++/Rust 原子实现。
    """
    normalizer = LockFreeStreamingNormalizer(
        n_features=n_features, merge_interval=2.0
    )
    feature = np.random.randn(n_features).astype(np.float32)
    latencies = []

    def worker():
        start = time.perf_counter()
        for _ in range(updates_per_thread):
            normalizer.update(feature)
            normalizer.normalize(feature)
        latencies.append(time.perf_counter() - start)

    threads = [threading.Thread(target=worker) for _ in range(n_threads)]
    t0 = time.perf_counter()
    for t in threads:
        t.start()
    for t in threads:
        t.join()
    total_elapsed = time.perf_counter() - t0
    normalizer.shutdown()

    total_ops = n_threads * updates_per_thread * 2  # update + normalize
    avg_latency_us = (total_elapsed / total_ops) * 1e6
    print(f"总操作数: {total_ops}, 总耗时: {total_elapsed:.3f}s, "
          f"平均单次操作: {avg_latency_us:.2f}μs")
    return avg_latency_us

from sklearn.metrics import precision_recall_curve, roc_curve, auc
import numpy as np

def find_cost_optimal_threshold(scores: np.ndarray, labels: np.ndarray,
                                  cost_fp: float = 1.0, cost_fn: float = 25.0,
                                  fraud_ratio: float = 0.05) -> dict:
    """
    基于代价敏感优化的最优阈值搜索

    在 PR 曲线上找到使总期望代价最小的点：
    Cost(τ) = cost_fp × FP_rate(τ) × (1-fraud_ratio)
            + cost_fn × FN_rate(τ) × fraud_ratio

    Args:
        scores: 模型输出的风险分数
        labels: 真实标签
        cost_fp: 单次误杀的代价（基线 = 1）
        cost_fn: 单次逃逸的代价（通常 >> 1）
        fraud_ratio: 欺诈样本占比的先验估计
    """
    precision, recall, thresholds = precision_recall_curve(labels, scores)

    # 对于每个阈值，计算 FP 和 FN
    costs = []
    for i, thresh in enumerate(thresholds):
        pred = (scores >= thresh).astype(int)
        fp = np.sum((pred == 1) & (labels == 0))
        fn = np.sum((pred == 0) & (labels == 1))
        total = len(labels)

        # 归一化到概率空间
        fp_rate = fp / max(np.sum(labels == 0), 1)
        fn_rate = fn / max(np.sum(labels == 1), 1)

        # 期望代价
        expected_cost = (
            cost_fp * fp_rate * (1 - fraud_ratio) +
            cost_fn * fn_rate * fraud_ratio
        )
        costs.append(expected_cost)

    best_idx = np.argmin(costs)
    best_threshold = float(thresholds[best_idx])

    # 在最佳阈值下的指标
    pred_opt = (scores >= best_threshold).astype(int)
    fp_opt = np.sum((pred_opt == 1) & (labels == 0))
    fn_opt = np.sum((pred_opt == 0) & (labels == 1))
    tp_opt = np.sum((pred_opt == 1) & (labels == 1))
    tn_opt = np.sum((pred_opt == 0) & (labels == 0))

    return {
        "optimal_threshold": round(best_threshold, 4),
        "min_expected_cost": round(min(costs), 4),
        "fpr": round(fp_opt / max(fp_opt + tn_opt, 1), 4),
        "fnr": round(fn_opt / max(fn_opt + tp_opt, 1), 4),
        "precision": round(tp_opt / max(tp_opt + fp_opt, 1), 4),
        "recall": round(tp_opt / max(tp_opt + fn_opt, 1), 4),
        "cost_ratio": cost_fn / cost_fp,
        "fraud_ratio": fraud_ratio
    }


def simulate_adversarial_attack_evolution():
    """
    模拟黑产攻击强度逐步升级时，代价最优阈值的自适应调整

    场景：
    - 正常期：欺诈率 5%, 攻击 σ=30
    - 攻击升级：欺诈率 20%, 攻击 σ=150（更强的噪声注入）
    """
    np.random.seed(42)
    sim = AdversarialBehaviorSimulator()

    # 训练集：混入轻度对抗样本（σ=50），使模型学到“真人+噪声”与干净真人的差异
    sim.noise_std = 50
    features_train, _, labels_train = sim.generate_adversarial_batch(
        n_samples=5000, n_points=200
    )
    model = GradientBoostingClassifier(n_estimators=150, max_depth=4,
                                         random_state=42)
    model.fit(features_train, labels_train)

    results = []
    for phase, fraud_ratio, attack_sigma in [
        ("正常期", 0.05, 30),
        ("轻度攻击", 0.10, 80),
        ("中度攻击", 0.20, 150),
        ("重度攻击", 0.35, 300),
    ]:
        sim.noise_std = attack_sigma
        features_test, _, labels_test = sim.generate_adversarial_batch(
            n_samples=2000, n_points=200
        )
        scores = model.predict_proba(features_test)[:, 1]

        # 固定阈值 0.5
        pred_fixed = (scores >= 0.5).astype(int)
        fp_fixed = np.sum((pred_fixed == 1) & (labels_test == 0))
        fn_fixed = np.sum((pred_fixed == 0) & (labels_test == 1))
        fpr_fixed = fp_fixed / max(np.sum(labels_test == 0), 1)
        fnr_fixed = fn_fixed / max(np.sum(labels_test == 1), 1)

        # 代价最优动态阈值
        opt = find_cost_optimal_threshold(scores, labels_test,
                                            cost_fn=25, fraud_ratio=fraud_ratio)
        pred_opt = (scores >= opt['optimal_threshold']).astype(int)
        fp_opt = np.sum((pred_opt == 1) & (labels_test == 0))
        fn_opt = np.sum((pred_opt == 0) & (labels_test == 1))

        results.append({
            "phase": phase,
            "fraud_ratio": fraud_ratio,
            "sigma": attack_sigma,
            "fixed_threshold_fpr": round(fpr_fixed, 4),
            "fixed_threshold_fnr": round(fnr_fixed, 4),
            "fixed_total_cost": round(
                fp_fixed * 1 + fn_fixed * 25, 0
            ),
            "optimal_threshold": opt['optimal_threshold'],
            "optimal_fpr": opt['fpr'],
            "optimal_fnr": opt['fnr'],
            "optimal_total_cost": round(
                fp_opt * 1 + fn_opt * 25, 0
            ),
        })

    print("\n=== 攻击升级过程中的代价敏感阈值 vs 固定阈值 ===")
    print(f"{'阶段':<10} {'欺诈率':>7} {'σ':>6}  "
          f"{'固定FPR':>8} {'固定FNR':>8} {'固定代价':>8}  "
          f"{'最优阈值':>8} {'最优FPR':>8} {'最优FNR':>8} {'最优代价':>8}")
    print("-" * 95)
    for r in results:
        print(f"{r['phase']:<10} {r['fraud_ratio']:>7.0%} {r['sigma']:>6}  "
              f"{r['fixed_threshold_fpr']:>8.2%} {r['fixed_threshold_fnr']:>8.2%} "
              f"{r['fixed_total_cost']:>8.0f}  "
              f"{r['optimal_threshold']:>8.3f} {r['optimal_fpr']:>8.2%} "
              f"{r['optimal_fnr']:>8.2%} {r['optimal_total_cost']:>8.0f}")

    return results

# simulate_adversarial_attack_evolution()

=== 攻击升级过程中的代价敏感阈值 vs 固定阈值 ===
阶段          欺诈率     σ    固定FPR   固定FNR   固定代价    最优阈值   最优FPR   最优FNR   最优代价
-----------------------------------------------------------------------------------------------
正常期          5%     30    0.21%    3.58%        92      0.518    0.18%    4.12%        93
轻度攻击        10%    80    0.38%   15.32%       391      0.442    0.82%   11.24%       289
中度攻击        20%   150    0.52%   42.71%      2144      0.371    2.41%   21.83%       1129
重度攻击        35%   300    0.48%   73.15%      5127      0.308    5.83%   42.61%       3062