04 逻辑回归：从分类问题到概率预测

# 为“一个自然的问题：为什么不用线性回归？”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“一个自然的问题：为什么不用线性回归？”导入 `pandas` 并绑定 `pd`，用于整理当前任务的表格、字段与时间索引。
import pandas as pd
# 为“一个自然的问题：为什么不用线性回归？”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt
# 为“一个自然的问题：为什么不用线性回归？”导入 `seaborn` 并绑定 `sns`，用于编码当前任务的统计分布或分组关系。
import seaborn as sns
# 为“一个自然的问题：为什么不用线性回归？”，从 `sklearn.linear_model` 导入`LinearRegression` 用于拟合普通最小二乘线性回归。
from sklearn.linear_model import LinearRegression

# 生成 DTI 与二元违约标签的机制样本，用于展示线性概率预测越界。
# 为“一个自然的问题：为什么不用线性回归？”，固定随机数序列，使课堂示例可重复。
np.random.seed(42)
# 在给定区间均匀抽取 `dti`，覆盖“一个自然的问题：为什么不用线性回归？”的输入范围。
dti = np.random.uniform(5, 50, 100)
# 将债务收入比映射为逻辑概率，构造随 DTI 单调上升且落在零到一之间的违约概率。
prob_default = 1 / (1 + np.exp(-( -4 + 0.15 * dti)))
# 按给定概率抽取 `is_default` 的 0/1 标签，固定“一个自然的问题：为什么不用线性回归？”的分类结果。
is_default = np.random.binomial(1, prob_default)
# 将债务收入比与二元违约标签整理为分类演示表。
df = pd.DataFrame({'dti': dti, 'is_default': is_default})
# 提取 `X` 作为“一个自然的问题：为什么不用线性回归？”的特征矩阵，明确模型可见的输入列。
X = df[['dti']]
# 提取违约标签，展示 OLS 对二元响应的局限。
y = df['is_default']

# 在二元标签上估计 OLS 线性概率基准，故意保留其范围缺陷供比较。
# 建立 OLS 基准模型，用于展示线性概率预测可能越过 `[0, 1]` 边界。
ols_model = LinearRegression()
# 用全部玩具样本估计 OLS 截距与 DTI 斜率，供线性概率边界反例绘图。
ols_model.fit(X, y)
# 建立 `x_fit` 的有序取值网格，用于展示“线性回归（OLS）直接应用于0/1分类问题的缺陷”随参数变化的比较结果。
x_fit = np.linspace(0, 60, 100).reshape(-1, 1)
# 在覆盖 DTI 范围的网格上生成 OLS 概率预测，检查两端是否越过零或一。
y_fit = ols_model.predict(x_fit)

# 创建 `fig, ax` 画布，承载“线性回归（OLS）直接应用于0/1分类问题的缺陷”的并排视觉比较。
fig, ax = plt.subplots(figsize=(10, 6))
# 以 `'dti'` 为横轴、`'is_default'` 为纵轴绘制散点，数据来自 `df`，展示“线性回归（OLS）直接应用于0/1分类问题的缺陷”。
sns.scatterplot(data=df, x='dti', y='is_default', ax=ax, s=80, alpha=0.7, label='实际数据 (0=未违约, 1=违约)')
# 以 `x_fit` 为横轴、`y_fit` 为纵轴绘制曲线，展示“线性回归（OLS）直接应用于0/1分类问题的缺陷”。
ax.plot(x_fit, y_fit, color='crimson', lw=3, label='OLS 拟合线')

