05 重采样方法 (Resampling Methods)

# 10次不同随机切分，每次计算1-5阶多项式的验证集MSE
fig, ax = plt.subplots(figsize=(10, 5))  # 创建画布
polynomial_degrees = range(1, 6)  # 多项式阶数从1到5

for split_iteration in range(10):  # 循环10次不同的随机切分
    np.random.seed(split_iteration)  # 每次使用不同的随机种子
    total_samples = len(feature_x)  # 总样本数
    shuffled_indices = np.random.permutation(total_samples)  # 生成随机排列的索引
    split_point = total_samples // 2  # 对半分割点
    train_indices = shuffled_indices[:split_point]  # 前一半作为训练集索引
    validation_indices = shuffled_indices[split_point:]  # 后一半作为验证集索引
    mse_per_degree = []  # 存储各阶多项式的MSE
    for degree in polynomial_degrees:  # 遍历每个多项式阶数
        poly_transformer = PolynomialFeatures(degree=degree, include_bias=False)  # 构造多项式特征转换器
        x_train_poly = poly_transformer.fit_transform(feature_x[train_indices])  # 对训练集生成多项式特征
        x_val_poly = poly_transformer.transform(feature_x[validation_indices])  # 对验证集生成同样的多项式特征
        model = LinearRegression().fit(x_train_poly, target_y[train_indices])  # 拟合线性回归模型
        val_predictions = model.predict(x_val_poly)  # 在验证集上预测
        mse_per_degree.append(mean_squared_error(target_y[validation_indices], val_predictions))  # 计算并记录MSE
    ax.plot(list(polynomial_degrees), mse_per_degree, marker='o', alpha=0.6)  # 绘制当次切分的MSE曲线

ax.set_xlabel('多项式阶数 (Polynomial Degree)', fontsize=12)  # X轴标签
ax.set_ylabel('验证集 MSE', fontsize=12)  # Y轴标签
ax.set_title('验证集方法的高方差：10次随机切分结果差异显著', fontsize=14)  # 标题
ax.grid(True, alpha=0.3)  # 添加网格线
plt.tight_layout()  # 调整布局
plt.show()  # 显示图形

# 右图：9次不同随机种子的10折CV
for cv_iteration in range(9):  # 重复9次10折CV
    kfold_mse_per_degree = []  # 存储当次实验各阶MSE
    for degree in range(1, 6):  # 遍历1到5阶多项式
        poly_transformer = PolynomialFeatures(degree=degree, include_bias=False)  # 多项式特征转换器
        x_poly_features = poly_transformer.fit_transform(feature_x)  # 生成多项式特征
        kfold_splitter = KFold(n_splits=10, shuffle=True, random_state=cv_iteration)  # 10折分割器，每次不同随机种子
        neg_mse_scores = cross_val_score(LinearRegression(), x_poly_features, target_y, cv=kfold_splitter, scoring='neg_mean_squared_error')  # 计算10折CV负MSE
        kfold_mse_per_degree.append(-neg_mse_scores.mean())  # 取负得到正MSE
    axes[1].plot(range(1, 6), kfold_mse_per_degree, marker='o', alpha=0.6)  # 绘制当次曲线

axes[1].set_xlabel('多项式阶数', fontsize=12)  # X轴标签
axes[1].set_ylabel('10折 CV MSE', fontsize=12)  # Y轴标签
axes[1].set_title('10折CV（9次不同切分，方差远低于验证集法）', fontsize=14)  # 标题
axes[1].grid(True, alpha=0.3)  # 网格线
plt.tight_layout()  # 调节布局间距
plt.show()  # 显示图表

