Sklearn-源码解析-书-v1-0-三-

Sklearn 源码解析（书）v1.0（三）

解释：MetadataRouter 维护一个 树形结构：根节点（元估计器）拥有子节点（子估计器），每个节点保存 请求 与 方法映射。在 process_routing 时根据该树生成实际传递给子对象的 kwargs。

生活类比：这就像一个家庭树——祖父母（根节点）掌握全家的资源分配规则（如谁可以使用车辆、谁负责买菜），而每个子女和孙辈（子节点）则根据这些规则申请特定资源。当某个孩子想用车时，他需要先按照家庭规则提出申请，祖父母根据树中的映射判断是否批准以及应该给哪辆车。

8.8.4 流程图：元数据路由全链路

flowchart TD A[MetaEstimator] -->|set_*_request| B[MetadataRouter] B -->|add_self_request| C[Self Request] B -->|add| D[Child Estimator 1] B -->|add| E[Child Estimator 2] A -->|fit(X, sample_weight=s)| F[process_routing] F -->|拆解请求| G[子 Estimator 1: fit(sample_weight=s)] F -->|拆解请求| H[子 Estimator 2: fit(sample_weight=s)] G -->|执行| I[子 Estimator 1 实际 fit] H -->|执行| J[子 Estimator 2 实际 fit]

---

8.9 文档字符串质量防线（src/sklearn/tests/test_docstring_parameters*.py、test_docstrings.py）

8.9.1 参数‑文档一致性检查（test_docstring_parameters.py）

源码路径：sklearn/tests/test_docstring_parameters.py

def test_docstring_parameters():
    for name in PUBLIC_MODULES:
        module = importlib.import_module(name)
        for cls_name, cls in inspect.getmembers(module, inspect.isclass):
            if cls.__module__.startswith("sklearn"):
                cdoc = docscrape.ClassDoc(cls)                     # ← 解析类文档
                errors = check_docstring_parameters(cls.__init__, cdoc)
                for method_name in cdoc.methods:                    # ← 检查每个方法
                    method = getattr(cls, method_name)
                    errors += check_docstring_parameters(method)
            # 同理检查函数

解释：check_docstring_parameters 对比 函数签名 与 numpydoc 参数段，确保每个参数都有准确描述、类型/默认值一致。错误码（如 RT02、GL01）在 test_docstrings.py 中进一步过滤，保持 友好的报错。

生活类比：这就像食品包装上的营养成分表——如果标签上写着“每100克含糖量5克”，但实际检测发现是8克，消费者就会感到被欺骗。同样，如果函数签名说接受一个 float 参数，但文档却写成了“整数”，用户在使用时可能会因为类型错误而调用失败。

8.9.2 类/函数文档一致性（test_docstring_parameters_consistency.py）

源码路径：sklearn/tests/test_docstring_parameters_consistency.py

def test_class_docstring_consistency(case):
    # 对同类家族（如 Bagging*）的公共参数进行交叉检查
    assert_docstring_consistency(**case)   # ← 统一正则匹配描述

解释：此测试确保 关联类（如 BaggingClassifier 与 BaggingRegressor）共享的参数（max_samples）文档保持 统一措辞，防止文档漂移。

生活类比：想象一个连锁超市对“牛奶”的描述——在北京店写的是“来自华北地区的鲜牛奶”，而在广州店却写成了“进口德牛奶”。如果顾客在不同城市看到矛盾的描述，会对品牌信任度产生怀疑。保持描述一致就如同确保所有分店都使用统一的供应链和标准描述。

8.9.3 全局 docstring 验证（test_docstrings.py）

源码路径：sklearn/tests/test_docstrings.py

def test_docstring():
    for Klass, method in get_all_methods():
        import_path = f"{Klass.__module__}.{Klass.__name__}"
        if method:
            import_path += f".{method}"
        res = numdoccpydoc.validate.validate(import_path)
        res["errors"] = list(filter_errors(res["errors"], method, Klass))
        if res["errors"]:
            raise ValueError(repr_errors(res, Klass, method))

解释：遍历所有公开类与方法，对 每个导入路径 进行 numpydoc.validate，并依据项目特定规则（忽略 RT02、GL01…）过滤已知不适用的错误。

生活类比：这就像质量检查员不仅要检查单个产品的标签是否正确，还要定期抽查整个仓库的货物，确保没有因批次更换或包装线故障导致的系统性错误。只有当所有产品都通过抽检，整个批次才能被视为合格。

8.9.4 流程图：文档校验管线

flowchart LR A[源码库] --> B[收集所有公开类/函数] B --> C[numpydoc.validate] C --> D{过滤已知错误码} D -->|通过| E[测试通过] D -->|失败| F[抛出 AssertionError]

---

8.10 配置线程安全与 OpenMP 检测（src/sklearn/tests/test_config.py、test_build.py）

8.10.1 配置上下文（test_config.py）

源码路径：sklearn/tests/test_config.py

def test_config_context():
    assert get_config() == {..."assume_finite": False,...}
    with config_context(assume_finite=True):
        assert get_config()["assume_finite"] is True      # ← 临时生效
    # 退出上下文后恢复原值
    assert get_config()["assume_finite"] is False

解释：config_context 使用 线程局部存储（threading.local），确保在 并行任务（Joblib、ThreadPoolExecutor）中不同线程的配置互不干扰。

生活类比：想象一个共享办公空间，每个人都可以调节自己工位的台灯亮度而不影响他人——一个人把灯调成阅读模式，另一个人却可以同时把灯调成会议模式，互不干扰。这种局部控制正是线程局部存储的体现。

8.10.2 线程安全验证（test_config_threadsafe_joblib、test_config_threadsafe）

源码路径：sklearn/tests/test_config.py

def set_assume_finite(assume_finite, sleep_duration):
    with config_context(assume_finite=assume_finite):
        time.sleep(sleep_duration)
        return get_config()["assume_finite"]
# 第 8 章 —— 并行调用后返回列表应保持调用顺序对应的值

解释：两个作业交叉执行，分别设定 assume_finite=False/True，最终返回值 不受对方影响，证明 全局配置是线程安全的。

生活类比：就像两个厨师在同一间厨房里做菜——一个在做需要高温快速翻炒的菜（对应 assume_finite=True），另一个在做需要慢火炖煮的汤（对应 assume_finite=False）。只要厨房的燃气阀和温度控制是独立的，两人就不会因为对方的操作而影响自己的菜肴。

8.10.3 OpenMP 编译检测（test_build.py）

源码路径：sklearn/tests/test_build.py

def test_openmp_parallelism_enabled():
    assert _openmp_parallelism_enabled(), """
        This test fails because scikit‑learn has been built without OpenMP.
        …
    """

解释：若编译时未链接 OpenMP，_openmp_parallelism_enabled() 返回 False，触发明确的 构建警告，提醒维护者在 高级安装文档 中提供正确的编译指令。

生活类比：这就像买了一辆标称为“四轮驱动”的越野车，但实际发现只有两轮带动力——如果不及时发现这个问题，在泥泞路面上可能会陷车。制造商应在交付前进行充分测试，并在使用手册中明确说明实际性能，以避免用户产生错误预期。

8.10.4 流程图：配置 → 并行执行 → 恢复

sequenceDiagram participant Main participant Worker1 participant Worker2 Main->>Worker1: config_context(assume_finite=False) Main->>Worker2: config_context(assume_finite=True) Worker1-->>Main: 返回 False Worker2-->>Main: 返回 True Note over Main: 配置在每个线程独立保存

---

8.11 设计取舍（一问一答）

Q1：为何不直接使用 @dataclass 来存储元数据请求？

A1：@dataclass 会在实例化时自动生成 __init__，但 元数据路由 需要 细粒度的运行时修改（如 set_fit_request(sample_weight="alias")、在 fit_transform 时自动合并请求）。手写的 MetadataRequest 与 MetadataRouter 能够 在每次调用前动态更新请求树，而 dataclass 难以实现这种 按需别名映射 与 自请求拷贝 的行为。

生活类比：这就像你不能用一个预先印好的固定菜单来应对每日变化的厨房需求——今天可能需要临时添加一道素食选项，明天又要根据过敏情况删除某种食材。只有手写的、可以随时修改的点菜系统才能灵活应对这种变化。

Q2：MetadataRouter 引入的复杂性是否值得？

A2：收益在于 显式、可验证的元数据流：每一次 fit、predict、transform 等都经过 统一路由检查，能够在 错误路径（未声明的 sample_weight）提前抛出 UnsetMetadataPassedError，防止模型在生产环境中因隐式参数遗漏而产生不可复现的行为。相对的，缺点是 API 增加 set_*_request 方法，使用者需要学习新概念。

生活类比：这相当于在工厂引入一个电子追踪系统来监控每个零件的来源和去向——虽然最初需要培训员工和维护系统，但一旦上线，就能快速定位哪批次产品因使用了错误的螺丝而导致故障，从而避免大规模召车。相对地，如果完全依赖人工目视检查，虽然看似简单，但一旦出问题就很难快速定位根源。

---

8.12 小结

本章通过 源码路径定位、逐行注释、后置解释 的方式，完整展现了 scikit‑learn 从 项目愿景、依赖管理、子包组织 到 BaseEstimator 克隆、元估计器委托、元数据路由、文档校验、配置线程安全 的全链路工程。

• README 的徽章与依赖声明为项目第一印象提供质量保证。

• 子包 __init__ 通过 re‑export 隐藏实现细节，保持 API 稳定。

• BaseEstimator.clone 与 __sklearn_clone__ 确保对象复制的可控性。

• 元估计器 通过 方法委托 与 数据验证下放 实现组合模型的透明行为。

• MetadataRouter 与 process_routing 构筑 元数据请求‑消费 框架，所有 sample_weight / metadata 必须显式声明。

• numpydoc 三级校验（参数约束、类/函数一致性、全局验证）守护文档质量。

• config_context 与 OpenMP 检测 提供运行时与编译时的安全保障。

后续阅读：第 9 章将深入 线性模型（ElasticNet、Ridge、LogisticRegression 等）的数值求解器与正则化路径，实现细节与实验指南。

8.13 生活类比

想象 scikit-learn 的工程体系是一座精密运作的医药研发中心： README 徽章矩阵 = 大楼门口的认证标牌（ISO质量认证、GMP药品生产规范），一眼展示可靠性 最小依赖版本声明 = 配方中的原料纯度要求（Python 3.11 = 纯度≥99.9%的试剂），低于标准立即拒收 子包 init.py = 各楼层的科室导览牌（linear_model 楼 = 骨科，svm 楼 = 心内科） test_common.py 全量体检 = 药品上市前的 III 期临床试验，每个估计器都必须通过数百项检查 克隆语义测试 = 药品留样制度（克隆品与原件成分完全相同，但独立存储、互不影响） 元数据路由测试 = 医院内部的标本流转系统（sample_weight 标签必须准确送达对应科室） 文档字符串验证 = 药品说明书审查（每个参数说明必须与配方成分精确对应） 就像医药研发需要从原料、生产、质检到说明书的全链条管控，scikit-learn 从 README 承诺、模块组织、克隆语义到测试套件构筑了完整的质量防线。

8.14 模块地图/架构图

README.md
├── 徽章矩阵配置（1-35行）
│   ├── Azure/CircleCI/Codecov 状态徽章
│   ├── PythonVersion/PyPI/DOI 元数据徽章
│   └── Benchmark 性能追踪徽章
├── 最小依赖版本声明（36-46行）
│   ├── |PythonMinVersion| replace:: 3.11
│   ├── |NumPyMinVersion| replace:: 1.24.1
│   └── |SciPyMinVersion| replace:: 1.10.0
├── 项目历史与定位（48-60行）
│   ├── 2007年 GSoC 项目起源
│   └── BSD 3-Clause 协议与志愿者社区
├── 安装指南（62-95行）
│   ├── 核心依赖与可选依赖
│   └── pip/conda 安装路径
└── 开发与贡献指南（97-160行）
    ├── git clone 快速上手
    ├── pytest sklearn 测试命令
    └── SKLEARN_SEED 环境变量
sklearn/linear_model/__init__.py
├── 模块 docstring（1行）
├── 13个私有模块导入（8-48行）
│   ├── _base: LinearRegression
│   ├── _coordinate_descent: ElasticNet/Lasso 家族
│   ├── _logistic: LogisticRegression(CV)
│   ├── _ridge: Ridge/RidgeCV/RidgeClassifier
│   └── _stochastic_gradient: SGDClassifier/Regressor
└── __all__ 白名单（50-90行）
sklearn/svm/__init__.py
├── _bounds: l1_min_c 边界工具
├── _classes: SVC/SVR/LinearSVC/OneClassSVM
└── __all__ 导出 8 个公共符号
sklearn/neural_network/__init__.py
├── _multilayer_perceptron: MLPClassifier/MLPRegressor
└── _rbm: BernoulliRBM
sklearn/feature_selection/__init__.py
├── _base: SelectorMixin 统一接口
├── _from_model: SelectFromModel
├── _rfe: RFE/RFECV 递归特征消除
├── _sequential: SequentialFeatureSelector
└── _univariate_selection: chi2/f_classif 等统计检验
sklearn/semi_supervised/__init__.py
├── _label_propagation: LabelPropagation/LabelSpreading
└── _self_training: SelfTrainingClassifier
sklearn/tests/test_base.py
├── 测试辅助类定义（30-95行）
│   ├── MyEstimator/K/T/ModifyInitParams 基础估计器
│   ├── Buggy/NoEstimator/VargEstimator 反模式
│   ├── NaNTag/NoNaNTag/OverrideTag/DiamondOverwriteTag/InheritDiamondOverwriteTag 标签系统
│   └── TreeBadVersion/TreeNoVersion/DontPickleAttributeMixin/SingleInheritanceEstimator 序列化辅助
├── clone 系列（100-430行）
│   ├── test_clone: 深度复制基本语义
│   ├── test_clone_2: 不复制自定义属性
│   ├── test_clone_buggy: Buggy/NoEstimator/VargEstimator/ModifyInitParams 反模式
│   ├── test_clone_empty_array/nan/dict/sparse_matrices/estimator_types: 特殊类型克隆
│   ├── test_clone_class_rather_than_instance: 类而非实例的报错
│   ├── test_clone_pandas_dataframe: DummyEstimator 克隆
│   ├── test_clone_protocol: FrozenEstimator.__sklearn_clone__ 协议
│   └── test_clone_keeps_output_config: 输出配置保留
├── 序列化版本追踪（210-448行）
│   ├── test_pickle_version_warning_is_not_raised_with_matching_version
│   ├── TreeBadVersion/TreeNoVersion 辅助类
│   ├── test_pickle_version_warning_is_issued_upon_different_version
│   ├── test_pickle_version_warning_is_issued_when_no_version_info_in_pickle
│   └── test_pickle_version_no_warning_is_issued_with_non_sklearn_estimator
├── getstate/setstate 重写（348-448行）
│   ├── DontPickleAttributeMixin.__getstate__/__setstate__
│   ├── MultiInheritanceEstimator/SingleInheritanceEstimator 辅助类
│   └── test_pickling_when_getstate_is_overwritten_by_mixin/..._by_mixin_outside_of_sklearn/...in_the_child_class
├── 类型判断（120-188行）
│   ├── test_is_classifier: SVC/GridSearchCV/Pipeline 穿透
│   ├── test_is_regressor: SVR/GridSearchCV/Pipeline 穿透
│   └── test_is_clusterer: KMeans/GridSearchCV/Pipeline 穿透
├── get_params/set_params（100-275行）
│   ├── test_get_params: deep=True/False 控制嵌套参数
│   ├── test_set_params: 双下划线语法与错误路径
│   ├── test_set_params_passes_all_parameters: 所有参数一起传递
│   ├── test_set_params_updates_valid_params: 更新合法参数
│   └── test_raises_on_get_params_non_attribute: 非属性参数报错
├── tag_inheritance（300-315行）
│   └── test_tag_inheritance: 菱形继承 MRO 标签合并规则
├── 输入验证（440-530行）
│   ├── test_n_features_in_validation/test_n_features_in_no_validation: _check_n_features
│   ├── test_feature_names_in/test_validate_data_skip_check_array/test_dataframe_protocol: validate_data 与特征名
│   └── NoOpTransformer/NoOpTransformer 辅助类
├── 表示层测试（485-540行）
│   ├── test_repr/test_str: 打印表示
│   ├── test_conditional_attrs_not_in_dir: __dir__ 条件属性
│   ├── test_repr_mimebundle_/test_repr_html_wraps: HTML 表示
│   └── make_estimator_with_param/DynamicEstimator: 动态估计器
├── 参数 HTML（525-615行）
│   ├── test_get_params_html: _get_params_html
│   ├── test_param_is_non_default/test_param_is_non_default_when_pandas_NA/test_param_is_default: 默认值检测
│   └── make_estimator_with_param/DynamicEstimator: 动态估计器模板
├── 元数据路由警告（480-530行）
│   ├── test_transformer_fit_transform_with_metadata_in_transform: CustomTransformer 辅助类
│   └── test_outlier_mixin_fit_predict_with_metadata_in_predict: CustomOutlierDetector 辅助类
└── 特殊方法测试（360-370行）
    └── test_score_sample_weight: 加权评分