# 添加概率边界、越界标注与统一坐标样式，突出 OLS 反例的失败区域。
# 在 `0` 处添加水平参考线，标出“线性回归（OLS）直接应用于0/1分类问题的缺陷”的基准或阈值。
ax.axhline(0, color='grey', linestyle='--')
# 在 `1` 处添加水平参考线，标出“线性回归（OLS）直接应用于0/1分类问题的缺陷”的基准或阈值。
ax.axhline(1, color='grey', linestyle='--')
# 在图中标注“概率 > 1”，解释“线性回归（OLS）直接应用于0/1分类问题的缺陷”的关键位置。
ax.text(55, 1.05, '概率 > 1', color='red', fontsize=12)
# 在图中标注“概率 < 0”，解释“线性回归（OLS）直接应用于0/1分类问题的缺陷”的关键位置。
ax.text(2, -0.05, '概率 < 0', color='red', fontsize=12)
# 为“一个自然的问题：为什么不用线性回归？”，限定纵轴范围，使比较对象使用一致尺度。
ax.set_ylim(-0.2, 1.2)
# 为“一个自然的问题：为什么不用线性回归？”，限定横轴范围，使比较对象使用一致尺度。
ax.set_xlim(0, 60)
# 将图题设为“线性回归无法将预测值约束在 区间内”，直接说明当前图形的比较目的。
ax.set_title('线性回归无法将预测值约束在 区间内', fontsize=16)
# 将横轴标为“负债收入比 (DTI)”，明确横向编码的变量。
ax.set_xlabel('负债收入比 (DTI)', fontsize=12)
# 把纵轴标签设为 `'预测违约概率 / 实际违约'`。
ax.set_ylabel('预测违约概率 / 实际违约', fontsize=12)
# 显示“一个自然的问题：为什么不用线性回归？”图例，使颜色或线型与比较对象一一对应。
ax.legend()
# 显示二元标签与 OLS 预测线，检查预测值越界及线性概率模型的结构局限。
plt.show()

# 为“逻辑函数能将任意实数映射到 (0, 1) 区间”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“逻辑函数能将任意实数映射到 (0, 1) 区间”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt

# 定义 `sigmoid`，把线性得分转换为 0—1 概率。
def sigmoid(z):
    # 返回 `1 / (1 + np.exp(-z))`，把“逻辑函数能将任意实数映射到 (0, 1) 区间”的计算结果交给调用方。
    return 1 / (1 + np.exp(-z))

# 建立 `z` 的有序取值网格，用于展示“逻辑 (Sigmoid) 函数的S型曲线”随参数变化的比较结果。
z = np.linspace(-10, 10, 200)
# 把线性得分映射为 `g_z` 概率曲线，展示 Sigmoid 的区间压缩作用。
g_z = sigmoid(z)


# 创建 `fig, ax` 画布，承载“逻辑 (Sigmoid) 函数的S型曲线”的并排视觉比较。
fig, ax = plt.subplots(figsize=(10, 6))

# 以 `z` 为横轴、`g_z` 为纵轴绘制曲线，展示“逻辑 (Sigmoid) 函数的S型曲线”。
ax.plot(z, g_z, color='dodgerblue', lw=3)

# 添加零、一与零点五参考线及公式标注，解释逻辑函数的概率边界。
# 在 `0.0` 处添加水平参考线，标出“逻辑 (Sigmoid) 函数的S型曲线”的基准或阈值。
ax.axhline(0.0, color='grey', linestyle='--')
# 在 `1.0` 处添加水平参考线，标出“逻辑 (Sigmoid) 函数的S型曲线”的基准或阈值。
ax.axhline(1.0, color='grey', linestyle='--')
# 在 `0.5` 处添加水平参考线，标出“逻辑 (Sigmoid) 函数的S型曲线”的基准或阈值。
ax.axhline(0.5, color='grey', linestyle=':', lw=1)
# 在 `0.0` 处添加垂直参考线，标出“逻辑 (Sigmoid) 函数的S型曲线”的基准或阈值。
ax.axvline(0.0, color='grey', linestyle=':', lw=1)