# 定义最优投资比例α的计算函数（马科维茨最小方差公式）
def compute_optimal_alpha(returns_data, sample_indices):  # 定义计算最优α的函数
    '''根据给定索引的收益率子集，计算使投资组合方差最小的α值'''
    x_returns = returns_data.iloc[sample_indices, 0].values  # 海康威视收益率子集
    y_returns = returns_data.iloc[sample_indices, 1].values  # 宁波银行收益率子集
    variance_x = np.var(x_returns, ddof=1)  # 海康威视收益率方差（无偏估计）
    variance_y = np.var(y_returns, ddof=1)  # 宁波银行收益率方差（无偏估计）
    covariance_xy = np.cov(x_returns, y_returns, ddof=1)[0, 1]  # 两资产协方差
    optimal_alpha = (variance_y - covariance_xy) / (variance_x + variance_y - 2 * covariance_xy)  # 马科维茨最优配比
    return optimal_alpha  # 返回最优投资比例

# 用完整样本计算α的点估计
observation_count = len(sample_returns)  # 样本观测数
original_alpha_estimate = compute_optimal_alpha(sample_returns, range(observation_count))  # 原始样本的α估计
print(f'原始α估计值: {original_alpha_estimate:.4f}')  # 输出α的点估计
print(f'解读: 约{original_alpha_estimate*100:.1f}%配置海康威视，{(1-original_alpha_estimate)*100:.1f}%配置宁波银行')  # 配比解读

# 生成模拟"真实"分布作为对比基准（教学用途：假设已知总体参数）
simulated_alpha_values = []  # 存储模拟分布的α
population_mean = sample_returns.mean()  # 用样本均值近似总体均值
population_cov = sample_returns.cov()  # 用样本协方差近似总体协方差

for _ in range(1000):  # 从"已知总体"中模拟1000次抽样
    simulated_data = np.random.multivariate_normal(population_mean, population_cov, observation_count)  # 从多元正态分布抽取新样本
    simulated_df = pd.DataFrame(simulated_data)  # 转换为DataFrame
    sim_alpha = compute_optimal_alpha(simulated_df, range(observation_count))  # 计算模拟样本的α
    simulated_alpha_values.append(sim_alpha)  # 追加记录

simulated_alpha_values = np.array(simulated_alpha_values)  # 转换为NumPy数组

# 绘制三张对比图
fig, axes = plt.subplots(1, 3, figsize=(18, 5))  # 创建1行3列子图

# 左图：模拟的"真实"抽样分布
axes[0].hist(simulated_alpha_values, bins=40, alpha=0.7, color='#3498DB', edgecolor='black')  # 绘制直方图
axes[0].axvline(np.mean(simulated_alpha_values), color='red', linestyle='--', label=f'均值: {np.mean(simulated_alpha_values):.3f}')  # 标记均值
axes[0].set_xlabel('α', fontsize=12)  # X轴标签
axes[0].set_title('模拟抽样分布\n(上帝视角基准)', fontsize=13)  # 标题
axes[0].legend()  # 图例
axes[0].grid(True, alpha=0.3)  # 网格

# 中图：Bootstrap分布
axes[1].hist(bootstrap_alpha_values, bins=40, alpha=0.7, color='#E74C3C', edgecolor='black')  # 绘制直方图
axes[1].axvline(original_alpha_estimate, color='blue', linewidth=2, label=f'原始估计: {original_alpha_estimate:.3f}')  # 标记原始估计
axes[1].set_xlabel('α', fontsize=12)  # X轴标签
axes[1].set_title(f'Bootstrap分布\n(SE = {bootstrap_se:.3f})', fontsize=13)  # 标题
axes[1].legend()  # 图例
axes[1].grid(True, alpha=0.3)  # 网格

# 右图：箱线图对比
axes[2].boxplot([simulated_alpha_values, bootstrap_alpha_values], labels=['模拟分布', 'Bootstrap'], patch_artist=True, boxprops=dict(facecolor='lightblue'))  # 箱线图
axes[2].set_title('分布对比', fontsize=13)  # 标题
axes[2].grid(True, alpha=0.3, axis='y')  # Y轴网格

plt.tight_layout()  # 调整布局
plt.show()  # 显示图表