sklearn/tests/test_common.py
├── test_all_estimator_no_base_class（44-51行）
├── test_get_check_estimator_ids（65-85行）
│   └── _sample_func/CallableEstimator 辅助
├── test_estimators（50-60行）
│   └── parametrize_with_checks 驱动全量测试
├── test_import_all_consistency（65-80行）
│   └── __all__ 中每个名字必须真实存在
├── test_root_import_all_completeness（82-95行）
│   └── 根包 __all__ 覆盖所有公开子模块
├── test_all_tests_are_importable（98-125行）
│   └── 每个子包必须有 tests 子包
├── test_class_support_removed（127-138行）
├── test_pandas_column_name_consistency（145-175行）
├── test_transformers_get_feature_names_out（190-205行）
│   └── _include_in_get_feature_names_out_check 过滤
├── test_estimators_get_feature_names_out_error（210-220行）
├── test_check_param_validation（225-235行）
├── test_set_output_transform（150-165行）
├── test_set_output_transform_configured（240-265行）
├── test_check_inplace_ensure_writeable（270-300行）
└── test_check_all_zero_sample_weights_error（305-320行）
sklearn/tests/test_public_functions.py
├── _get_func_info（10-30行）
├── _check_function_param_validation（30-85行）
│   ├── BadType 强制触发 InvalidParameterError
│   └── generate_valid/invalid_param_val 自动边界生成
├── PARAM_VALIDATION_FUNCTION_LIST（90-230行）
│   └── 200+ 公开函数的约束元数据
├── test_function_param_validation（230-245行）
└── test_class_wrapper_param_validation（270-290行）
sklearn/tests/test_metaestimators.py
├── DelegatorData 配置类（30-55行）
├── DELEGATING_METAESTIMATORS 列表（55-90行）
├── test_metaestimator_delegation（55-100行）
│   ├── SubEstimator @hides 动态方法隐藏
│   ├── SubEstimator.__init__/fit/_check_fit/inverse_transform/transform/predict/predict_proba/predict_log_proba/decision_function/score
│   └── 拟合前后委托调用验证
├── _get_instance_with_pipeline（170-210行）
├── _generate_meta_estimator_instances_with_pipeline（213-240行）
├── DATA_VALIDATION_META_ESTIMATORS_TO_IGNORE/DATA_VALIDATION_META_ESTIMATORS 列表
└── test_meta_estimators_delegate_data_validation（170-200行）
sklearn/tests/test_metaestimators_metadata_routing.py
├── METAESTIMATORS 列表（60-280行）
│   ├── estimator_name/routing_methods/preserves_metadata
│   └── method_mapping 描述 caller→callee 映射
├── get_init_args（300-370行）
├── filter_metadata_in_routing_methods（373-400行）
├── set_requests（403-440行）
├── test_unsupported_estimators_get_metadata_routing（443-450行）
├── test_unsupported_estimators_fit_with_metadata（453-465行）
├── test_registry_copy（468-475行）
├── test_default_request（478-492行）
├── test_error_on_missing_requests_for_sub_estimator（190-230行）
├── test_setting_request_on_sub_estimator_removes_error（235-280行）
├── test_non_consuming_estimator_works（495-530行）
├── test_metadata_is_routed_correctly_to_scorer（300-330行）
├── test_metadata_is_routed_correctly_to_splitter（335-360行）
└── test_metadata_routed_to_group_splitter（363-390行）
sklearn/tests/metadata_routing_common.py
├── record_metadata/record_metadata_not_default（25-44行）
├── check_recorded_metadata（47-85行）
├── assert_request_is_empty/assert_request_equal（88-122行）
├── _Registry: 追踪克隆后丢失引用的实例（125-134行）
├── ConsumingRegressor（150-188行）
│   ├── __init__/partial_fit/fit/predict/score
├── NonConsumingClassifier（192-225行）
│   ├── __init__/fit/partial_fit/decision_function/predict/predict_proba/predict_log_proba
├── NonConsumingRegressor（228-237行）
├── ConsumingClassifier（255-325行）
│   ├── __init__/partial_fit/fit/predict/predict_proba/predict_log_proba/decision_function/score
├── ConsumingClassifierWithoutPredictProba（330-340行）
├── ConsumingClassifierWithoutPredictLogProba（345-355行）
├── ConsumingClassifierWithOnlyPredict（360-372行）
├── ConsumingTransformer（370-408行）
│   ├── __init__/fit/transform/fit_transform/inverse_transform
├── ConsumingNoFitTransformTransformer（411-426行）
├── ConsumingScorer（429-444行）
│   ├── __init__/_score
├── ConsumingSplitter（447-473行）
│   ├── __init__/split/get_n_splits/_iter_test_indices
├── ConsumingSplitterInheritingFromGroupKFold（476-478行）
├── MetaRegressor（479-494行）
│   ├── __init__/fit/get_metadata_routing
├── WeightedMetaRegressor（497-530行）
│   ├── __init__/fit/predict/get_metadata_routing
├── WeightedMetaClassifier（533-563行）
│   ├── __init__/fit/get_metadata_routing
└── MetaTransformer（566-590行）
    ├── __init__/fit/transform/get_metadata_routing
sklearn/tests/test_metadata_routing.py
├── SimplePipeline（48-100行）
│   ├── __init__/fit/predict/get_metadata_routing
├── test_assert_request_is_empty（60-90行）
├── test_estimator_puts_self_in_registry（60-65行）
├── test_request_type_is_alias（95-110行）
├── test_request_type_is_valid（113-130行）
├── test_default_requests（115-140行）
│   └── OddEstimator.fit
├── test_default_request_override（160-185行）
│   └── Base/class_1/Class_1.split
├── test_process_routing_invalid_method（193-195行）
├── test_process_routing_invalid_object（80-90行）
├── test_process_routing_empty_params_get_with_default（198-210行）
├── test_simple_metadata_routing（195-230行）
├── test_nested_routing（100-140行）
├── test_nested_routing_conflict（140-165行）
├── test_invalid_metadata（250-275行）
├── test_get_metadata_routing（250-275行）
│   └── TestDefaults.fit/score/predict
├── test_setting_default_requests（350-390行）
│   ├── ExplicitRequest/ExplicitRequestOverwrite/ImplicitRequest/ImplicitRequestRemoval.fit
├── test_removing_non_existing_param_raises（175-190行）
├── test_method_metadata_request（210-230行）
├── test_get_routing_for_object（250-275行）
│   └── Consumer.fit
├── test_metadata_request_consumes_method（280-300行）
├── test_metadata_router_consumes_method（305-330行）
├── test_metaestimator_warnings（345-360行）
├── test_estimator_warnings（365-380行）
├── test_string_representations（400-430行）
├── test_validations（435-460行）
├── test_methodmapping（465-485行）
├── test_metadatarouter_add_self_request（250-270行）
├── test_metadata_routing_add（510-545行）
├── test_metadata_routing_get_param_names（550-585行）
├── test_method_generation（590-630行）
│   ├── SimpleEstimator 第一版：无 sample_weight 参数的方法集
│   └── SimpleEstimator 第二版：有 sample_weight 参数的方法集
├── test_composite_methods（310-340行）
├── test_no_feature_flag_raises_error（640-650行）
├── test_none_metadata_passed（655-660行）
├── test_no_metadata_always_works（665-685行）
├── test_unsetmetadatapassederror_correct（690-710行）
├── test_unsetmetadatapassederror_correct_for_composite_methods（715-740行）
└── test_unbound_set_methods_work（745-770行）
sklearn/tests/test_docstring_parameters.py
├── PUBLIC_MODULES 发现（30-45行）
├── test_docstring_parameters（50-95行）
├── _construct_searchcv_instance（100-110行）
├── _construct_compose_pipeline_instance（113-125行）
├── _construct_sparse_coder（128-140行）
├── test_fit_docstring_attributes（130-200行）
└── _get_all_fitted_attributes（250-275行）
sklearn/tests/test_docstring_parameters_consistency.py
├── CLASS_DOCSTRING_CONSISTENCY_CASES（10-30行）
├── FUNCTION_DOCSTRING_CONSISTENCY_CASES（32-80行）
├── test_class_docstring_consistency（60-70行）
└── test_function_docstring_consistency（75-85行）
sklearn/tests/test_docstrings.py
├── get_all_methods（15-40行）
├── get_all_functions_names（42-50行）
├── filter_errors（42-70行）
├── repr_errors（80-115行）
├── test_function_docstring（125-135行）
├── test_docstring（140-155行）
└── __main__（158-185行）
sklearn/tests/test_config.py
├── test_config_context（15-45行）
├── test_config_context_exception（48-60行）
├── test_set_config（50-70行）
├── set_assume_finite（73-80行）
├── test_config_threadsafe_joblib（50-70行）
├── test_config_threadsafe（75-95行）
└── test_config_array_api_dispatch_error_scipy（98-115行）
sklearn/tests/test_build.py
└── test_openmp_parallelism_enabled（10-30行）
sklearn/tests/test_init.py
└── test_import_skl（15-20行）
sklearn/tests/test_check_build.py
└── test_raise_build_error（10-15行）
sklearn/tests/test_min_dependencies_readme.py
├── TOY_* 玩具数据（15-80行）
├── test_min_dependencies_readme（55-75行）
├── extract_packages_and_pyproject_tags（80-100行）
├── check_pyproject_sections（85-130行）
├── test_min_dependencies_pyproject_toml（135-150行）
├── test_check_matching_pyproject_section（160-175行）
└── test_check_non_matching_pyproject_section（180-210行）
sklearn/tests/__init__.py
└── 空文件（包标记）

以上地图列出本章源码模块及其职责，后文将按数据流逐一解析。

8.15 设计取舍（一问一答）

为什么采用当前方案而不是更复杂的替代方案？ 本章实现优先保证与既有 API 的一致性、可维护性与运行效率。这意味着在少数极端场景下，调用者需要自行在灵活性、内存与速度之间做取舍，换取默认路径的清晰与稳定。

第 9 章 —— 概率校准 —— 让模型输出“可信的置信度”

9.1 学习目标

• 难度：★★★☆☆（3/5）

• 预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础

• 理解概率校准的必要性：分类器原始输出不一定是真实概率，需要后处理修正

• 掌握 CalibratedClassifierCV 的交叉验证管线设计：ensemble 模式与 non-ensemble 模式的区别与适用场景

• 深入三种校准方法（Sigmoid、Isotonic、Temperature）的数学原理与源码实现

• 理解校准器工厂 _fit_calibrator 如何根据方法类型分拣不同的校准策略

• 掌握 Array API 兼容策略：为什么只有 Temperature Scaling 支持跨后端计算

• 理解校准曲线（reliability diagram）的计算原理：分箱统计与可靠性评估

• 能够阅读并分析校准质量测试防线，理解 Brier 分数、标签重映射不变性等验证方法

9.2 生活类比

想象概率校准是一位翻译官的「语义校正」工作：基分类器如同一位说「方言」的演讲者，他的分数（decision_function）或概率（predict_proba）可能带有系统性偏差；校准器则是翻译官，将方言转译为标准的「概率语言」（真实置信度）。其中，Sigmoid 校准是线性翻译器，通过调整斜率（温度）和截距（偏移）来对齐语义；Isotonic 校准是非线性翻译器，不预设线性关系，完全根据数据「逐段」学习映射；Temperature Scaling 则只用一个「温度旋钮」统一控制所有类别的翻译语气，特别适合多分类场景。再想象 CalibratedClassifierCV 的交叉验证管线是一家「质检-校准」双工位流水线：ensemble 模式下，多条独立流水线并行工作，每条流水线用不同的数据子集训练「质检员」（基分类器）和「校准员」（校准器），最终取各流水线输出的平均；non-ensemble 模式下，只用一条流水线，但通过「交叉验证预测」模拟出每个样本的「无偏差质检报告」，然后用这份报告训练唯一的校准员；而 FrozenEstimator 则是已经培训好的质检员直接进入校准工位，所有数据都用来培训校准员。最后，想象校准曲线（reliability diagram）是一张「天气预报准确性报告卡」：分箱将预测概率划分为若干「置信度区间」，如 0-20%、20%-40% 等；prob_pred（x 轴）是每个区间内预报员的平均预测概率，prob_true（y 轴）是每个区间内实际下雨的频率；完美校准落在对角线 y=x 上，预报说 70% 下雨，实际 70% 的时间真的下了雨；校准不足/过度则表现为曲线偏离对角线，需要引入校准器来修正。

9.3 源码地图

sklearn/calibration.py

├── class CalibratedClassifierCV (line 107-560)

│ ├── init() # 初始化参数：estimator、method、cv、n_jobs、ensemble

│ ├── _get_estimator() # 默认估计器解析（LinearSVC）

│ ├── fit() # 拟合校准模型，支持 ensemble/非 ensemble 两种模式

│ ├── predict_proba() # 校准概率预测：对各校准分类器取算术平均

│ ├── predict() # 类别预测：argmax 概率最大类别

│ ├── get_metadata_routing() # 元数据路由配置

│ └── sklearn_tags() # Array API 支持标记

├── _fit_classifier_calibrator_pair() (line 562-620) # ensemble 模式核心：训练单个校准对

├── _fit_calibrator() (line 623-667) # 校准器工厂：按方法分拣校准工具箱

├── class _CalibratedClassifier (line 669-780)

│ ├── init() # 管线封装：合并基分类器与校准器

│ └── predict_proba() # 校准概率计算：OvR 归一化或 Temperature softmax

├── _sigmoid_calibration() (line 783-867) # Platt Scaling 核心：L-BFGS-B 优化

├── _convert_to_logits() (line 870-918) # 统一 logits 转换：概率还原或 (−x, x) 构造

├── class _SigmoidCalibration (line 921-953) # Sigmoid 回归模型

│ ├── fit() # 拟合斜率 a_ 和截距 b_

│ └── predict() # expit(−(a*T + b)) 预测

├── class _TemperatureScaling (line 956-1060) # 温度缩放模型

│ ├── fit() # minimize_scalar 一维优化 log_beta

│ ├── predict() # softmax(beta * logits) 预测

│ └── sklearn_tags() # 标记输入形状约束

├── calibration_curve() (line 1063-1165) # 校准曲线计算：分箱统计与可靠性图谱

└── class CalibrationDisplay (line 1168-1383) # 可视化呈现：可靠性图谱绘图面板

├── init() # 存储曲线数据与元信息

├── plot() # 绘制校准曲线与参考对角线

├── from_estimator() # 从估计器构造（响应值提取）

└── from_predictions() # 从原始预测构造（校准曲线计算）

sklearn/tests/test_calibration.py

├── test_calibration_method_raises() # 无效校准方法参数验证

├── test_calibration() # 校准核心验证：Brier 分数改进、标签重映射

├── test_calibration_default_estimator() # 默认估计器验证（LinearSVC）

├── test_calibration_cv_splitter() # CV 分割器验证

├── test_calibration_cv_nfold() # 折数超限与 LeaveOneOut 禁止验证

├── test_sample_weight() # 样本权重对校准概率的影响

├── test_parallel_execution() # 并行训练与串行训练一致性

├── test_calibration_multiclass() # 多分类校准：概率和为1、Brier 分数改进

├── test_calibration_zero_probability() # 零概率回退：均匀分布

├── test_calibration_frozen() # FrozenEstimator 校准：全量数据用于校准

├── test_calibration_ensemble_false() # non-ensemble 模式手动验证

├── test_sigmoid_calibration() # Sigmoid 校准数值验证

├── test_temperature_scaling() # 温度缩放：log loss 改进、ROC AUC 不变性

├── test_temperature_scaling_input_validation() # 温度缩放输入形状验证

├── test_calibration_curve() # 校准曲线计算：分箱统计验证

├── test_calibration_nan_imputer() # 含 NaN 输入校准

├── test_calibration_prob_sum() # 概率和归一化验证

├── test_calibration_less_classes() # 训练集缺失类别时的校准行为

├── test_calibration_accepts_ndarray() # n 维数组输入验证

├── test_calibration_dict_pipeline() # 字典数据流水线校准

├── test_calibration_attributes() # 属性传播验证

├── test_calibration_inconsistent_prefit_n_features_in() # 预拟合特征数不一致验证

├── test_calibration_votingclassifier() # VotingClassifier 支持验证

├── test_calibration_display_compute() # CalibrationDisplay 与 calibration_curve 一致性

├── test_plot_calibration_curve_pipeline() # 流水线支持验证

├── test_calibration_display_default_labels() # 显示默认标签行为

├── test_calibration_display_label_class_plot() # plot 方法的标签覆盖

├── test_calibration_display_name_multiple_calls() # 多次调用 plot 的 name 覆盖

├── test_calibration_display_ref_line() # 参考线只绘制一次

├── test_calibration_curve_pos_label_error_str() # 字符串标签缺少 pos_label 时的错误

├── test_calibration_curve_pos_label() # 显式 pos_label 的校准曲线行为

├── test_calibration_display_kwargs() # matplotlib 参数别名处理

├── test_calibration_display_pos_label() # pos_label 在显示中的行为

├── test_calibrated_classifier_cv_double_sample_weights_equivalence() # 样本权重与数据重复的等价性

├── test_calibration_with_fit_params() # fit_params 传递验证

├── test_calibration_with_sample_weight_estimator() # sample_weight 传递验证

├── test_calibration_without_sample_weight_estimator() # 估计器不支持 sample_weight 时的警告

├── test_calibration_with_non_sample_aligned_fit_param() # 非样本对齐 fit 参数验证

├── test_calibrated_classifier_cv_works_with_large_confidence_scores() # 大置信度分数下 Sigmoid 与 Isotonic 一致性

├── test_sigmoid_calibration_max_abs_prediction_threshold() # 缩放阈值验证

├── test_float32_predict_proba() # float32 预测概率 dtype 一致性

├── test_error_less_class_samples_than_folds() # 字符串标签 + 折数大于样本数

├── test_temperature_scaling_array_api_compliance() # Array API 兼容性：跨后端温度缩放

└── test_temperature_scaling_array_api_with_str_y_estimator_not_prefit() # 字符串标签 + Array API 温度缩放

9.4 概率校准框架总览 —— 认识这台“概率置信度修正仪”

9.4.1 为什么需要概率校准？