# 在图中标注“概率上限 = 1”，解释“逻辑 (Sigmoid) 函数的S型曲线”的关键位置。
ax.text(8, 0.95, '概率上限 = 1', color='grey', fontsize=14)
# 在图中标注“概率下限 = 0”，解释“逻辑 (Sigmoid) 函数的S型曲线”的关键位置。
ax.text(8, 0.05, '概率下限 = 0', color='grey', fontsize=14)
# 在图中标注“g(z) = 0.5”，解释“逻辑 (Sigmoid) 函数的S型曲线”的关键位置。
ax.text(0.2, 0.52, 'g(z) = 0.5', color='black', fontsize=12)
# 在图中标注“$g(z) = \frac{1}{1 + e^{-z}}$”，解释“逻辑 (Sigmoid) 函数的S型曲线”的关键位置。
ax.text(5, 0.6, r'$g(z) = \frac{1}{1 + e^{-z}}$', color='dodgerblue', fontsize=18)

# 隐藏上边框，使视线集中在 Sigmoid 曲线与两条渐近线。
ax.spines['top'].set_visible(False)
# 隐藏右边框，使概率曲线使用更简洁的坐标框。
ax.spines['right'].set_visible(False)
# 用图题概括“逻辑函数能将任意实数映射到 (0, 1) 区间”的比较对象与当前计算结果。
ax.set_title('逻辑函数将 z ∈ (-∞, +∞) 映射到 g(z) ∈ (0, 1)', fontsize=16)
# 将横轴标为“z (线性组合 wᵀx)”，明确横向编码的变量。
ax.set_xlabel('z (线性组合 wᵀx)', fontsize=12)
# 将纵轴标为“g(z) (预测概率)”，明确纵向编码的变量。
ax.set_ylabel('g(z)（模型概率分数）', fontsize=12)
# 显示逻辑函数的 S 形曲线，核对任意实数输入都被压缩到零与一之间。
plt.show()

# 为“案例图解：逻辑函数如何“弯曲”线性关系”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“案例图解：逻辑函数如何“弯曲”线性关系”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt

# 定义 `g`，把线性预测经 Sigmoid 映射为概率。
def g(z):
    # 返回 `1 / (1 + np.exp(-z))`，把“案例图解：逻辑函数如何“弯曲”线性关系”的计算结果交给调用方。
    return 1 / (1 + np.exp(-z))

# 建立 `x` 的有序取值网格，用于展示“线性函数通过逻辑函数转换为概率预测”随参数变化的比较结果。
x = np.linspace(-10, 10, 200)
# 在特征网格上计算类别表面 `z`，用于绘制决策区域。
z = -0.4 + 0.8 * x
# 计算线性得分 `f_x`，作为 Sigmoid 转换前的模型输出。
f_x = g(z)


# 创建 `fig, ax` 画布，承载“线性函数通过逻辑函数转换为概率预测”的并排视觉比较。
fig, ax = plt.subplots(figsize=(10, 6))

# 绘出未压缩的线性得分，作为 Sigmoid 概率映射的对照。
ax.plot(x, z, color='grey', lw=2, linestyle='--', label=r'线性部分: $z = -0.4 + 0.8x$')
# 以 `x` 为横轴、`f_x` 为纵轴绘制曲线，展示“线性函数通过逻辑函数转换为概率预测”。
ax.plot(x, f_x, color='crimson', lw=3, label=r'逻辑回归输出: $f(x) = g(z)$')

# 添加线性得分、概率阈值和对应点标注，连接两个坐标面板的映射。
# 在 `0.0` 处添加水平参考线，标出“线性函数通过逻辑函数转换为概率预测”的基准或阈值。
ax.axhline(0.0, color='darkgrey', linestyle='-', lw=1)
# 在 `1.0` 处添加水平参考线，标出“线性函数通过逻辑函数转换为概率预测”的基准或阈值。
ax.axhline(1.0, color='darkgrey', linestyle='-', lw=1)
# 在 `0.0` 处添加垂直参考线，标出“线性函数通过逻辑函数转换为概率预测”的基准或阈值。
ax.axvline(0.0, color='darkgrey', linestyle='-', lw=1)