许多分类器输出的 predict_proba 并非真实的概率估计（如 SVM、朴素贝叶斯），需要进行后处理修正。校准的核心思想是学习一个从原始分数到校准概率的单调映射。常用的校准方法包括 Sigmoid（Platt Scaling）、Isotonic（保序回归）和 Temperature Scaling。

9.4.2 CalibratedClassifierCV 的交叉验证管线设计

ensemble=True 时，每折训练一个基分类器 + 校准器对，预测时取各折概率的平均；ensemble=False 时，用 cross_val_predict 获取无偏预测来拟合校准器，最终使用全量数据训练的基分类器。已拟合分类器可通过 FrozenEstimator 包装后直接校准，此时所有数据都用于校准。

9.4.3 响应方法的选择策略

优先使用 decision_function（如果存在），否则回退到 predict_proba。二分类场景下，decision_function 输出一维数组，需重塑为 (n_samples, 1) 以匹配校准器输入格式。

9.4.4 核心类型定义：CalibratedClassifierCV

源码路径：sklearn/calibration.py - CalibratedClassifierCV（第107-560行）

class CalibratedClassifierCV(ClassifierMixin, MetaEstimatorMixin, BaseEstimator):
    """Calibrate probabilities using isotonic, sigmoid, or temperature scaling."""

    _parameter_constraints: dict = {
        "estimator": [
            HasMethods(["fit", "predict_proba"]),
            HasMethods(["fit", "decision_function"]),
            None,
        ],
        "method": [StrOptions({"isotonic", "sigmoid", "temperature"})],
        "cv": ["cv_object"],
        "n_jobs": [Integral, None],
        "ensemble": ["boolean", StrOptions({"auto"})],
    }

    def __init__(
        self,
        estimator=None,
        *,
        method="sigmoid",
        cv=None,
        n_jobs=None,
        ensemble="auto",
    ):
        self.estimator = estimator
        self.method = method
        self.cv = cv
        self.n_jobs = n_jobs
        self.ensemble = ensemble

这段代码定义了校准主类的参数约束与初始化逻辑。estimator 接受具有 fit/predict_proba 或 fit/decision_function 的分类器，默认使用 LinearSVC；method 选择三种校准策略之一；cv 控制交叉验证分割；ensemble 决定是否使用集成模式，"auto" 会在检测到 FrozenEstimator 时自动关闭集成。

9.4.5 逐行解析关键函数：_get_estimator()

源码路径：sklearn/calibration.py - CalibratedClassifierCV._get_estimator()（第335-345行）

    def _get_estimator(self):
        """Resolve which estimator to return (default is LinearSVC)"""
        if self.estimator is None:
            # we want all classifiers that don't expose a random_state
            # to be deterministic (and we don't want to expose this one).
            estimator = LinearSVC(random_state=0)
            if _routing_enabled():
                estimator.set_fit_request(sample_weight=True)
        else:
            estimator = self.estimator

        return estimator

这段代码实现了默认估计器的解析逻辑：当用户未提供 estimator 时，返回一个固定随机种子的 LinearSVC 并开启 sample_weight 路由支持，保证可复现性与元数据路由兼容。

9.4.6 逐行解析关键函数：fit() —— ensemble 分支

源码路径：sklearn/calibration.py - CalibratedClassifierCV.fit()（第347-490行，重点 467-480 行）

        if _ensemble:
            parallel = Parallel(n_jobs=self.n_jobs)
            self.calibrated_classifiers_ = parallel(
                delayed(_fit_classifier_calibrator_pair)(
                    clone(estimator),
                    X,
                    y,
                    train=train,
                    test=test,
                    method=self.method,
                    classes=self.classes_,
                    xp=xp,
                    sample_weight=sample_weight,
                    fit_params=routed_params.estimator.fit,
                )
                for train, test in cv.split(X, y, **routed_params.splitter.split)
            )

这段代码展示了 ensemble 模式的并行训练核心：使用 joblib.Parallel 并行执行各折的 _fit_classifier_calibrator_pair，每折克隆独立的估计器实例，按 train/test 索引切分数据与 fit_params，最后收集所有校准对到 calibrated_classifiers_ 列表。

9.4.7 逐行解析关键函数：fit() —— non-ensemble 分支

源码路径：sklearn/calibration.py - CalibratedClassifierCV.fit()（第482-510行）

        else:
            this_estimator = clone(estimator)
            method_name = _check_response_method(
                this_estimator,
                ["decision_function", "predict_proba"],
            ).__name__
            predictions = cross_val_predict(
                estimator=this_estimator,
                X=X,
                y=y,
                cv=cv,
                method=method_name,
                n_jobs=self.n_jobs,
                params=routed_params.estimator.fit,
            )
            if self.classes_.shape[0] == 2:
                # Ensure shape (n_samples, 1) in the binary case
                if method_name == "predict_proba":
                    # Select the probability column of the positive class
                    predictions = _process_predict_proba(
                        y_pred=predictions,
                        target_type="binary",
                        classes=self.classes_,
                        pos_label=self.classes_[1],
                    )
                predictions = predictions.reshape(-1, 1)

            if sample_weight is not None:
                # Check that the sample_weight dtype is consistent with the
                # predictions to avoid unintentional upcasts.
                sample_weight = _check_sample_weight(
                    sample_weight, predictions, dtype=predictions.dtype
                )

            this_estimator.fit(X, y, **routed_params.estimator.fit)
            # Note: Here we don't pass on fit_params because the supported
            # calibrators don't support fit_params anyway
            calibrated_classifier = _fit_calibrator(
                this_estimator,
                predictions,
                y,
                self.classes_,
                self.method,
                xp=xp,
                sample_weight=sample_weight,
            )
            self.calibrated_classifiers_.append(calibrated_classifier)

这段代码实现了 non-ensemble 模式：先用 cross_val_predict 获得每个样本在未见过折上的无偏预测，二分类时提取正类概率并重塑为列向量；然后用全量数据重新拟合基础分类器，最后用交叉验证预测与全量标签拟合唯一校准器，calibrated_classifiers_ 长度固定为 1。

9.4.8 逐行解析关键函数：predict_proba() 与 predict()

源码路径：sklearn/calibration.py - CalibratedClassifierCV.predict_proba()（第512-528行）与 predict()（第530-545行）

    def predict_proba(self, X):
        check_is_fitted(self)
        # Compute the arithmetic mean of the predictions of the calibrated
        # classifiers
        xp, _, device_ = get_namespace_and_device(X)
        mean_proba = xp.zeros((_num_samples(X), self.classes_.shape[0]), device=device_)
        for calibrated_classifier in self.calibrated_classifiers_:
            proba = calibrated_classifier.predict_proba(X)
            mean_proba += proba

        mean_proba /= len(self.calibrated_classifiers_)

        return mean_proba

    def predict(self, X):
        xp, _ = get_namespace(X)
        check_is_fitted(self)
        class_indices = xp.argmax(self.predict_proba(X), axis=1)
        if isinstance(self.classes_[0], str):
            class_indices = _convert_to_numpy(class_indices, xp=xp)

        return self.classes_[class_indices]

这段代码实现了预测阶段的集成逻辑：predict_proba 对所有校准分类器的概率输出取算术平均；predict 则在平均概率上取 argmax 得到类别索引，字符串标签时需转回 NumPy 再索引 classes_。

9.4.9 完整流程图：CalibratedClassifierCV.fit() 双模式决策

flowchart TD A[fit(X, y, sample_weight, **fit_params)] --> B{ensemble 模式?} B -->|True| C[并行: 对每个 CV 折] C --> C1[clone(estimator)] C1 --> C2[train 索引拟合基分类器] C2 --> C3[test 索引获取原始预测] C3 --> C4[_fit_calibrator 拟合校准器] C4 --> C5[收集到 calibrated_classifiers_] B -->|False| D[串行: cross_val_predict] D --> D1[获得无偏预测 predictions] D1 --> D2[二分类提取正类概率 reshape] D2 --> D3[全量数据重新拟合基分类器] D3 --> D4[_fit_calibrator 拟合唯一校准器] D4 --> D5[calibrated_classifiers_ 长度=1] C5 & D5 --> E[提取 n_features_in_, feature_names_in_] E --> F[返回 self]

9.5 ensemble 模式 —— 训练多条“校准流水线”并取平均

9.5.1 为什么需要 ensemble 模式？

每个交叉验证折都独立训练基分类器和校准器，降低单折预测偏差对校准结果的影响。预测时对多个校准分类器的概率输出取算术平均，提升概率估计的稳定性。

9.5.2 _fit_classifier_calibrator_pair 的工作流程

从 X, y 中按 train/test 索引切分出训练子集和校准子集；先用训练子集拟合克隆的基础估计器，再在测试子集上获取原始预测；将原始预测和真实标签交给 _fit_calibrator 拟合校准器。

9.5.3 并行执行的细节

使用 Parallel 和 delayed 并行训练各折的 estimator-calibrator 对；fit_params 通过 _check_method_params 按 train 索引筛选，确保每个子估计器只接收对应折的训练参数；sample_weight 也按 test 索引筛选取出校准集的样本权重。

9.5.4 逐行解析关键函数：_fit_classifier_calibrator_pair()

源码路径：sklearn/calibration.py - _fit_classifier_calibrator_pair()（第562-620行）

def _fit_classifier_calibrator_pair(
    estimator,
    X,
    y,
    train,
    test,
    method,
    classes,
    xp,
    sample_weight=None,
    fit_params=None,
):
    """Fit a classifier/calibration pair on a given train/test split."""
    fit_params_train = _check_method_params(X, params=fit_params, indices=train)
    X_train, y_train = _safe_indexing(X, train), _safe_indexing(y, train)
    X_test, y_test = _safe_indexing(X, test), _safe_indexing(y, test)

    estimator.fit(X_train, y_train, **fit_params_train)

    predictions, _ = _get_response_values(
        estimator,
        X_test,
        response_method=["decision_function", "predict_proba"],
    )
    if predictions.ndim == 1:
        # Reshape binary output from `(n_samples,)` to `(n_samples, 1)`
        predictions = predictions.reshape(-1, 1)

    if sample_weight is not None:
        sample_weight = _check_sample_weight(sample_weight, X, dtype=predictions.dtype)
        sw_test = _safe_indexing(sample_weight, test)
    else:
        sw_test = None
    calibrated_classifier = _fit_calibrator(
        estimator,
        predictions,
        y_test,
        classes,
        method,
        xp=xp,
        sample_weight=sw_test,
    )
    return calibrated_classifier

这段代码完整展示了单折训练流程：按索引切分数据与 fit_params，拟合克隆估计器，获取测试集原始预测（优先 decision_function），二分类时重塑为列向量，切分校准集样本权重，最后调用 _fit_calibrator 生成校准管线对象。

9.6 非 ensemble 模式 —— 借用 cross_val_predict 的“无偏预测”

9.6.1 为什么需要 non-ensemble 模式？

适用于基分类器训练成本高的场景（如大规模数据），只需训练一个最终分类器。通过 cross_val_predict 获得每个样本在“未见过”折上的预测，避免用训练集自身预测导致的过度乐观。

9.6.2 二分类场景的特殊处理

当基础分类器提供 predict_proba 时，只提取正类概率列（pos_label=self.classes_[1]）；将提取的概率重塑为 (n_samples, 1)，与 decision_function 的输出格式保持一致。

9.6.3 最终模型的组装

先用全量数据重新拟合基础分类器；再用交叉验证预测和全量标签拟合唯一校准器；最终 calibrated_classifiers_ 列表长度为 1。

9.6.4 逐行解析关键函数：fit() non-ensemble 分支（复习）

源码路径：sklearn/calibration.py - CalibratedClassifierCV.fit()（第482-510行）

上一节已详细解析，此处不再赘述。核心区别在于：ensemble 训练多条流水线取平均，non-ensemble 只训练一条流水线但用交叉验证预测作为校准输入。

9.7 校准器工厂 —— 通过 _fit_calibrator 分拣“校准工具箱”

9.7.1 Sigmoid 和 Isotonic 的 OvR 策略

使用 label_binarize 将多分类标签转为 one-hot 编码矩阵 Y；为每个正类列独立拟合一个校准器，形成 OvR（One-vs-Rest）结构；二分类时只需一个校准器，对应正类的原始预测。

9.7.2 Temperature Scaling 的特殊路径

多分类场景天然支持：直接对 logits 或概率做温度缩放；二分类且 predict_proba 输出单列时，需拼接 1 - predictions 构造两列概率；温度校准器只需一个，对所有类共享同一个温度参数。

9.7.3 _CalibratedClassifier 管线封装

将基础分类器和校准器列表组合为一个类管线对象；保存 classes 和方法信息，供 predict_proba 阶段使用。

9.7.4 逐行解析关键函数：_fit_calibrator()

源码路径：sklearn/calibration.py - _fit_calibrator()（第623-667行）

def _fit_calibrator(clf, predictions, y, classes, method, xp, sample_weight=None):
    """Fit calibrator(s) and return a `_CalibratedClassifier` instance."""
    calibrators = []

    if method in ("isotonic", "sigmoid"):
        Y = label_binarize(y, classes=classes)
        label_encoder = LabelEncoder().fit(classes)
        pos_class_indices = label_encoder.transform(clf.classes_)
        for class_idx, this_pred in zip(pos_class_indices, predictions.T):
            if method == "isotonic":
                calibrator = IsotonicRegression(out_of_bounds="clip")
            else:  # "sigmoid"
                calibrator = _SigmoidCalibration()
            calibrator.fit(this_pred, Y[:, class_idx], sample_weight)
            calibrators.append(calibrator)
    elif method == "temperature":
        if classes.shape[0] == 2 and predictions.shape[-1] == 1:
            response_method_name = _check_response_method(
                clf,
                ["decision_function", "predict_proba"],
            ).__name__
            if response_method_name == "predict_proba":
                predictions = xp.concat([1 - predictions, predictions], axis=1)
        calibrator = _TemperatureScaling()
        calibrator.fit(predictions, y, sample_weight)
        calibrators.append(calibrator)

    pipeline = _CalibratedClassifier(clf, calibrators, method=method, classes=classes)
    return pipeline

这段代码是校准器工厂的核心：Sigmoid/Isotonic 走 OvR 循环，逐类拟合校准器；Temperature 走单校准器路径，二分类 predict_proba 时先拼接两列。最终将分类器与校准器列表封装进 _CalibratedClassifier。

9.7.5 核心类型定义：_CalibratedClassifier

源码路径：sklearn/calibration.py - _CalibratedClassifier（第669-780行）

class _CalibratedClassifier:
    """Pipeline-like chaining a fitted classifier and its fitted calibrators."""

    def __init__(self, estimator, calibrators, *, classes, method="sigmoid"):
        self.estimator = estimator
        self.calibrators = calibrators
        self.classes = classes
        self.method = method

这段代码定义了校准管线的轻量级容器，仅保存基分类器、校准器列表、类别标签与校准方法，供后续 predict_proba 调用。

9.8 概率预测与归一化 —— 校准结果的“最终装配线”

9.8.1 Sigmoid/Isotonic 的预测与归一化

二分类时，正类校准概率 = calibrator.predict(predictions)，负类概率 = 1 - 正类概率；多分类时，各 OvR 校准器的输出可能不归一化，需要除以行和；当所有校准器对某样本都返回零概率时，回退到均匀分布（1/n_classes），避免除零错误。

9.8.2 Temperature 的预测路径

二分类且 predict_proba 时，先拼接两列概率；直接调用温度校准器的 predict，内部进行 logits 转换和 softmax 缩放。

9.8.3 数值稳定性保护

将微小的概率上溢（如 1.0 < proba <= 1.0 + 1e-5）截断为 1.0；在 CalibratedClassifierCV.predict_proba 中对多个校准分类器的概率取均值；预测类别时，使用 argmax 选择概率最大的类别索引。

9.8.4 逐行解析关键函数：_CalibratedClassifier.predict_proba()

源码路径：sklearn/calibration.py - _CalibratedClassifier.predict_proba()（第698-780行）

    def predict_proba(self, X):
        predictions, _ = _get_response_values(
            self.estimator,
            X,
            response_method=["decision_function", "predict_proba"],
        )
        if predictions.ndim == 1:
            predictions = predictions.reshape(-1, 1)

        n_classes = self.classes.shape[0]

        proba = np.zeros((_num_samples(X), n_classes))

        if self.method in ("sigmoid", "isotonic"):
            label_encoder = LabelEncoder().fit(self.classes)
            pos_class_indices = label_encoder.transform(self.estimator.classes_)
            for class_idx, this_pred, calibrator in zip(
                pos_class_indices, predictions.T, self.calibrators
            ):
                if n_classes == 2:
                    class_idx += 1
                proba[:, class_idx] = calibrator.predict(this_pred)
            # Normalize the probabilities
            if n_classes == 2:
                proba[:, 0] = 1.0 - proba[:, 1]
            else:
                denominator = np.sum(proba, axis=1)[:, np.newaxis]
                uniform_proba = np.full_like(proba, 1 / n_classes)
                proba = np.divide(
                    proba, denominator, out=uniform_proba, where=denominator != 0
                )
        elif self.method == "temperature":
            xp, _ = get_namespace(predictions)
            if n_classes == 2 and predictions.shape[-1] == 1:
                response_method_name = _check_response_method(
                    self.estimator,
                    ["decision_function", "predict_proba"],
                ).__name__
                if response_method_name == "predict_proba":
                    predictions = xp.concat([1 - predictions, predictions], axis=1)
            proba = self.calibrators[0].predict(predictions)

        # Deal with cases where the predicted probability minimally exceeds 1.0
        proba[(1.0 < proba) & (proba <= 1.0 + 1e-5)] = 1.0

        return proba