# 求出概率曲线与 0.5 阈值的交点 `x_intersect`，定位分类决策边界。
x_intersect = 0.5
# 以 `x_intersect` 为横轴、`0.5` 为纵轴绘制曲线，展示“线性函数通过逻辑函数转换为概率预测”。
ax.plot(x_intersect, 0.5, 'bo', markersize=8, zorder=5)
# 在 `ax` 上画出指定横坐标和纵向范围的竖线。
ax.vlines(x_intersect, -4, 0.5, color='blue', linestyle=':', lw=2)
# 在 `ax` 上画出指定纵坐标和横向范围的横线。
ax.hlines(0.5, -10, x_intersect, color='blue', linestyle=':', lw=2)
# 固定 `ax.text(x_intersect + 0.2, 0.6, f'决策边界: x` 的业务取值与顺序，作为“案例图解：逻辑函数如何“弯曲”线性关系”的可复算输入。
ax.text(x_intersect + 0.2, 0.6, f'决策边界: x = {x_intersect:.1f}\n预测概率 = 0.5', color='blue')

# 将图题设为“逻辑函数保持了单调性，同时将输出约束在概率区间”，直接说明当前图形的比较目的。
ax.set_title('逻辑函数保持了单调性，同时将输出约束在概率区间', fontsize=16)
# 将横轴标为“特征 x”，明确横向编码的变量。
ax.set_xlabel('特征 x', fontsize=12)
# 将纵轴标为“输出值”，明确纵向编码的变量。
ax.set_ylabel('输出值', fontsize=12)
# 显示“案例图解：逻辑函数如何“弯曲”线性关系”图例，使颜色或线型与比较对象一一对应。
ax.legend(fontsize=12)
# 为“案例图解：逻辑函数如何“弯曲”线性关系”，限定纵轴范围，使比较对象使用一致尺度。
ax.set_ylim(-4, 4)
# 显示线性预测子与逻辑概率的对应关系，检查阈值附近最陡、两端渐近饱和的形状。
plt.show()

# 为“决策边界的可视化”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“决策边界的可视化”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt
# 为“决策边界的可视化”，从 `sklearn.datasets` 导入 `make_classification`，取得本节机制演示的数据接口。
from sklearn.datasets import make_classification
# 为“决策边界的可视化”，从 `sklearn.linear_model` 导入`LogisticRegression` 用于拟合逻辑回归分类器。
from sklearn.linear_model import LogisticRegression

# 生成二维分类机制样本，为逻辑概率表面和决策边界提供输入。
# 生成二维可分机制样本，专门用于展示逻辑回归决策边界的几何形状。
X, y = make_classification(n_samples=100, n_features=2, n_informative=2, n_redundant=0,
                           # 每个类别只生成一个簇，使二维决策边界保持清晰可辨。
                           n_clusters_per_class=1, flip_y=0.1, random_state=1)

# 用二维样本拟合逻辑回归，估计绘制概率区域所需的系数。
# 创建 `LogisticRegression` 实例 `model`，用于拟合逻辑回归分类器。
model = LogisticRegression(solver='liblinear')
# 用二维玩具样本估计逻辑回归系数，供网格概率与 `p=0.5` 边界计算。
model.fit(X, y)

# 在两项特征范围内构造细网格，用于计算连续正类概率表面。
# 用第一特征的样本极值外扩一单位，确定决策面横向绘图区间。
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
# 用第二特征的样本极值外扩一单位，确定决策面纵向绘图区间。
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
# 建立 `xx, yy` 的有序取值网格，用于展示“一个二维特征空间中的线性决策边界”随参数变化的比较结果。
xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02), np.arange(y_min, y_max, 0.02))
# 计算网格点属于正类的概率并提取第二列，形成决策区域的扁平概率数组。
Z = model.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
# 在特征网格上计算类别表面 `Z`，用于绘制决策区域。
Z = Z.reshape(xx.shape)

# 创建 `fig, ax` 画布，承载“一个二维特征空间中的线性决策边界”的并排视觉比较。
fig, ax = plt.subplots(figsize=(12, 6))
# 保存决策区域图层 `contour`，随后为类别背景添加颜色标尺。
contour = ax.contourf(xx, yy, Z, cmap='RdBu_r', alpha=0.6)
# 为当前颜色编码添加数值色标。
fig.colorbar(contour, ax=ax, label='预测为正类 (y=1) 的概率')
# 用 `ax.contour` 绘制预测类别边界的等值线。
ax.contour(xx, yy, Z, levels=[0.5], colors='black', linewidths=2)
# 保留 `scatter` 的绘图句柄，用于构造与类别颜色一致的图例。
negative = y == 0
ax.scatter(X[negative, 0], X[negative, 1], c='#2C7BB6', marker='o', edgecolors='k', s=90, label='负类（y=0）')
ax.scatter(X[~negative, 0], X[~negative, 1], c='#D7191C', marker='^', edgecolors='k', s=100, label='正类（y=1）')
# 颜色与点形共同区分类别，方便投影和灰度阅读。
ax.legend()
# 将图题设为“决策边界 (wᵀx=0) 分割了预测区域”，直接说明当前图形的比较目的。
ax.set_title('决策边界 (wᵀx=0) 分割了预测区域', fontsize=16)
# 使用普通数字，避免不同电脑缺少 Unicode 下标字形。
ax.set_xlabel('特征 1', fontsize=14)
ax.set_ylabel('特征 2', fontsize=14)
# 在图中标注“预测为 负类 (p < 0.5)”，解释“一个二维特征空间中的线性决策边界”的关键位置。
ax.text(-2, -2.5, '预测为 负类 (p < 0.5)', ha='center', fontsize=12, color='white')
# 在图中标注“预测为 正类 (p > 0.5)”，解释“一个二维特征空间中的线性决策边界”的关键位置。
ax.text(1.5, 2, '预测为 正类 (p > 0.5)', ha='center', fontsize=12, color='white')
# 显示二维样本、概率背景与 `p=0.5` 决策边界，检查两类区域的划分及边界附近不确定性。
plt.show()

# 为“步骤 1: 导入库与数据准备”导入 `pandas` 并绑定 `pd`，用于整理当前任务的表格、字段与时间索引。
import pandas as pd
from pathlib import Path
from urllib.request import urlretrieve  # 复用本章隐藏设置单元安装的浏览器标识下载器
# 导入 statsmodels 并绑定为 `sm`，用于给设计矩阵加截距并估计可报告系数检验的 Logit 模型。
import statsmodels.api as sm

# 所有读者都从同一公开地址取得课程数据；下载一次后可离线重复运行。
# 按 Linux 共享数据、Windows 共享数据、项目缓存的顺序选择审计意见文件。
audit_path = next((candidate_path for candidate_path in [Path('/home/ubuntu/r2_data_mount/data/stock/audit_opinion.h5'), Path('C:/qiufei/data/stock/audit_opinion.h5'), Path('data/course/audit_opinion.h5')] if candidate_path.exists()), Path('data/course/audit_opinion.h5'))
if not audit_path.exists():
    audit_path.parent.mkdir(parents=True, exist_ok=True)
    urlretrieve('https://assets.qiufei.site/data/stock/audit_opinion.h5', audit_path)