这段代码实现了三种校准方法的预测分支：Sigmoid/Isotonic 遍历 OvR 校准器填充概率矩阵，二分类直接补全负类，多分类归一化并处理零概率回退；Temperature 统一转 logits 再 softmax，二分类 predict_proba 时先拼接。最后做数值裁剪保证概率不超过 1.0。

9.8.5 完整数据流图：从原始预测到校准概率

flowchart TD A[原始预测 predictions] --> B{校准方法} B -->|Sigmoid/Isotonic| C[OvR 循环: 每类一个校准器] C --> C1[calibrator.predict(this_pred)] C1 --> C2[填充 proba[:, class_idx]] C2 --> C3{n_classes == 2?} C3 -->|是| C4[proba[:, 0] = 1 - proba[:, 1]] C3 -->|否| C5[行归一化 / 零概率回退均匀分布] B -->|Temperature| D[_convert_to_logits -> logits] D --> D1[softmax(beta * logits)] C4 & C5 & D1 --> E[裁剪微小超出 1.0 的值] E --> F[返回校准概率]

9.9 Sigmoid 校准 —— Platt Scaling 的“数值稳定化改造”

9.9.1 Platt 的原始思想

假设原始预测值 F 与正类概率之间满足 sigmoid 关系：P(y=1|F) = 1/(1+exp(A*F+B))；通过最小化负对数似然来拟合参数 A（斜率）和 B（截距）。

9.9.2 大分数值的缩放预处理

当 max(|F|) >= 30 时，将 F 除以 max(|F|) 进行缩放，避免 sigmoid 函数饱和；缩放不影响最终结果，因为线性模型对特征缩放具有不变性。

9.9.3 贝叶斯先验与目标修正

使用 Platt 提出的先验修正：T_+ = (N_+ + 1)/(N_+ + 2), T_- = 1/(N_- + 2)；先验计数支持 sample_weight 加权求和；目标值 T 代替硬 0/1 标签可以缓解过拟合。

9.9.4 L-BFGS-B 优化

使用 scipy.optimize.minimize 的 L-BFGS-B 方法求解；初始值 AB0 = [0, log((prior0+1)/(prior1+1))]，截距初始化为先验对数几率；梯度通过 HalfBinomialLoss 的 loss_gradient 计算，保证 dtype 一致性。

9.9.5 逐行解析关键函数：_sigmoid_calibration()

源码路径：sklearn/calibration.py - _sigmoid_calibration()（第783-867行）

def _sigmoid_calibration(
    predictions, y, sample_weight=None, max_abs_prediction_threshold=30
):
    """Probability Calibration with sigmoid method (Platt 2000)"""
    predictions = column_or_1d(predictions)
    y = column_or_1d(y)

    F = predictions  # F follows Platt's notations

    scale_constant = 1.0
    max_prediction = np.max(np.abs(F))

    # If the predictions have large values we scale them in order to bring
    # them within a suitable range. This has no effect on the final
    # (prediction) result because linear models like Logisitic Regression
    # without a penalty are invariant to multiplying the features by a
    # constant.
    if max_prediction >= max_abs_prediction_threshold:
        scale_constant = max_prediction
        F = F / scale_constant

    # Bayesian priors (see Platt end of section 2.2):
    mask_negative_samples = y <= 0
    if sample_weight is not None:
        prior0 = (sample_weight[mask_negative_samples]).sum()
        prior1 = (sample_weight[~mask_negative_samples]).sum()
    else:
        prior0 = float(np.sum(mask_negative_samples))
        prior1 = y.shape[0] - prior0
    T = np.zeros_like(y, dtype=predictions.dtype)
    T[y > 0] = (prior1 + 1.0) / (prior1 + 2.0)
    T[y <= 0] = 1.0 / (prior0 + 2.0)

    bin_loss = HalfBinomialLoss()

    def loss_grad(AB):
        raw_prediction = -(AB[0] * F + AB[1]).astype(dtype=predictions.dtype)
        l, g = bin_loss.loss_gradient(
            y_true=T,
            raw_prediction=raw_prediction,
            sample_weight=sample_weight,
        )
        loss = l.sum()
        grad = np.asarray([-g @ F, -g.sum()], dtype=np.float64)
        return loss, grad

    AB0 = np.array([0.0, log((prior0 + 1.0) / (prior1 + 1.0))])

    opt_result = minimize(
        loss_grad,
        AB0,
        method="L-BFGS-B",
        jac=True,
        options={
            "gtol": 1e-6,
            "ftol": 64 * np.finfo(float).eps,
        },
    )
    AB_ = opt_result.x

    return AB_[0] / scale_constant, AB_[1]

这段代码完整实现了 Platt Scaling：先对大数值特征做缩放预处理，计算带先验的软目标 T，构造 HalfBinomialLoss 的损失梯度函数，用 L-BFGS-B 优化斜率和截距，最后将斜率按缩放常数还原。

9.9.6 核心类型定义：_SigmoidCalibration

源码路径：sklearn/calibration.py - _SigmoidCalibration（第921-953行）

class _SigmoidCalibration(RegressorMixin, BaseEstimator):
    """Sigmoid regression model."""

    def fit(self, X, y, sample_weight=None):
        X = column_or_1d(X)
        y = column_or_1d(y)
        X, y = indexable(X, y)

        self.a_, self.b_ = _sigmoid_calibration(X, y, sample_weight)
        return self

    def predict(self, T):
        T = column_or_1d(T)
        return expit(-(self.a_ * T + self.b_))

这段代码将 _sigmoid_calibration 封装为兼容 scikit-learn 接口的回归器：fit 调用核心优化得到 a_、b_；predict 直接套用 sigmoid 公式 expit(-(a*T + b))。

9.10 Temperature Scaling —— 一条温度参数“熨平”全部类别的概率分布

9.10.1 温度缩放的数学原理

对 logits 除以温度 T（或乘以逆温度 beta），使概率分布更陡峭或更平滑；温度 T > 1 时概率分布更平滑，T < 1 时更尖锐；优化目标是最小化负对数似然（即 log loss），参数只有 log_beta 一个标量。

9.10.2 logits 输入的统一转换

decision_function 输出一维时，转换为 (-x, x) 的两列 logits；predict_proba 输出满足概率性质（非负且行和为 1）时，用 log(p + eps) 还原为 logits；使用 eps=1e-12 防止 log(0) 产生无穷大。

9.10.3 一维标量优化

使用 scipy.optimize.minimize_scalar 在 [-10, 10] 区间内搜索最优 log_beta；优化 tolerance 设为 64 * eps，保证 float32/float64 下的精度；优化成功后，beta_ = exp(log_beta_opt)。

9.10.4 逐行解析关键函数：_convert_to_logits()

源码路径：sklearn/calibration.py - _convert_to_logits()（第870-918行）

def _convert_to_logits(decision_values, eps=1e-12, xp=None):
    """Convert decision_function values to 2D and predict_proba values to logits."""
    xp, _, device_ = get_namespace_and_device(decision_values, xp=xp)
    decision_values = check_array(
        decision_values, dtype=[xp.float64, xp.float32], ensure_2d=False
    )
    if (decision_values.ndim == 2) and (decision_values.shape[1] > 1):
        entries_zero_to_one = xp.all((decision_values >= 0) & (decision_values <= 1))
        row_sums_to_one = xp.all(
            xpx.isclose(
                xp.sum(decision_values, axis=1),
                xp.asarray(1.0, device=device_, dtype=decision_values.dtype),
            )
        )

        if entries_zero_to_one and row_sums_to_one:
            logits = xp.log(decision_values + eps)
        else:
            logits = decision_values

    elif (decision_values.ndim == 2) and (decision_values.shape[1] == 1):
        logits = xp.concat([-decision_values, decision_values], axis=1)

    elif decision_values.ndim == 1:
        decision_values = xp.reshape(decision_values, (-1, 1))
        logits = xp.concat([-decision_values, decision_values], axis=1)

    return logits

这段代码统一了 logits 转换逻辑：二维多列输入先判断是否为概率（非负且行和≈1），是则取对数还原 logits，否则视为已是 logits；二维单列或一维输入则构造 (-x, x) 两列 logits。全程使用 Array API 命名空间 xp 保证跨后端兼容。

9.10.5 逐行解析关键函数：_TemperatureScaling.fit()

源码路径：sklearn/calibration.py - _TemperatureScaling.fit()（第956-1020行）

    def fit(self, X, y, sample_weight=None):
        xp, _, xp_device = get_namespace_and_device(X, y)
        X, y = indexable(X, y)
        check_consistent_length(X, y)
        logits = _convert_to_logits(X)  # guarantees xp.float64 or xp.float32

        dtype_ = logits.dtype
        labels = column_or_1d(y, dtype=dtype_)

        if sample_weight is not None:
            sample_weight = _check_sample_weight(sample_weight, labels, dtype=dtype_)

        is_numpy_namespace = _is_numpy_namespace(xp)
        multinomial_loss = (
            HalfMultinomialLoss(n_classes=logits.shape[1])
            if is_numpy_namespace
            else HalfMultinomialLossArrayAPI(n_classes=logits.shape[1])
        )

        def log_loss(log_beta=0.0):
            log_beta = xp.asarray(log_beta, dtype=dtype_, device=xp_device)
            raw_prediction = xp.exp(log_beta) * logits
            return multinomial_loss(
                labels,
                raw_prediction,
                sample_weight,
                xp=xp,
            )

        xatol = 64 * xp.finfo(dtype_).eps
        log_beta_minimizer = minimize_scalar(
            log_loss,
            bounds=(-10.0, 10.0),
            options={
                "xatol": xatol,
            },
        )

        if not log_beta_minimizer.success:
            raise RuntimeError(
                "Temperature scaling fails to optimize during calibration. "
                f"Reason from `scipy.optimize.minimize_scalar`: "
                f"{log_beta_minimizer.message}"
            )

        self.beta_ = xp.exp(
            xp.asarray(log_beta_minimizer.x, dtype=dtype_, device=xp_device)
        )

        return self

这段代码实现了温度缩放的拟合：统一转 logits，根据后端选择 HalfMultinomialLoss 或 HalfMultinomialLossArrayAPI，构造以 log_beta 为变量的 log loss 函数，用 minimize_scalar 在 [-10, 10] 搜索最优值，最后指数变换得到 beta_。

9.10.6 逐行解析关键函数：_TemperatureScaling.predict()

源码路径：sklearn/calibration.py - _TemperatureScaling.predict()（第1022-1045行）

    def predict(self, X):
        logits = _convert_to_logits(X)
        return softmax(self.beta_ * logits)

这段代码极其简洁：转 logits，乘以学到的 beta_，再做 softmax 得到校准概率。

9.10.7 核心类型定义：_TemperatureScaling.sklearn_tags()

源码路径：sklearn/calibration.py - _TemperatureScaling.__sklearn_tags__()（第1047-1052行）

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.input_tags.one_d_array = True
        tags.input_tags.two_d_array = False
        return tags

这段代码声明该校准器接受一维数组输入（二分类 decision_function），不接受二维数组，供标签系统校验使用。

9.11 Array API 兼容策略 —— 让温度校准“跨后端旅行”

9.11.1 为什么只有 temperature 支持 Array API？

Sigmoid 校准依赖 scipy.optimize.minimize，Isotonic 校准依赖 sklearn.isotonic，这两者仅支持 NumPy；Temperature Scaling 只需一维标量优化和 softmax 计算，可以完全用 Array API 命名空间实现。

9.11.2 损失函数的多态选择

NumPy 后端使用 HalfMultinomialLoss（Cython 实现的高性能版本）；其他 Array API 后端使用 HalfMultinomialLossArrayAPI（纯 Python + xp 命名空间实现）。

9.11.3 dtype 与 device 的透传

log_beta 被转换为与 logits 相同的 dtype 和 device；raw_prediction = exp(log_beta) * logits 确保计算不发生隐式类型提升；__sklearn_tags__ 标记 array_api_support = estimator_tags.array_api_support and method == "temperature"。

9.11.4 逐行解析关键函数：CalibratedClassifierCV.sklearn_tags()

源码路径：sklearn/calibration.py - CalibratedClassifierCV.__sklearn_tags__()（第552-560行）

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        estimator_tags = get_tags(self._get_estimator())
        tags.input_tags.sparse = estimator_tags.input_tags.sparse
        tags.array_api_support = (
            estimator_tags.array_api_support and self.method == "temperature"
        )
        return tags

这段代码精准控制了 Array API 支持标记：仅当基础估计器支持 Array API 且校准方法为 temperature 时，整个校准器才宣称支持 Array API。

9.11.5 逐行解析关键函数：CalibratedClassifierCV.get_metadata_routing()

源码路径：sklearn/calibration.py - CalibratedClassifierCV.get_metadata_routing()（第548-565行）

    def get_metadata_routing(self):
        router = (
            MetadataRouter(owner=self)
            .add_self_request(self)
            .add(
                estimator=self._get_estimator(),
                method_mapping=MethodMapping().add(caller="fit", callee="fit"),
            )
            .add(
                splitter=self.cv,
                method_mapping=MethodMapping().add(caller="fit", callee="split"),
            )
        )
        return router

这段代码声明了元数据路由拓扑：self 接收自身请求，estimator 的 fit 方法接收调用者的 fit 请求，splitter（即 cv）的 split 方法接收调用者的 fit 请求。

9.12 校准曲线计算 —— 用“分箱统计”画出可靠性图谱

9.12.1 校准曲线的核心思想

将预测概率 [0,1] 区间划分为若干 bins；每个 bin 内计算平均预测概率（x 轴）和实际正类比例（y 轴）；完美校准时曲线应为对角线 y=x。

9.12.2 两种分箱策略

uniform：等宽分箱，使用 np.linspace(0, 1, n_bins+1) 确定边界；quantile：等频分箱，使用 np.percentile 确定边界，保证每个 bin 样本数接近。

9.12.3 np.searchsorted 与 np.bincount 的高效分箱

binids = np.searchsorted(bins[1:-1], y_prob) 将样本分配到对应 bin；用 np.bincount 分别累计每个 bin 的样本数、正类数和预测概率和；剔除空 bin（bin_total != 0），避免除零。

9.12.4 输入验证

检查 y_prob 是否在 [0,1] 范围内；只支持二分类（len(np.unique(y_true)) <= 2）；pos_label 用于将多格式标签（如 {-1,1} 或字符串）统一转为布尔正类掩码。

9.12.5 逐行解析关键函数：calibration_curve()

源码路径：sklearn/calibration.py - calibration_curve()（第1063-1165行）

def calibration_curve(
    y_true,
    y_prob,
    *,
    pos_label=None,
    n_bins=5,
    strategy="uniform",
):
    y_true = column_or_1d(y_true)
    y_prob = column_or_1d(y_prob)
    check_consistent_length(y_true, y_prob)
    pos_label = _check_pos_label_consistency(pos_label, y_true)

    if y_prob.min() < 0 or y_prob.max() > 1:
        raise ValueError("y_prob has values outside [0, 1].")

    labels = np.unique(y_true)
    if len(labels) > 2:
        raise ValueError(
            f"Only binary classification is supported. Provided labels {labels}."
        )
    y_true = y_true == pos_label

    if strategy == "quantile":
        quantiles = np.linspace(0, 1, n_bins + 1)
        bins = np.percentile(y_prob, quantiles * 100)
    elif strategy == "uniform":
        bins = np.linspace(0.0, 1.0, n_bins + 1)
    else:
        raise ValueError(
            "Invalid entry to 'strategy' input. Strategy "
            "must be either 'quantile' or 'uniform'."
        )

    binids = np.searchsorted(bins[1:-1], y_prob)

    bin_sums = np.bincount(binids, weights=y_prob, minlength=len(bins))
    bin_true = np.bincount(binids, weights=y_true, minlength=len(bins))
    bin_total = np.bincount(binids, minlength=len(bins))

    nonzero = bin_total != 0
    prob_true = bin_true[nonzero] / bin_total[nonzero]
    prob_pred = bin_sums[nonzero] / bin_total[nonzero]

    return prob_true, prob_pred

这段代码实现了校准曲线的核心计算：输入检查与标签二值化，按策略生成分箱边界，searchsorted 分箱，bincount 加权累计求均值，剔除空 bin 返回有效坐标点。

9.13 CalibrationDisplay 可视化 —— 组装“可靠性图谱的绘图面板”

9.13.1 类继承与混入

继承 _BinaryClassifierCurveDisplayMixin，获得 _validate_plot_params、_validate_and_get_response_values 等通用方法；混入提供了从估计器或原始预测两路构造的验证逻辑。

9.13.2 from_estimator 构造器

使用 _validate_and_get_response_values 从估计器获取 predict_proba 输出；自动确定正类标签（pos_label）和曲线名称（name）；委托给 from_predictions 完成后续计算和绘图。

9.13.3 plot 方法的关键细节

默认曲线样式：marker="s", linestyle="-"；绘制“Perfectly calibrated”参考对角线，且只绘制一次（避免多曲线叠加时重复图例）；坐标轴标签包含正类信息，如 Mean predicted probability (Positive class: 1)；图例始终显示在右下角（loc="lower right"）。

9.13.4 逐行解析关键函数：CalibrationDisplay.plot()

源码路径：sklearn/calibration.py - CalibrationDisplay.plot()（第1178-1217行）

    def plot(self, *, ax=None, name=None, ref_line=True, **kwargs):
        self.ax_, self.figure_, name = self._validate_plot_params(ax=ax, name=name)

        info_pos_label = (
            f"(Positive class: {self.pos_label})" if self.pos_label is not None else ""
        )

        default_line_kwargs = {"marker": "s", "linestyle": "-"}
        if name is not None:
            default_line_kwargs["label"] = name
        line_kwargs = _validate_style_kwargs(default_line_kwargs, kwargs)

        ref_line_label = "Perfectly calibrated"
        existing_ref_line = ref_line_label in self.ax_.get_legend_handles_labels()[1]
        if ref_line and not existing_ref_line:
            self.ax_.plot([0, 1], [0, 1], "k:", label=ref_line_label)
        self.line_ = self.ax_.plot(self.prob_pred, self.prob_true, **line_kwargs)[0]

        self.ax_.legend(loc="lower right")

        xlabel = f"Mean predicted probability {info_pos_label}"
        ylabel = f"Fraction of positives {info_pos_label}"
        self.ax_.set(xlabel=xlabel, ylabel=ylabel)

        return self

这段代码实现了绘图逻辑：验证坐标轴与曲线名，构建带正类信息的轴标签，绘制校准曲线与参考对角线（去重），设置图例位置，返回自身支持链式调用。

9.14 设计中的取舍

9.14.1 为什么 Sigmoid 和 Isotonic 采用 OvR 而 Temperature 天然支持多分类？

Sigmoid/Isotonic 本质上是二分类校准器，将多分类拆解为多个二分类问题（OvR）是自然的扩展方式，但会导致各类校准概率和不为 1，需要事后归一化；Temperature Scaling 直接在 logits 空间操作，通过单一温度参数缩放所有类别的 logits 再 softmax，天然保持概率和为 1，且参数量极少（仅 1 个标量），不易过拟合。

9.14.2 ensemble 与 non-ensemble 的 trade-off 是什么？

ensemble 模式训练 n_folds 个基分类器，计算成本高但预测更稳健（集成平均），适合基分类器训练快、数据量适中的场景；non-ensemble 只训练 1 个最终分类器，计算成本低，但校准器依赖 cross_val_predict 的无偏预测，适合基分类器训练慢、数据量大的场景。"auto" 默认开启 ensemble，仅在 FrozenEstimator 时关闭，体现了“默认稳健、显式优化”的设计哲学。

9.14.3 为什么 Temperature Scaling 的 beta_ 优化区间限制在 [-10, 10]？

对数温度 log_beta 在 [-10, 10] 对应温度 T = 1/beta 约在 [4.5e-5, 2.2e4]，足以覆盖从极度平滑到极度尖锐的概率分布；更宽的区间会导致数值溢出（exp(10) ≈ 22026 乘以 logits 可能溢出 float32）或优化困难，当前区间是数值稳定性与表达能力的平衡。

9.15 动手练习

• 练习 1：阅读 CalibratedClassifierCV.fit() 的 ensemble 分支

  阅读 sklearn/calibration.py 第440-480行，理解 ensemble 模式下 Parallel 如何并行训练校准对：

  1. delayed(_fit_classifier_calibrator_pair) 的调用方式

  2. train 和 test 索引如何从 cv.split() 获得

  3. fit_params 如何通过 routed_params.estimator.fit 传递

  回答问题：

  • 为什么每个折需要 clone(estimator) 而不是直接使用同一个实例？

  • sample_weight 在 ensemble 模式下如何被切分为校准子集的权重？

• 练习 2：分析 _fit_calibrator 的 OvR 策略

  阅读 sklearn/calibration.py 第623-667行，理解 _fit_calibrator 如何为不同校准方法分拣策略：

  1. Sigmoid/Isotonic 使用 label_binarize 将多分类标签转为 one-hot 矩阵

  2. 二分类时只拟合一个校准器（对应正类），预测时负类概率 = 1 - 正类概率

  3. Temperature 在二分类 + predict_proba 单列输出时如何拼接两列概率

  回答问题：

  • 多分类 Sigmoid/Isotonic 的 OvR 校准器输出为何需要归一化？

  • Temperature Scaling 为什么只用一个校准器？

• 练习 3：阅读 _convert_to_logits 的输入检测逻辑

  阅读 sklearn/calibration.py 第870-918行，理解 _convert_to_logits 如何区分 decision_function 和 predict_proba 输出：

  1. 如何判断输入是概率输出？条件是什么？

  2. 一维 decision_function 输出如何转换为 (-x, x) 两列？

  3. eps=1e-12 在 log 转换中的作用是什么？

  回答问题：

  • 如果 predict_proba 输出的概率不满足行和为1，_convert_to_logits 会如何处理？

  • 为什么概率还原为 logits 时需要加 eps？

• 练习 4：分析 calibration_curve 的分箱统计实现

  阅读 sklearn/calibration.py 第1063-1165行，理解 calibration_curve 的分箱计算流程：

  1. uniform 策略下 np.linspace(0, 1, n_bins+1) 如何确定边界

  2. quantile 策略下 np.percentile 如何确定边界

  3. np.searchsorted 和 np.bincount 如何高效完成样本分箱

  回答问题：

  • 为什么 bin_sums 和 bin_total 使用 minlength=len(bins)？

  • 如何剔除空 bin？为什么要剔除？

• 练习 5：分析 _CalibratedClassifier.predict_proba 的归一化与回退策略

  阅读 sklearn/calibration.py 第698-780行，理解 _CalibratedClassifier.predict_proba 的预测流程：

  1. Sigmoid/Isotonic 多分类时如何归一化概率

  2. 当所有校准器输出零概率时，如何回退到均匀分布

  3. Temperature 方法如何调用 softmax 得到概率

  回答问题：

  • 二分类时为什么不需要除分母归一化？

  • 边缘情况：所有校准器输出零概率时使用均匀分布的意义是什么？

• 练习 6：分析测试中 Brier 分数改进验证与标签重映射不变性

  阅读 sklearn/tests/test_calibration.py 第70-120行，理解校准核心测试的验证逻辑：

  1. brier_score_loss(y_test, prob_pos_clf) > brier_score_loss(y_test, prob_pos_cal_clf) 断言的含义

  2. 标签从 [0,1] 映射到 [1,2]、[-1,1]、[1,0] 时，概率输出的一致性验证

  3. Sigmoid 与 Isotonic 在标签反转不变性上的差异

  回答问题：

  • 为什么 Sigmoid 对标签反转满足 prob_new = 1 - prob_old，而 Isotonic 不保证？

  • Brier 分数改进验证为什么是评估校准质量的关键指标？

9.16 本章小结

这一章中我们学习/了解/讨论了概率校准的完整体系。首先，我们理解了为什么需要校准：许多分类器的原始输出并非真实概率，必须通过后处理映射修正。其次，我们深入剖析了 CalibratedClassifierCV 的双模式交叉验证管线：ensemble 模式并行训练多条校准流水线取平均，non-ensemble 模式借用 cross_val_predict 获得无偏预测仅训练单一校准器。接着，我们详细解读了三种校准方法的核心实现：Sigmoid（Platt Scaling）通过 L-BFGS-B 优化带先验的二项对数似然，Isotonic 复用保序回归实现非线性映射，Temperature Scaling 仅用一维标量优化实现多分类原生支持。然后，我们探讨了 _fit_calibrator 工厂函数如何根据方法分拣 OvR 与单校准器两种策略，以及 _CalibratedClassifier 如何封装预测管线并处理归一化与零概率回退。最后，我们了解了校准曲线的分箱统计原理与 CalibrationDisplay 的可视化机制，并分析了测试防线如何从 Brier 分数、标签不变性、并行一致性、Array API 兼容性等多维度守护校准质量。

本章我们一起学习了以下概念：

| 概念 | 解释 |

|------|------|

| CalibratedClassifierCV | 校准主类：通过交叉验证训练基分类器与校准器，支持 ensemble 与 non-ensemble 两种模式 |

| ensemble=True/False | 决定校准器拟合方式：ensemble 训练多条流水线取平均，non-ensemble 用 cross_val_predict 获得无偏预测 |

| _fit_classifier_calibrator_pair | ensemble 模式核心：在 train 子集拟合基分类器，在 test 子集获取原始预测拟合校准器 |

| cross_val_predict | non-ensemble 模式核心：通过交叉验证获得每个样本在未见过折上的无偏预测 |

| _fit_calibrator | 校准器工厂：Sigmoid/Isotonic 使用 OvR 策略，Temperature 天然支持多分类 |

| _SigmoidCalibration | Platt Scaling：假设 sigmoid 关系，用 L-BFGS-B 优化拟合斜率 a_ 和截距 b_ |

| _TemperatureScaling | 温度缩放：单一温度参数 T 控制概率分布平滑度，一维标量优化 log_beta |

| _convert_to_logits | 统一 logits 转换：概率输出用 log(p+eps) 还原，一维决策函数构造 (-x, x) 两列 |

| calibration_curve | 校准曲线计算：将预测概率分箱，统计每箱的平均预测概率与实际正类比例 |

| CalibrationDisplay | 可视化类：绘制可靠性图谱，支持从估计器或原始预测两路构造 |

| _CalibratedClassifier | 管线封装：组合基分类器与校准器列表，提供 predict_proba 实现 OvR 归一化 |

| FrozenEstimator | 冻结已拟合分类器：跳过基分类器训练，所有数据用于校准 |

| Array API 支持 | 仅 Temperature Scaling 支持跨后端：softmax 和标量优化可用 Array API 命名空间实现 |

| get_metadata_routing | 元数据路由：声明 self、estimator、splitter 之间的元数据传递关系 |

| __sklearn_tags__ | 标签系统：标记稀疏输入支持与 Array API 兼容性 |

| test_calibration | 核心测试：Brier 分数改进与标签重映射不变性验证 |

| test_sample_weight | 样本权重测试：验证权重对校准概率的影响 |

| test_parallel_execution | 并行一致性测试：验证 n_jobs=2 与 n_jobs=1 概率输出一致 |

| test_calibration_multiclass | 多分类校准测试：验证概率和为 1 与 Brier 分数改进 |

| test_calibrated_classifier_cv_works_with_large_confidence_scores | 大置信度分数测试：验证 Sigmoid 与 Isotonic 的 ROC AUC 一致性 |

| test_temperature_scaling_array_api_compliance | Array API 合规性测试：验证温度缩放在不同后端的一致性 |

下一章中，我们将学习线性模型大家族 —— 探索“可解释性的主战场”，从普通最小二乘、坐标下降、随机梯度下降到逻辑回归、岭回归与贝叶斯方法，揭示线性方法在工程中的广度与深度。

第 10 章 —— 线性模型大家族 —— 探索“可解释性的主战场”

10.1 学习目标

• 难度：★★★☆☆（3/5）

• 预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础

• 理解 LinearRegression 在 dense/sparse/positive 三种模式下的求解路径与数值稳定性策略

• 掌握 _preprocess_data 与 _rescale_data 对数据中心化、样本权重缩放和稀疏偏移的处理逻辑

• 深入 _pre_fit 中 Gram 矩阵预计算的触发条件、校验机制与 Xy 的计算方式

• 了解 LinearClassifierMixin 与 SparseCoefMixin 如何支撑分类与稀疏系数模型

• 能阅读并分析 LinearModelLoss 的统一损失/梯度/Hessian 接口及其在 GLM 与 Logistic 中的复用

• 掌握 BayesianRidge/ARDRegression 的证据最大化迭代与后验协方差计算

• 了解 HuberRegressor、QuantileRegressor、RANSACRegressor、TheilSenRegressor 的抗离群点机制

10.2 生活类比

想象线性模型是一条精密的数据流水线：_preprocess_data 相当于预处理车间，每个零件（样本）先被称重（sample_weight），再校准基准线（fit_intercept），确保后续加工更精准。稠密数据直接减去加权均值；稀疏数据只能记录均值（X_offset），避免破坏稀疏结构。Gram 预计算相当于事先算好所有零件之间的“配合度矩阵”，后续组装时无需重复测量，显著加速 L1/L0 正则化求解。LinearRegression 相当于最简单的直线装配机，按照最小二乘规则寻找最佳切割角度，支持正系数约束（positive=True）或稀疏求解（lsqr）。LinearModelLoss 相当于统一的误差仪表盘，无论是平方损失还是对数损失，都能输出当前参数下的损失值、梯度和曲率（Hessian），为各种线性模型提供统一的数值后端。BayesianRidge 相当于带“置信度”的装配师，不仅给出最佳切割角度，还给出切割角度的不确定范围（后验协方差），让模型在预测时能够量化不确定性。HuberRegressor 相当于容忍毛刺的质检员，对明显偏离的样本降低惩罚权重，避免单个坏零件毁掉整条产线。QuantileRegressor 相当于分位数裁判，不只关注平均误差，而是确保“90% 的样本落在预测线以下”等概率承诺。RANSACRegressor 相当于随机抽样投票机，反复抽取最小样本集拟合，看哪个模型能获得最多“内点”支持。TheilSenRegressor 相当于鲁棒斜率计算器，对所有成对样本的斜率取空间中位数，天然抗离群点。这一系列组件协同工作，构成了线性模型家族强大且灵活的核心引擎。

10.3 源码地图

sklearn/linear_model/_base.py
├── make_dataset()                        # 稠密/稀疏数据集抽象，稀疏数据 intercept_decay=0.01
├── _preprocess_data()                   # 数据居中、样本权重缩放、稀疏 X_offset 统计
├── _rescale_data()                      # sqrt(sample_weight) 逐样本缩放 X 和 y
├── LinearModel                           # 线性模型基类
│   ├── _decision_function()             # X @ coef + intercept 决策函数
│   ├── predict()                        # 回归预测接口
│   └── _set_intercept()                 # 由 offset 和 coef 反解 intercept_
├── LinearClassifierMixin                # 线性分类器混入
│   ├── decision_function()              # sparse-safe 决策分数
│   ├── predict()                        # argmax 或符号判定预测
│   └── _predict_proba_lr()              # OvR logistic 概率估计
├── SparseCoefMixin                       # coef_ 稀疏化/densify 混入
│   ├── densify()                        # coef_ 转回 dense ndarray
│   └── sparsify()                       # coef_ 转为 CSR 稀疏矩阵
├── LinearRegression                      # 普通最小二乘
│   ├── __init__()                       # 初始化 positive/tol/n_jobs 等
│   └── fit()                            # dense→lstsq / sparse→lsqr / positive→nnls
├── _check_precomputed_gram_matrix()      # 校验用户传入 Gram 矩阵的一致性
└── _pre_fit()                            # L1/L0 模型共用的 Gram/Xy 预计算入口

10.4 线性模型基座与数据预处理 —— 构筑“回归引擎的中央调度室”

10.4.1 关键概念

LinearModel 基类提供了三大核心能力：_decision_function 根据 coef_ 与 intercept_ 计算 X @ coef + intercept，支持稀疏输入；predict 直接复用决策函数完成回归预测；_set_intercept 在 fit_intercept=True 时通过 X_offset 与 y_offset 反解截距，兼容稠密与稀疏情形。make_dataset 为稠密/稀疏数据创建统一的 Dataset 抽象，稀疏数据使用 intercept_decay=0.01 防止截距震荡。数据居中与缩放 _preprocess_data 对稠密数据直接减去加权均值；稀疏数据仅记录均值（X_offset），防止稀疏结构被破坏。若 rescale_with_sw=True，随后调用 _rescale_data 通过 sqrt(sample_weight) 对 X 与 y 进行按样本加权的等价缩放。_pre_fit 在 L1/L0 正则化模型中预计算 Gram 矩阵与 Xy，并在必要时校验用户传入的 Gram 矩阵一致性。

10.4.2 代码解析

flowchart TD A[make_dataset] --> B{X 稀疏?} B -->|是| C[CSRDataset + intercept_decay=0.01] B -->|否| D[ArrayDataset + intercept_decay=1.0] C --> E[返回 dataset, intercept_decay] D --> E

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: make_dataset
# 第 10 章 —— 行号: 65-115
def make_dataset(X, y, sample_weight, random_state=None):
    """Create ``Dataset`` abstraction for sparse and dense inputs."""
    rng = check_random_state(random_state)
    seed = rng.randint(1, np.iinfo(np.int32).max)

    if X.dtype == np.float32:
        CSRData = CSRDataset32
        ArrayData = ArrayDataset32
    else:
        CSRData = CSRDataset64
        ArrayData = ArrayDataset64

    if sp.issparse(X):
        dataset = CSRData(X.data, X.indptr, X.indices, y, sample_weight, seed=seed)
        intercept_decay = SPARSE_INTERCEPT_DECAY  # 0.01
    else:
        X = np.ascontiguousarray(X)
        dataset = ArrayData(X, y, sample_weight, seed=seed)
        intercept_decay = 1.0
    return dataset, intercept_decay

make_dataset 根据输入类型选择 CSRDataset 或 ArrayDataset，稀疏数据使用较小的 intercept_decay 以稳定坐标下降法中的截距更新。