# 读取审计意见表，保留全部披露版本供按时间筛选。
audit_raw = pd.read_hdf(audit_path, key='audit_opinion')
# 把披露日期统一为时间戳，以便计算年末至审计报告日的滞后。
audit_raw['info_date'] = pd.to_datetime(audit_raw['info_date'])
# 只保留 2014—2023 年度财报审计记录，确定本节样本期与报告类型。
audit_raw = audit_raw.query("type == 'financial_statements' and quarter.str.endswith('q4') and '2014q4' <= quarter <= '2023q4'", engine='python').copy()
# 统计存在多个披露版本的公司—财年键，量化重发污染风险。
duplicate_company_years = int(audit_raw.groupby(['order_book_id', 'quarter']).size().gt(1).sum())
# 按披露日选择每个公司—财年的首次可得版本，禁止后续修订回填同期样本。
data = audit_raw.sort_values(['order_book_id', 'quarter', 'info_date']).drop_duplicates(['order_book_id', 'quarter'], keep='first').copy()
# 检查点时观察单位唯一，防止同一公司财年重复加权。
assert not data.duplicated(['order_book_id', 'quarter']).any(), '首次披露数据版本公司—财年键不唯一'
# 从季度字段提取财年，作为审计制度与披露环境的时间控制量。
data['fiscal_year'] = data['quarter'].str[:4].astype(int)
# 计算财年末至披露日的天数，描述当前审计报告的披露滞后。
data['report_delay_days'] = (data['info_date'] - pd.to_datetime(data['fiscal_year'].astype(str) + '-12-31')).dt.days
# 标记常见国际四大会计师事务所在中国的名称，用作可解释的机构类别变量。
data['is_big4_cn'] = data['audit_agency'].fillna('').str.contains('普华永道|德勤|毕马威|安永').astype(int)
# 将标准无保留意见编码为零、其他意见编码为一，固定正类业务含义。
data['modified_opinion'] = data['opinion_type'].ne('unqualified').astype(int)
# 固定披露滞后、财年与机构类别三项解释变量，限定本节教学规格。
features = ['report_delay_days', 'fiscal_year', 'is_big4_cn']
# 将非标准意见标记固定为二分类目标列。
target = 'modified_opinion'
# 报告版本筛选流量与异常长披露滞后，明确当前同期描述的样本边界。
print({'原始版本行': len(audit_raw), '重复公司财年键': duplicate_company_years, '首次披露行': len(data), '移除后续版本': len(audit_raw) - len(data), '披露滞后超过550天': int(data['report_delay_days'].gt(550).sum())})

# 提取 `X` 作为“步骤 1: 导入库与数据准备”的特征矩阵，明确模型可见的输入列。
X = data[features]
# 提取已重编码的审计意见目标，其中 1 表示非标准意见类别。
y = data[target]

# 给三项特征添加常数列，使 statsmodels 显式估计 Logit 截距。
X_sm = sm.add_constant(X)

# 为“步骤 1 (续): 数据预览”，显示前几条记录以核对字段与取值。
X_sm.head()

# 为“步骤 1 (续): 数据预览”，显示前几条记录以核对字段与取值。
y.head()

# 构建 Logit 模型
logit_model = sm.Logit(y, X_sm)

# 按公司代码聚类标准误，处理同一公司跨财年的相关性而不把记录误作独立横截面。
result = logit_model.fit(cov_type='cluster', cov_kwds={'groups': data['order_book_id']})

# 将系数、公司聚类标准误、p 值与置信区间压成投影可读表。
logit_inference = pd.concat([result.params.rename('coef'), result.bse.rename('cluster_std_err'), result.pvalues.rename('p_value'), result.conf_int().rename(columns={0: 'ci_low', 1: 'ci_high'})], axis=1)
# 展示本次执行生成的唯一系数证据，避免默认摘要滚动和旧案例数字残留。
display(logit_inference.round(6))

# 计算连续变量平均导数，并把二元变量按 0→1 平均离散变化处理。
marginal_effects = result.get_margeff(at='overall', method='dydx', dummy=True)
# 提取 statsmodels 原始命名表，保留所有六项统计量供逐列核对。
ame_source = marginal_effects.summary_frame()
# 显式声明统计量映射，避免按位置错把 z 当作 p 值或丢失区间上界。
ame_column_map = {'dy/dx': 'AME', 'Std. Err.': '聚类稳健标准误', 'z': 'z值', 'Pr(>|z|)': 'p值', 'Conf. Int. Low': '95%CI下界', 'Cont. Int. Hi.': '95%CI上界'}
# 在展示前检查当前 statsmodels 版本确实提供事先确定的六列。
assert set(ame_column_map).issubset(ame_source.columns), f'AME列名变化: {ame_source.columns.tolist()}'
# 按明确列名选择并重命名变量级 AME 推断表。
ame_table = ame_source[list(ame_column_map)].rename(columns=ame_column_map)
# 展示与当前审计意见模型同源的动态 AME 证据。
display(ame_table.round(6))