flowchart TD A[输入 X, y, sample_weight] --> B{check_input?} B -->|True| C[validate_data] B -->|False| D[仅类型转换与可选拷贝] C --> E{fit_intercept?} D --> E E -->|True| F{X 是稀疏?} E -->|False| G[X_offset=0, y_offset=0] F -->|是| H[mean_variance_axis 计算 X_offset] F -->|否| I[_average 计算 X_offset 并 X -= X_offset] H --> J[_average 计算 y_offset 并 y -= y_offset] I --> J G --> K[X_scale = 1] J --> K K --> L{sample_weight 且 rescale_with_sw?} L -->|是| M[_rescale_data 缩放 X, y] L -->|否| N[sample_weight_sqrt = None] M --> O[返回 X, y, X_offset, y_offset, X_scale, sample_weight_sqrt] N --> O

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: _preprocess_data
# 第 10 章 —— 行号: 92-188
def _preprocess_data(
    X, y, *, fit_intercept, copy=True, sample_weight=None,
    check_input=True, rescale_with_sw=True,
):
    xp, _, device_ = get_namespace_and_device(X, y, sample_weight)
    n_samples, n_features = X.shape
    X_is_sparse = sp.issparse(X)

    if check_input:
        X = check_array(
            X, copy=copy, accept_sparse=["csr", "csc"],
            dtype=supported_float_dtypes(xp)
        )
        y = check_array(y, dtype=X.dtype, copy=True, ensure_2d=False)
    else:
        y = xp.astype(y, X.dtype)
        if copy:
            if X_is_sparse:
                X = X.copy()
            else:
                X = _asarray_with_order(X, order="K", copy=True, xp=xp)

    if fit_intercept:
        if X_is_sparse:
            X_offset, X_var = mean_variance_axis(X, axis=0, weights=sample_weight)
        else:
            X_offset = _average(X, axis=0, weights=sample_weight, xp=xp)
            X_offset = xp.astype(X_offset, X.dtype, copy=False)
            X -= X_offset

        y_offset = _average(y, axis=0, weights=sample_weight, xp=xp)
        y -= y_offset
    else:
        X_offset = xp.zeros(n_features, dtype=X.dtype, device=device_)
        y_offset = xp.asarray(0.0, dtype=X.dtype, device=device_) if y.ndim == 1 else xp.zeros(y.shape[1], dtype=X.dtype, device=device_)

    X_scale = xp.ones(n_features, dtype=X.dtype, device=device_)

    if sample_weight is not None and rescale_with_sw:
        X, y, sample_weight_sqrt = _rescale_data(X, y, sample_weight, inplace=True)
    else:
        sample_weight_sqrt = None
    return X, y, X_offset, y_offset, X_scale, sample_weight_sqrt

关键点：稀疏数据只记录均值，不做实际减法；随后若需要权重缩放，_rescale_data 会构造对角矩阵 sqrt(sample_weight) 并左乘稀疏矩阵，实现等价的加权最小二乘。

flowchart TD A[输入 X, y, sample_weight] --> B{稀疏?} B -->|是| C[构造对角稀疏矩阵 sw_matrix] B -->|否| D[inplace 判断] C --> E[X = safe_sparse_dot(sw_matrix, X)] C --> F[y = safe_sparse_dot(sw_matrix, y)] D -->|True| G[X *= sample_weight_sqrt[:, None]] D -->|False| H[X = X * sample_weight_sqrt[:, None]] G --> I[返回 X, y, sample_weight_sqrt] H --> I F --> I

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: _rescale_data
# 第 10 章 —— 行号: 193-240
def _rescale_data(X, y, sample_weight, inplace=False):
    xp, _ = get_namespace(X, y, sample_weight)
    n_samples = X.shape[0]
    sample_weight_sqrt = xp.sqrt(sample_weight)

    if sp.issparse(X) or sp.issparse(y):
        sw_matrix = sparse.dia_matrix((sample_weight_sqrt, 0), shape=(n_samples, n_samples))

    if sp.issparse(X):
        X = safe_sparse_dot(sw_matrix, X)
    else:
        if inplace:
            X *= sample_weight_sqrt[:, None]
        else:
            X = X * sample_weight_sqrt[:, None]

    if sp.issparse(y):
        y = safe_sparse_dot(sw_matrix, y)
    else:
        if inplace:
            if y.ndim == 1:
                y *= sample_weight_sqrt
            else:
                y *= sample_weight_sqrt[:, None]
        else:
            if y.ndim == 1:
                y = y * sample_weight_sqrt
            else:
                y = y * sample_weight_sqrt[:, None]
    return X, y, sample_weight_sqrt

稀疏路径：使用 safe_sparse_dot 将对角矩阵左乘稀疏 X，保持稀疏性。

稠密路径：若 inplace=True（仅当在 _preprocess_data 中使用）直接原位乘法，避免额外拷贝。

此函数将加权最小二乘问题 ||S^{1/2}(y - Xw)||^2 转化为普通最小二乘 ||y' - X'w||^2，其中 X' = S^{1/2}X, y' = S^{1/2}y，使得后续求解器无需显式处理样本权重。

flowchart TD A[LinearModel._decision_function] --> B[check_is_fitted] B --> C[validate_data accept_sparse=csr/csc/coo] C --> D{coef_.ndim == 1?} D -->|是| E[return X @ coef_ + intercept_] D -->|否| F[return X @ coef_.T + intercept_]

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: LinearModel._decision_function
# 第 10 章 —— 行号: 270-280
def _decision_function(self, X):
    check_is_fitted(self)
    X = validate_data(self, X, accept_sparse=["csr", "csc", "coo"], reset=False)
    coef_ = self.coef_
    if coef_.ndim == 1:
        return X @ coef_ + self.intercept_
    else:
        return X @ coef_.T + self.intercept_

该方法实现了 X 与系数的矩阵乘法，并在单/多输出情形下分别使用 coef_ 或 coef_.T，确保稀疏矩阵也能高效计算。

flowchart TD A[LinearModel._set_intercept] --> B{fit_intercept?} B -->|否| C[intercept_ = 0.0] B -->|是| D[X_scale 给定?] D -->|是| E[coef_ /= X_scale] D -->|否| F[不缩放] E --> G{coef_.ndim == 1?} F --> G G -->|是| H[intercept_ = y_offset - X_offset @ coef_] G -->|否| I[intercept_ = y_offset - X_offset @ coef_.T]

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: LinearModel._set_intercept
# 第 10 章 —— 行号: 282-298
def _set_intercept(self, X_offset, y_offset, X_scale=None):
    xp, _ = get_namespace(X_offset, y_offset, X_scale)
    if self.fit_intercept:
        self.coef_ = xp.astype(self.coef_, X_offset.dtype, copy=False)
        if X_scale is not None:
            self.coef_ = xp.divide(self.coef_, X_scale)
        if self.coef_.ndim == 1:
            self.intercept_ = y_offset - X_offset @ self.coef_
        else:
            self.intercept_ = y_offset - X_offset @ self.coef_.T
    else:
        self.intercept_ = 0.0

_set_intercept 根据是否进行了中心化（X_offset）以及是否有特征缩放（X_scale）反解截距，兼容单/多输出。

flowchart TD A[_pre_fit] --> B{sparse X?} B -->|是| C[copy=False, precompute=False, rescale_with_sw=False] B -->|否| D[rescale_with_sw=True] C --> E[_preprocess_data] D --> E E --> F{precompute 为数组?} F -->|是| G[check_gram 校验] F -->|否| H[跳过校验] G --> I{fit_intercept 且 X_offset 非零?} I -->|是| J[警告并重置 precompute=auto] H --> K[precompute 自动/显式判断] I --> K J --> K K --> L{precompute is True?} L -->|是| M[分配 precompute 并计算 X.T @ X] L -->|否| N[precompute 非数组时 Xy=None] M --> O{precompute 是数组且 Xy 为 None?} N --> O O -->|是| P[计算 Xy = X.T @ y] O -->|否| Q[保持 Xy] P --> R[返回 X, y, X_offset, y_offset, X_scale, precompute, Xy] Q --> R

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: _pre_fit
# 第 10 章 —— 行号: 615-710
def _pre_fit(
    X, y, Xy, precompute, fit_intercept, copy, check_gram=True, sample_weight=None
):
    n_samples, n_features = X.shape

    if sparse.issparse(X):
        copy = False
        precompute = False
        rescale_with_sw = False
    else:
        rescale_with_sw = True

    X, y, X_offset, y_offset, X_scale, _ = _preprocess_data(
        X, y, fit_intercept=fit_intercept, copy=copy,
        sample_weight=sample_weight, check_input=False,
        rescale_with_sw=rescale_with_sw,
    )

    if hasattr(precompute, "__array__"):
        if fit_intercept and not np.allclose(X_offset, np.zeros(n_features)):
            warnings.warn(
                "Gram matrix was provided but X was centered to fit "
                "intercept: recomputing Gram matrix.", UserWarning)
            precompute = "auto"
            Xy = None
        elif check_gram:
            _check_precomputed_gram_matrix(X, precompute, X_offset, X_scale)

    if isinstance(precompute, str) and precompute == "auto":
        precompute = n_samples > n_features

    if precompute is True:
        precompute = np.empty(shape=(n_features, n_features), dtype=X.dtype, order="C")
        np.dot(X.T, X, out=precompute)

    if not hasattr(precompute, "__array__"):
        Xy = None

    if hasattr(precompute, "__array__") and Xy is None:
        common_dtype = np.result_type(X.dtype, y.dtype)
        if y.ndim == 1:
            Xy = np.empty(shape=n_features, dtype=common_dtype, order="C")
            np.dot(X.T, y, out=Xy)
        else:
            n_targets = y.shape[1]
            Xy = np.empty(shape=(n_features, n_targets), dtype=common_dtype, order="F")
            np.dot(y.T, X, out=Xy.T)

    return X, y, X_offset, y_offset, X_scale, precompute, Xy

_pre_fit 是 Lasso/ElasticNet 等 L1 正则化模型的统一预处理入口。它根据数据稀疏性决定是否预计算 Gram 矩阵，并对用户提供的 Gram 矩阵进行一致性校验（_check_precomputed_gram_matrix）。

flowchart TD A[_check_precomputed_gram_matrix] --> B[选取 f1=floor(n/2), f2=min(f1+1, n-1)] B --> C[计算 v1 = (X[:,f1]-X_off[f1])*X_scale[f1]] C --> D[计算 v2 = (X[:,f2]-X_off[f2])*X_scale[f2]] D --> E[expected = v1 @ v2] E --> F[actual = precompute[f1,f2]] F --> G[rtol 根据 dtype 自适应] G --> H{np.isclose?} H -->|否| I[抛出 ValueError] H -->|是| J[校验通过]

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: _check_precomputed_gram_matrix
# 第 10 章 —— 行号: 555-610
def _check_precomputed_gram_matrix(
    X, precompute, X_offset, X_scale, rtol=None, atol=1e-5
):
    n_features = X.shape[1]
    f1 = n_features // 2
    f2 = min(f1 + 1, n_features - 1)

    v1 = (X[:, f1] - X_offset[f1]) * X_scale[f1]
    v2 = (X[:, f2] - X_offset[f2]) * X_scale[f2]

    expected = np.dot(v1, v2)
    actual = precompute[f1, f2]

    dtypes = [precompute.dtype, expected.dtype]
    if rtol is None:
        rtols = [1e-4 if dtype == np.float32 else 1e-7 for dtype in dtypes]
        rtol = max(rtols)

    if not np.isclose(expected, actual, rtol=rtol, atol=atol):
        raise ValueError(
            "Gram matrix passed in via 'precompute' parameter "
            "did not pass validation when a single element was "
            "checked - please check that it was computed "
            f"properly. For element ({f1},{f2}) we computed "
            f"{expected} but the user-supplied value was "
            f"{actual}."
        )

该函数仅校验 Gram 矩阵的单个非对角元素，以极低开销捕获明显的计算错误。

10.5 小结

这套预处理逻辑为后续求解器（lstsq、lsqr、nnls、cd 等）提供了 统一且数值稳定的输入，并且在稀疏场景下最大限度保持稀疏结构。

10.6 LinearRegression 与分类混入 —— 感受“最小二乘的快与稳”

10.6.1 关键概念

LinearRegression 提供三条求解路径：positive=True 时使用 scipy.optimize.nnls（仅单目标、稠密）；X 稀疏时构造 LinearOperator 并调用 scipy.sparse.linalg.lsqr（atol 与 btol 由 tol 控制）；稠密且非正约束时直接调用 scipy.linalg.lstsq，并根据机器精度设定 cond = max(shape) * eps 防止小奇异值放大误差。稀疏 LSQR 采用自定义的 matvec/rmatvec，在有权重时自动减去均值的影响。LinearClassifierMixin 提供 decision_function（稀疏安全）与 predict（基于符号或 argmax），以及 OvR 多类概率估计实现 _predict_proba_lr。SparseCoefMixin 在 L1 正则化模型（如 Lasso）训练后可将稀疏系数转回 dense（densify）或保持 CSR（sparsify），便于后续预测或模型持久化。

10.6.2 代码解析

flowchart TD A[LinearRegression.fit] --> B{positive?} B -->|True| C[nnls 逐目标并行求解] B -->|False| D{X 稀疏?} D -->|是| E[构造 LinearOperator matvec/rmatvec] E --> F{has_sw?} F -->|是| G[matvec 含 sample_weight_sqrt] F -->|否| H[matvec 仅减 X_offset] G --> I[lsqr 求解 单/多目标并行] H --> I D -->|否| J[lstsq 稠密求解 cond=max(shape)*eps] C --> K[_set_intercept] I --> K J --> K K --> L[返回 self]

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: LinearRegression.fit
# 第 10 章 —— 行号: 465-590
def fit(self, X, y, sample_weight=None):
    n_jobs_ = self.n_jobs
    accept_sparse = False if self.positive else ["csr", "csc", "coo"]

    X, y = validate_data(
        self, X, y, accept_sparse=accept_sparse,
        y_numeric=True, multi_output=True, force_writeable=True,
    )
    has_sw = sample_weight is not None
    if has_sw:
        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype, ensure_non_negative=True)

    copy_X_in_preprocess_data = self.copy_X and not sp.issparse(X)
    X, y, X_offset, y_offset, _, sample_weight_sqrt = _preprocess_data(
        X, y, fit_intercept=self.fit_intercept, copy=copy_X_in_preprocess_data,
        sample_weight=sample_weight,
    )

    if self.positive:
        if y.ndim < 2:
            self.coef_ = optimize.nnls(X, y)[0]
        else:
            outs = Parallel(n_jobs=n_jobs_)(
                delayed(optimize.nnls)(X, y[:, j]) for j in range(y.shape[1])
            )
            self.coef_ = np.vstack([out[0] for out in outs])
    elif sp.issparse(X):
        if has_sw:
            def matvec(b):
                return X.dot(b) - sample_weight_sqrt * b.dot(X_offset)
            def rmatvec(b):
                return X.T.dot(b) - X_offset * b.dot(sample_weight_sqrt)
        else:
            def matvec(b):
                return X.dot(b) - b.dot(X_offset)
            def rmatvec(b):
                return X.T.dot(b) - X_offset * b.sum()
        X_centered = sparse.linalg.LinearOperator(
            shape=X.shape, matvec=matvec, rmatvec=rmatvec
        )
        if y.ndim < 2:
            self.coef_ = lsqr(X_centered, y, atol=self.tol, btol=self.tol)[0]
        else:
            outs = Parallel(n_jobs=n_jobs_)(
                delayed(lsqr)(X_centered, y[:, j].ravel(),
                              atol=self.tol, btol=self.tol) for j in range(y.shape[1]))
            self.coef_ = np.vstack([out[0] for out in outs])
    else:
        cond = max(X.shape) * np.finfo(X.dtype).eps
        self.coef_, _, self.rank_, self.singular_ = linalg.lstsq(X, y, cond=cond)
        self.coef_ = self.coef_.T

    if y.ndim == 1:
        self.coef_ = np.ravel(self.coef_)
    self._set_intercept(X_offset, y_offset)
    return self

• 正系数路径：仅在 positive=True 且 y 为单目标时使用 nnls（非负最小二乘）。

• 稀疏路径：自定义 LinearOperator 能够在 LSQR 中直接考虑 中心化偏置（X_offset）和 样本权重（sample_weight_sqrt），而不需要显式构造稠密矩阵。

• 稠密路径：使用 LAPACK lstsq，cond 与机器 epsilon 结合抑制奇异值导致的数值不稳定。

flowchart TD A[LinearClassifierMixin.decision_function] --> B[check_is_fitted] B --> C[validate_data accept_sparse=csr] C --> D[coef_T = coef_.T if 2D else coef_] D --> E[scores = safe_sparse_dot(X, coef_T) + intercept_] E --> F{scores 为二维且第二维为1?} F -->|是| G[reshape(-1)] F -->|否| H[保持原形状] G --> I[return] H --> I

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: LinearClassifierMixin.decision_function
# 第 10 章 —— 行号: 350-370
def decision_function(self, X):
    check_is_fitted(self)
    xp, _ = get_namespace(X)
    X = validate_data(self, X, accept_sparse="csr", reset=False)
    coef_T = self.coef_.T if self.coef_.ndim == 2 else self.coef_
    scores = safe_sparse_dot(X, coef_T, dense_output=True) + self.intercept_
    return (
        xp.reshape(scores, (-1,))
        if (scores.ndim > 1 and scores.shape[1] == 1)
        else scores
    )

通过 safe_sparse_dot 实现稀疏安全的矩阵乘法，输出 一维 或 二维 置信分数。

flowchart TD A[LinearClassifierMixin._predict_proba_lr] --> B[decision_function 得到 prob] B --> C[_expit(prob) 计算 sigmoid] C --> D{prob.ndim == 1?} D -->|是| E[stack [1-prob, prob]] D -->|否| F[OvR 归一化] F --> G[处理全零概率情况] G --> H[return prob / prob_sum]

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: LinearClassifierMixin._predict_proba_lr
# 第 10 章 —— 行号: 395-415
def _predict_proba_lr(self, X):
    xp, _ = get_namespace(X)
    prob = self.decision_function(X)
    prob = _expit(prob, out=prob, xp=xp)
    if prob.ndim == 1:
        return xp.stack([1 - prob, prob], axis=1)
    else:
        prob_sum = prob.sum(axis=1)
        all_zero = prob_sum == 0
        if xp.any(all_zero):
            prob[all_zero, :] = 1
            prob_sum[all_zero] = prob.shape[1]
        prob /= xp.reshape(prob_sum, (prob.shape[0], -1))
        return prob

实现 OvR (One-vs-Rest) Logistic 回归 的概率估计：二分类直接用 sigmoid，多分类对各类 sigmoid 结果归一化，并处理全零概率的极端情况。

flowchart TD A[SparseCoefMixin.densify] --> B[check_is_fitted] B --> C{coef_ 是稀疏?} C -->|是| D[coef_ = coef_.toarray()] C -->|否| E[无操作] D --> F[return self] E --> F

# 第 10 章 —— 文件: sklearn/linear_model/_base.py
# 第 10 章 —— 方法: SparseCoefMixin.densify / sparsify
# 第 10 章 —— 行号: 420-470
def densify(self):
    msg = "Estimator, %(name)s, must be fitted before densifying."
    check_is_fitted(self, msg=msg)
    if sp.issparse(self.coef_):
        self.coef_ = self.coef_.toarray()
    return self

def sparsify(self):
    msg = "Estimator, %(name)s, must be fitted before sparsifying."
    check_is_fitted(self, msg=msg)
    self.coef_ = sp.csr_matrix(self.coef_)
    return self

densify 与 sparsify 提供稀疏/稠密系数的双向转换，L1 正则化模型训练后可根据下游任务选择内存/速度最优的表示。

10.6.3 小结

LinearRegression 通过 三分支求解 兼顾 正系数约束、稀疏求解 与 数值稳健的稠密求解；分类混入与稀疏系数混入为 线性分类器 与 L1 正则化模型 提供统一且高效的接口。

10.7 LinearModelLoss 与数值核心 —— 解读“损失函数的万能插座”

10.7.1 关键概念

LinearModelLoss 统一封装任意基损失（如 HalfBinomialLoss、HalfMultinomialLoss、HalfPoissonLoss），提供 loss、loss_gradient、gradient、gradient_hessian、gradient_hessian_product 五大入口。系数/截距拆分由 weight_intercept 与 weight_intercept_raw 完成，兼容多类与单类情形。Hessian 核心计算 sandwich_dot 实现 X.T @ diag(W) @ X，在稀疏情况下使用 sparse.dia_matrix 与 safe_sparse_dot，在稠密情况下直接利用 BLAS GEMM。梯度与 Hessian 均除以 样本权重总和（sw_sum），保证 加权平均 的数值一致性。L2 正则化仅加在 权重部分（不含截距）的对角线上。

10.7.2 代码解析

flowchart TD A[LinearModelLoss 统一接口] --> B[base_loss: HalfBinomial/HalfMultinomial/HalfPoisson] A --> C[fit_intercept] A --> D[weight_intercept: 拆分 coef 为 weights + intercept] A --> E[weight_intercept_raw: 同时计算 raw_prediction] A --> F[sandwich_dot: X.T @ diag(W) @ X] F --> G{稀疏?} G -->|是| H[sparse.dia_matrix + safe_sparse_dot] G -->|否| I[W[:, None] * X 然后 X.T @ WX] D --> J[loss/loss_gradient/gradient/gradient_hessian/gradient_hessian_product] J --> K[除以 sw_sum 做加权平均] J --> L[加 L2 正则化到权重对角线]

# 第 10 章 —— 文件: sklearn/linear_model/_linear_loss.py
# 第 10 章 —— 方法: sandwich_dot
# 第 10 章 —— 行号: 20-40
def sandwich_dot(X, W):
    n_samples = X.shape[0]
    if sparse.issparse(X):
        return safe_sparse_dot(
            X.T,
            sparse.dia_matrix((W, 0), shape=(n_samples, n_samples)) @ X,
            dense_output=True,
        )
    else:
        WX = W[:, None] * X
        return X.T @ WX

此函数是 Hessian 计算的瓶颈，对稠密矩阵使用 W[:,None]*X（向量化）避免显式构造对角矩阵，提升缓存友好性；对稀疏矩阵通过 safe_sparse_dot 与 dia_matrix 实现高效乘积。

flowchart TD A[LinearModelLoss.weight_intercept] --> B{is_multiclass?} B -->|否| C{fit_intercept?} C -->|是| D[intercept=coef[-1], weights=coef[:-1]] C -->|否| E[intercept=0, weights=coef] B -->|是| F{coef.ndim == 1?} F -->|是| G[reshape (n_classes, -1) order=F] F -->|否| H[直接使用] G --> I{fit_intercept?} H --> I I -->|是| J[intercept=weights[:,-1], weights=weights[:,:-1]] I -->|否| K[intercept=0]

# 第 10 章 —— 文件: sklearn/linear_model/_linear_loss.py
# 第 10 章 —— 方法: LinearModelLoss.weight_intercept
# 第 10 章 —— 行号: 100-135
def weight_intercept(self, coef):
    if not self.base_loss.is_multiclass:
        if self.fit_intercept:
            intercept = coef[-1]
            weights = coef[:-1]
        else:
            intercept = 0.0
            weights = coef
    else:
        if coef.ndim == 1:
            weights = coef.reshape((self.base_loss.n_classes, -1), order="F")
        else:
            weights = coef
        if self.fit_intercept:
            intercept = weights[:, -1]
            weights = weights[:, :-1]
        else:
            intercept = 0.0
    return weights, intercept

weight_intercept 将扁平化的系数向量拆分为权重与截距，多类时遵循 Fortran-order (列主序) 布局，使得同一特征在所有类别中的系数在内存中连续，加速矩阵-向量乘积。

flowchart TD A[LinearModelLoss.weight_intercept_raw] --> B[调用 weight_intercept 得到 weights, intercept] B --> C{is_multiclass?} C -->|否| D[raw = X @ weights + intercept] C -->|是| E[raw = X @ weights.T + intercept] D --> F[return weights, intercept, raw] E --> F

# 第 10 章 —— 文件: sklearn/linear_model/_linear_loss.py
# 第 10 章 —— 方法: LinearModelLoss.weight_intercept_raw
# 第 10 章 —— 行号: 137-160
def weight_intercept_raw(self, coef, X):
    weights, intercept = self.weight_intercept(coef)
    xp, _, device_ = get_namespace_and_device(X)
    weights_xp = xp.asarray(weights, dtype=X.dtype, device=device_)
    intercept_xp = xp.asarray(intercept, dtype=X.dtype, device=device_)
    if not self.base_loss.is_multiclass:
        raw_prediction = X @ weights_xp + intercept_xp
    else:
        raw_prediction = X @ weights_xp.T + intercept_xp
    return weights, intercept, raw_prediction

同时返回 权重、截距与原始预测值，避免在损失/梯度/Hessian 计算中重复矩阵乘法。

flowchart TD A[gradient_hessian] --> B{raw_prediction 给定?} B -->|否| C[weight_intercept_raw 计算 weights, intercept, raw_prediction] B -->|是| D[weight_intercept 仅拆分] C --> E[base_loss.gradient_hessian 得到 pointwise grad/hess] D --> E E --> F[grad_pointwise /= sw_sum] E --> G[hess_pointwise /= sw_sum] F --> H{多类?} G --> H H -->|否| I[grad[:n_features] = X.T @ grad_pointwise + l2*weights] H -->|是| J[遍历类别 填充块矩阵] I --> K[hess[:n_features, :n_features] = sandwich_dot(X, hess_pointwise)] J --> K K --> L[对角线加 L2 正则化] L --> M{fit_intercept?} M -->|是| N[计算 Xh = X.T @ hess_pointwise 填充截距交叉块] M -->|否| O[返回 grad, hess, warning] N --> O

# 第 10 章 —— 文件: sklearn/linear_model/_linear_loss.py
# 第 10 章 —— 方法: LinearModelLoss.gradient_hessian
# 第 10 章 —— 行号: 290-430
def gradient_hessian(self, coef, X, y, sample_weight=None,
                    l2_reg_strength=0.0, n_threads=1,
                    gradient_out=None, hessian_out=None,
                    raw_prediction=None):
    (n_samples, n_features), n_classes = X.shape, self.base_loss.n_classes
    n_dof = n_features + int(self.fit_intercept)
    if raw_prediction is None:
        weights, intercept, raw_prediction = self.weight_intercept_raw(coef, X)
    else:
        weights, intercept = self.weight_intercept(coef)

    grad_pointwise, hess_pointwise = self.base_loss.gradient_hessian(
        y_true=y, raw_prediction=raw_prediction,
        sample_weight=sample_weight, n_threads=n_threads,
    )
    sw_sum = n_samples if sample_weight is None else np.sum(sample_weight)
    grad_pointwise /= sw_sum
    hess_pointwise /= sw_sum

    if not self.base_loss.is_multiclass:
        if gradient_out is None:
            grad = np.empty_like(coef, dtype=weights.dtype)
        else:
            grad = gradient_out
        grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights
        if self.fit_intercept:
            grad[-1] = grad_pointwise.sum()

        if hessian_out is None:
            hess = np.empty((n_dof, n_dof), dtype=weights.dtype)
        else:
            hess = hessian_out

        hessian_warning = (
            np.average(hess_pointwise <= 0, weights=sample_weight) > 0.25
        )
        hess_pointwise = np.abs(hess_pointwise)

        if hessian_warning:
            return grad, hess, hessian_warning

        hess[:n_features, :n_features] = sandwich_dot(X, hess_pointwise)

        if l2_reg_strength > 0:
            order = "C" if hess.flags.c_contiguous else "F"
            hess.reshape(-1, order=order)[: (n_features * n_dof) : (n_dof + 1)] += l2_reg_strength

        if self.fit_intercept:
            Xh = X.T @ hess_pointwise
            hess[:-1, -1] = Xh
            hess[-1, :-1] = Xh
            hess[-1, -1] = hess_pointwise.sum()
        return grad, hess, hessian_warning
    else:
        # 多分类情形：遍历每个类别，使用 sandwich_dot 填充块矩阵
        # ...

• 点wise 梯度/ Hessian 通过基损失计算后 除以样本权重总和，实现加权均值。

• 稠密情形：sandwich_dot 完成 X' * diag(h) * X，随后在对角线上加入 L2 正则化（仅权重部分）。

• 截距块：通过 Xh = X.T @ hess_pointwise 填充与截距相关的行列，保持 完整二次形式 ([X,1]' @ diag(h) @ [X,1])。

flowchart TD A[gradient_hessian_product] --> B{is_multiclass?} B -->|否| C[weight_intercept_raw 得到 weights, intercept, raw] C --> D[base_loss.gradient_hessian 得到 grad_pt, hess_pt] D --> E[grad = X.T @ grad_pt + l2*weights] E --> F{fit_intercept?} F -->|是| G[grad[-1] = sum(grad_pt)] F -->|否| H[跳过] G --> I[预计算 hX = diag(hess_pt) @ X] I --> J[预计算 hX_sum, hess_sum] J --> K[定义 hessp(s): X.T @ hX @ s[:n_f] + l2*s[:n_f] + s[-1]*hX_sum (含截距项)] K --> L[return grad, hessp] B -->|是| M[base_loss.gradient_proba 得到 grad_pt, proba] M --> N[定义 hessp(s): 基于 proba 计算全 Hessian 向量积] N --> L

# 第 10 章 —— 文件: sklearn/linear_model/_linear_loss.py
# 第 10 章 —— 方法: LinearModelLoss.gradient_hessian_product (单类部分)
# 第 10 章 —— 行号: 480-540
def gradient_hessian_product(self, coef, X, y, sample_weight=None,
                            l2_reg_strength=0.0, n_threads=1):
    # ... 前置代码省略 ...
    grad_pointwise, hess_pointwise = self.base_loss.gradient_hessian(
        y_true=y, raw_prediction=raw_prediction,
        sample_weight=sample_weight, n_threads=n_threads,
    )
    grad_pointwise /= sw_sum
    hess_pointwise /= sw_sum

    grad = np.empty_like(coef, dtype=weights.dtype)
    grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights
    if self.fit_intercept:
        grad[-1] = grad_pointwise.sum()

    hessian_sum = hess_pointwise.sum()
    if sparse.issparse(X):
        hX = sparse.dia_matrix((hess_pointwise, 0), shape=(n_samples, n_samples)) @ X
    else:
        hX = hess_pointwise[:, np.newaxis] * X

    if self.fit_intercept:
        hX_sum = np.squeeze(np.asarray(hX.sum(axis=0)))
        hX_sum = np.atleast_1d(hX_sum)

    def hessp(s):
        ret = np.empty_like(s)
        if sparse.issparse(X):
            ret[:n_features] = X.T @ (hX @ s[:n_features])
        else:
            ret[:n_features] = np.linalg.multi_dot([X.T, hX, s[:n_features]])
        ret[:n_features] += l2_reg_strength * s[:n_features]
        if self.fit_intercept:
            ret[:n_features] += s[-1] * hX_sum
            ret[-1] = hX_sum @ s[:n_features] + hessian_sum * s[-1]
        return ret
    return grad, hessp

gradient_hessian_product 预计算 hX = diag(hess) @ X 及相关汇总统计量，使得 Hessian-向量积 hessp(s) 无需显式构造完整 Hessian 矩阵，极大降低 牛顿法/共轭梯度法 的内存与计算开销。

10.7.3 小结

LinearModelLoss 将 损失 + 正则化 抽象为统一的 数值算子，为 GLM（如 Logistic、Poisson）以及 线性回归 提供可靠的梯度与 Hessian 计算基石。

10.8 贝叶斯线性回归 —— 触摸“会说话的不确定性”

10.8.1 关键概念

BayesianRidge 与 ARDRegression 都采用 证据最大化 迭代估计超参数。BayesianRidge 假设所有权重共享同一精度 lambda_，而 ARDRegression 为每个特征学习独立的精度 lambda_（自动相关性判定）。两者均通过 SVD 分解 加速后验协方差计算：当 n_samples > n_features 时使用 特征空间 的 Vh 进行低维逆；当 n_samples < n_features 时使用 样本空间 的 U 配合 Woodbury 恒等式 仅逆 n_samples × n_samples 矩阵。迭代中通过 γ（有效自由度）更新 alpha_ 与 lambda_，并可选计算 对数边际似然 监控收敛。最终得到 后验均值 coef_ 与 后验协方差 sigma_，支持 predict(return_std=True) 输出预测不确定性。

10.8.2 代码解析

flowchart TD A[BayesianRidge.fit] --> B[_preprocess_data 中心化缩放] B --> C[SVD 分解 X] C --> D{n_samples > n_features?} D -->|是| E[full_matrices=False, Vh (K,N)] D -->|否| F[full_matrices=True, U (M,M), Vh_full (N,N)] E --> G[证据最大化迭代循环] F --> G G --> H[_update_coef_: 后验均值] H --> I[_log_marginal_likelihood 计算证据] I --> J[gamma = sum(alpha*eig / (lambda + alpha*eig))] J --> K[lambda = (gamma + 2*lambda_1) / (||coef||^2 + 2*lambda_2)] K --> L[alpha = (sw_sum - gamma + 2*alpha_1) / (sse + 2*alpha_2)] L --> M{收敛?} M -->|否| G M -->|是| N[最终 _update_coef_] N --> O[计算 sigma_: Vh_full.T @ Vh_full / (alpha*eig_full + lambda)] O --> P[_set_intercept]

# 第 10 章 —— 文件: sklearn/linear_model/_bayes.py
# 第 10 章 —— 方法: BayesianRidge.fit (核心初始化)
# 第 10 章 —— 行号: 200-310
X, y, X_offset_, y_offset_, X_scale_, _ = _preprocess_data(
    X, y, fit_intercept=self.fit_intercept, copy=self.copy_X,
    sample_weight=sample_weight, rescale_with_sw=True,
)
self.X_offset_ = X_offset_
self.X_scale_ = X_scale_

U, S, Vh_full = linalg.svd(X, full_matrices=(n_samples < n_features))
K = len(S)
eigen_vals_ = S**2
eigen_vals_full = np.zeros(n_features, dtype=dtype)
eigen_vals_full[:K] = eigen_vals_
Vh = Vh_full[:K, :]

full_matrices 参数根据样本与特征的大小动态决定，保证 最小维度 的 SVD，从而节省计算资源。

flowchart TD A[_update_coef_] --> B{n_samples > n_features?} B -->|是| C[coef = Vh.T @ (Vh / (eig + lambda/alpha)) @ XT_y] B -->|否| D[coef = X.T @ (U / (eig + lambda/alpha)) @ U.T @ y] C --> E[sse = ||y - X @ coef||^2] D --> E E --> F[return coef, sse]

# 第 10 章 —— 文件: sklearn/linear_model/_bayes.py
# 第 10 章 —— 方法: BayesianRidge._update_coef_
# 第 10 章 —— 行号: 325-355
def _update_coef_(self, X, y, n_samples, n_features, XT_y, U, Vh, eigen_vals_, alpha_, lambda_):
    if n_samples > n_features:
        coef_ = np.linalg.multi_dot(
            [Vh.T,
             Vh / (eigen_vals_ + lambda_ / alpha_)[:, np.newaxis],
             XT_y]
        )
    else:
        coef_ = np.linalg.multi_dot(
            [X.T,
             U / (eigen_vals_ + lambda_ / alpha_)[None, :],
             U.T,
             y]
        )
    sse_ = np.sum((y - np.dot(X, coef_)) ** 2)
    return coef_, sse_

• 当 样本多 时，利用 Vh（特征空间）进行 稀疏逆；

• 当 特征多 时，利用 U（样本空间）与 Woodbury 进行 低维逆。

这种分支策略使得算法在高维（特征远多于样本）与低维场景下均保持高效。

flowchart TD A[_log_marginal_likelihood] --> B{n_samples > n_features?} B -->|是| C[logdet_sigma = -sum(log(lambda + alpha*eig))] B -->|否| D[logdet_sigma = -sum(log(lambda + alpha*eig_full[:n_samples]))] C --> E[score = lambda_1*log(lambda) - lambda_2*lambda + alpha_1*log(alpha) - alpha_2*alpha] D --> E E --> F[score += 0.5*(n_features*log(lambda) + sw_sum*log(alpha) - alpha*sse - lambda*||coef||^2 + logdet_sigma - sw_sum*log(2pi))] F --> G[return score]

# 第 10 章 —— 文件: sklearn/linear_model/_bayes.py
# 第 10 章 —— 方法: BayesianRidge._log_marginal_likelihood
# 第 10 章 —— 行号: 360-395
def _log_marginal_likelihood(self, n_samples, n_features, sw_sum,
                             eigen_vals, alpha_, lambda_, coef, sse):
    if n_samples > n_features:
        logdet_sigma = -np.sum(np.log(lambda_ + alpha_ * eigen_vals))
    else:
        logdet_sigma = np.full(n_features, lambda_, dtype=np.array(lambda_).dtype)
        logdet_sigma[:n_samples] += alpha_ * eigen_vals
        logdet_sigma = -np.sum(np.log(logdet_sigma))

    score = lambda_1 * log(lambda_) - lambda_2 * lambda_
    score += alpha_1 * log(alpha_) - alpha_2 * alpha_
    score += 0.5 * (
        n_features * log(lambda_)
        + sw_sum * log(alpha_)
        - alpha_ * sse
        - lambda_ * np.sum(coef**2)
        + logdet_sigma
        - sw_sum * log(2 * np.pi)
    )
    return score

计算 后验协方差的行列式（logdet_sigma）与 数据拟合误差（sse），形成 证据下界（Marginal Likelihood）供模型选择。超参数更新规则直接源自该目标函数的驻点条件。

flowchart TD A[ARDRegression.fit] --> B[_preprocess_data] B --> C[初始化 alpha_, lambda_ (向量), keep_lambda] C --> D{n_samples >= n_features?} D -->|是| E[_update_sigma: 直接逆 (n_features x n_features)] D -->|否| F[_update_sigma_woodbury: Woodbury 逆 (n_samples x n_samples)] E --> G[迭代循环] F --> G G --> H[update_coeff: coef = alpha * sigma @ X_keep.T @ y] H --> I[计算 sse, gamma = 1 - lambda_keep * diag(sigma)] I --> J[lambda_keep = (gamma + 2*lambda_1) / (coef^2 + 2*lambda_2)] J --> K[alpha = (n_samples - sum(gamma) + 2*alpha_1) / (sse + 2*alpha_2)] K --> L[剪枝: keep_lambda = lambda < threshold_lambda] L --> M{收敛?} M -->|否| G M -->|是| N[最终 update_sigma & update_coeff] N --> O[_set_intercept]

# 第 10 章 —— 文件: sklearn/linear_model/_bayes.py
# 第 10 章 —— 方法: ARDRegression._update_sigma_woodbury
# 第 10 章 —— 行号: 560-575
def _update_sigma_woodbury(self, X, alpha_, lambda_, keep_lambda):
    n_samples = X.shape[0]
    X_keep = X[:, keep_lambda]
    inv_lambda = 1 / lambda_[keep_lambda].reshape(1, -1)
    sigma_ = pinvh(
        np.eye(n_samples, dtype=X.dtype) / alpha_
        + np.dot(X_keep * inv_lambda, X_keep.T)
    )
    sigma_ = np.dot(sigma_, X_keep * inv_lambda)
    sigma_ = -np.dot(inv_lambda.reshape(-1, 1) * X_keep.T, sigma_)
    sigma_[np.diag_indices(sigma_.shape[1])] += 1.0 / lambda_[keep_lambda]
    return sigma_

ARDRegression 在 n_samples < n_features 时使用 Woodbury 公式，仅需逆 n_samples × n_samples 矩阵，避免了对大规模特征协方差矩阵的直接求逆。

flowchart TD A[BayesianRidge.predict / ARDRegression.predict] --> B[_decision_function 得到 y_mean] B --> C{return_std?} C -->|否| D[return y_mean] C -->|是| E[计算 sigmas_squared_data = sum((X @ sigma_) * X, axis=1)] E --> F[y_std = sqrt(sigmas_squared_data + 1/alpha_)] F --> G[return y_mean, y_std]

# 第 10 章 —— 文件: sklearn/linear_model/_bayes.py
# 第 10 章 —— 方法: BayesianRidge.predict (含 return_std)
# 第 10 章 —— 行号: 400-420
def predict(self, X, return_std=False):
    y_mean = self._decision_function(X)
    if not return_std:
        return y_mean
    else:
        sigmas_squared_data = (np.dot(X, self.sigma_) * X).sum(axis=1)
        y_std = np.sqrt(sigmas_squared_data + (1.0 / self.alpha_))
        return y_mean, y_std

预测标准差包含两部分：模型参数不确定性（X @ sigma_ @ X.T 对角线）与 观测噪声（1/alpha_）。

10.8.3 小结

贝叶斯回归通过 SVD/ Woodbury 高效计算后验分布，并在每一步使用 γ 与 证据下界 更新超参数，实现 自动正则化 与 不确定性评估。

10.9 对比四种鲁棒回归器的核心机制 —— 探索“抗噪与保序的特种兵”

下表总结了四种鲁棒回归器的核心机制对比：

| 模型 | 目标 | 关键实现 | 典型使用场景 |

|------|------|----------|--------------|

| HuberRegressor | 同时优化 二次（小残差）与 线性（大残差）损失 | _huber_loss_and_gradient 通过 sigma 统一尺度；使用 L‑BFGS‑B 求解并对 sigma 加下界约束 | 噪声中混有少量强离群点的回归 |

| QuantileRegressor | 最小化 Pinball loss（分位数） | 将问题转化为 线性规划（scipy.optimize.linprog），变量分解为正负 slack；对截距不加惩罚 | 需要显式控制 分位数预测（如风险价值） |

| RANSACRegressor | 在大量离群点中寻找 内点子集 | 随机抽样 → 拟合 → 依据残差阈值划分内外点；动态计算所需最大迭代次数 _dynamic_max_trials | 计算机视觉中 模型估计（如基础矩阵） |

| TheilSenRegressor | 对所有成对斜率取 空间中位数 | 对每个子样本集合求最小二乘 → 通过 _spatial_median 计算中位数；支持随机子采样以控制计算复杂度 | 需要 高鲁棒性 且对 多维斜率 进行稳健估计的回归 |

10.9.1 代码速览（取自各实现）

flowchart TD A[HuberRegressor.fit] --> B[构造初始参数向量 coef+intercept+sigma] B --> C[设置 bounds: sigma >= eps] C --> D[optimize.minimize L-BFGS-B] D --> E[_huber_loss_and_gradient 计算 loss 和 grad] E --> F[outliers_mask = |residual| > epsilon*sigma] F --> G[loss = n_samples*sigma + squared_loss + outlier_loss + alpha*||w||^2] G --> H[grad: squared_part + linear_part + penalty + sigma_grad + intercept_grad] H --> D

# 第 10 章 —— HuberRegressor._fit (核心 L‑BFGS‑B 调用)
opt_res = optimize.minimize(
    _huber_loss_and_gradient,
    parameters,
    method="L-BFGS-B",
    jac=True,
    args=(X, y, self.epsilon, self.alpha, sample_weight),
    options={"maxiter": self.max_iter, "gtol": self.tol, **_get_additional_lbfgs_options_dict("iprint", -1)},
    bounds=bounds,
)

HuberRegressor 将损失函数与梯度打包为单一可调用对象，利用 L-BFGS-B 的边界约束保证 sigma > 0，从而实现自适应尺度估计。bounds 对最后一个参数（sigma）设置下界 eps，防止除零与数值溢出。

flowchart TD A[QuantileRegressor.fit] --> B[过滤零权重样本] B --> C[alpha *= sum(sample_weight)] C --> D{fit_intercept?} D -->|是| E[n_params = n_features + 1] D -->|否| F[n_params = n_features] E --> G[构造目标向量 c: 2*n_params 个 alpha, n_samples 个 quantile*sw, n_samples 个 (1-quantile)*sw] F --> G G --> H[截距对应的 c[0] 和 c[n_params] 置 0] H --> I[构造等式约束 A_eq x = b_eq] I --> J{solver 选择} J -->|highs| K[稀疏 CSC 矩阵构造] J -->|其他| L[稠密矩阵构造] K --> M[linprog 求解] L --> M M --> N[解向量拆分: params = s - t] N --> O[提取 coef_ 和 intercept_]

# 第 10 章 —— QuantileRegressor 线性规划变量构造
c = np.concatenate([np.full(2 * n_params, fill_value=alpha),
                    sample_weight * self.quantile,
                    sample_weight * (1 - self.quantile)])
# 第 10 章 —— 对截距不加惩罚
if self.fit_intercept:
    c[0] = 0
    c[n_params] = 0

QuantileRegressor 将 Pinball 损失转化为 线性规划 标准形式。目标向量 c 中前 2*n_params 项对应系数的正负松弛变量（L1 惩罚），后 2*n_samples 项对应残差的正负松弛（Pinball 损失）。alpha 乘以 sum(sample_weight) 是为了使正则化强度与样本权重一致。对截距不加惩罚（c[0]=c[n_params]=0）体现了截距不应被正则化的建模假设。

flowchart TD A[RANSACRegressor.fit] --> B[确定 min_samples] B --> C[初始化 best_score=-inf] C --> D[循环 n_trials < max_trials] D --> E[随机抽取 min_samples 个样本] E --> F{is_data_valid?} F -->|否| G[n_skips_invalid_data_++] F -->|是| H[estimator.fit 子集] H --> I{is_model_valid?} I -->|否| J[n_skips_invalid_model_++] I -->|是| K[全样本计算残差] K --> L[inlier_mask = 残差 <= residual_threshold] L --> M{n_inliers < n_inliers_best?} M -->|是| N[n_skips_no_inliers_++] M -->|否| O[评分 inlier 子集] O --> P{同 inlier 但 score 更差?} P -->|是| Q[continue] P -->|否| R[更新 best 模型] R --> S[_dynamic_max_trials 更新 max_trials] S --> T{满足 stop_n_inliers 或 stop_score?} T -->|是| U[break] T -->|否| D D --> V[用所有 inlier 重拟合最终模型]

# 第 10 章 —— RANSACRegressor 动态计算所需迭代次数
def _dynamic_max_trials(n_inliers, n_samples, min_samples, probability):
    inlier_ratio = n_inliers / float(n_samples)
    nom = max(_EPSILON, 1 - probability)
    denom = max(_EPSILON, 1 - inlier_ratio**min_samples)
    return abs(float(np.ceil(np.log(nom) / np.log(denom))))

RANSACRegressor 的核心是 随机采样一致性。_dynamic_max_trials 根据当前内点比例动态调整剩余迭代上限：内点越多，所需尝试次数呈指数级下降。stop_probability 与 max_trials 交互决定最终迭代上限，保证以指定置信度至少采样到一次全内点子集。

flowchart TD A[TheilSenRegressor.fit] --> B[_check_subparams 确定 n_subsamples] B --> C{n_choose_k <= max_subpopulation?} C -->|是| D[全组合 combinations] C -->|否| E[随机采样 n_subpopulation 个子集] D --> F[并行 _lstsq 每个子集最小二乘] E --> F F --> G[堆叠所有 weights] G --> H[_spatial_median 迭代 Weiszfeld 算法] H --> I[提取中位数 coefs] I --> J[拆分 intercept_ 与 coef_]

# 第 10 章 —— TheilSenRegressor 子样本最小二乘
weights = _lstsq(X, y, index_list[job], self.fit_intercept)
# 第 10 章 —— 空间中位数求解
n_iter_, coefs = _spatial_median(weights, max_iter=self.max_iter, tol=self.tol)

TheilSenRegressor 对每个子样本集（大小 n_subsamples）求最小二乘解，得到一组系数向量；再对这些向量求 空间中位数（L1 中位数），通过 修正 Weiszfeld 算法 迭代求解。当组合数爆炸时（binom(n, k) > max_subpopulation），退化为随机子采样，以控制计算与内存开销。

flowchart TD A[_spatial_median] --> B{n_features == 1?} B -->|是| C[return 1, median(X)] B -->|否| D[初始化 spatial_median_old = mean(X)] D --> E[循环 max_iter] E --> F[_modified_weiszfeld_step] F --> G{||new - old||^2 < tol?} G -->|是| H[break] G -->|否| I[old = new] I --> E E --> J[返回 n_iter, spatial_median]

# 第 10 章 —— 文件: sklearn/linear_model/_theil_sen.py
# 第 10 章 —— 方法: _spatial_median / _modified_weiszfeld_step
# 第 10 章 —— 行号: 50-110
def _modified_weiszfeld_step(X, x_old):
    diff = X - x_old
    diff_norm = np.sqrt(np.sum(diff**2, axis=1))
    mask = diff_norm >= _EPSILON
    is_x_old_in_X = int(mask.sum() < X.shape[0])
    diff = diff[mask]
    diff_norm = diff_norm[mask][:, np.newaxis]
    quotient_norm = linalg.norm(np.sum(diff / diff_norm, axis=0))
    if quotient_norm > _EPSILON:
        new_direction = np.sum(X[mask, :] / diff_norm, axis=0) / np.sum(1 / diff_norm, axis=0)
    else:
        new_direction = 1.0
        quotient_norm = 1.0
    return (
        max(0.0, 1.0 - is_x_old_in_X / quotient_norm) * new_direction
        + min(1.0, is_x_old_in_X / quotient_norm) * x_old
    )

def _spatial_median(X, max_iter=300, tol=1.0e-3):
    if X.shape[1] == 1:
        return 1, np.median(X.ravel(), keepdims=True)
    tol **= 2
    spatial_median_old = np.mean(X, axis=0)
    for n_iter in range(max_iter):
        spatial_median = _modified_weiszfeld_step(X, spatial_median_old)
        if np.sum((spatial_median_old - spatial_median) ** 2) < tol:
            break
        spatial_median_old = spatial_median
    else:
        warnings.warn("Maximum number of iterations reached in spatial median", ConvergenceWarning)
    return n_iter, spatial_median

_spatial_median 通过 修正 Weiszfeld 算法 迭代求解 L1 中位数。当当前估计落在某个数据点上时（is_x_old_in_X=1），算法会在该点与加权平均方向之间进行凸组合，保证收敛性。

10.9.2 小结

四种鲁棒回归器分别从 损失函数设计、优化问题转化、随机采样共识、空间中位数 四个不同角度解决离群点干扰问题，适用于不同的数据分布与业务需求。

10.10 设计中的取舍

LinearRegression 三路径分离

为什么 LinearRegression 的稀疏路径不支持 positive=True？正系数约束 (nnls) 只能在 稠密矩阵 上求解，因为 scipy.optimize.nnls 需要完整的系数矩阵进行非负最小二乘。稀疏 LSQR 采用迭代求解方式，没有直接的非负约束实现，若强制正系数会导致收敛不确定性。因此设计上选择在稠密情形下提供 positive 支持，而稀疏情形使用标准 LSQR。

LinearRegression 中 tol 参数仅对稀疏 LSQR 有效，为什么不统一到 dense 求解？稠密 lstsq 已经内部使用 机器精度 与 条件数阈值 控制数值误差，额外的容忍度对其影响微乎其微。稀疏 LSQR 采用迭代方法，需要用户显式提供停止准则 atol 与 btol，tol 正好映射到这两个阈值，保持 API 简洁而不破坏已有稠密实现的行为。

稀疏数据的隐式中心化

_preprocess_data 为什么在稀疏情况下不实际中心化 X？对稀疏矩阵直接减去均值会导致 稀疏度急剧下降（几乎所有零元素变为非零），从而失去稀疏优势并显著增加内存需求。_preprocess_data 只记录 X_offset，随后在求解器（如 LSQR）通过自定义 LinearOperator 在乘法阶段 隐式完成中心化，既保持稀疏结构，又实现所需的数学等价性。

LinearModelLoss 统一接口的设计权衡

为什么不直接在基损失类中实现正则化？将 L2 正则化放在 LinearModelLoss 层而非基损失类中，实现了 关注点分离：基损失类只关注单样本的似然函数特性，而线性模型特有的参数结构（权重+截距）与正则化逻辑由上层统一管理。这种设计使得同一基损失可复用于 非线性模型（如 GBM 中的基学习器）而不携带线性模型特有的正则化假设。

BayesianRidge/ARDRegression SVD 与 Woodbury 的选型权衡

BayesianRidge 根据 n_samples > n_features 动态选择 SVD 模式：样本多时用 Vh（特征空间），特征多时用 U + Woodbury（样本空间）。ARDRegression 同理，但在特征多时显式调用 _update_sigma_woodbury。这种 维度自适应 策略避免了对大矩阵的直接求逆，将计算复杂度从 O(min(M,N)^2 * max(M,N)) 降至 O(min(M,N)^3) 级别，是高维贝叶斯推断的关键工程技巧。

HuberRegressor 联合优化 scale 与 coef 的数值稳定性权衡

HuberRegressor 将 sigma 作为优化变量与系数联合优化，而非交替更新。联合优化避免了交替更新可能的振荡，但引入了 非凸性 与 尺度敏感性。为此，bounds 对 sigma 设置严格下界（eps*10），且初始化 sigma=1，配合 L-BFGS-B 的拟牛顿方向搜索，在实践中表现稳健。损失函数中的 n_samples * sigma 项保证了尺度参数有明确的物理意义（M-estimator 的伴随尺度估计）。

QuantileRegressor 线性规划变量拆分的内存/速度权衡

QuantileRegressor 将每个系数拆分为正负两个非负变量（s - t），将每个残差拆分为正负两个松弛变量（u - v），变量总数达到 2*n_params + 2*n_samples。虽然增加了变量维度，但将 非光滑的 Pinball 损失 转化为 标准线性规划，可直接调用高度优化的 scipy.optimize.linprog（HiGHS 求解器），在中小规模问题上通常比次梯度法更快、更稳定。对截距不加惩罚（c[0]=c[n_params]=0）减少了两个变量的正则化压力，符合统计学直觉。

RANSAC 动态迭代与固定迭代的权衡