示例

mport numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import warnings

warnings.filterwarnings('ignore')


def fast_vif_calculation(X_train, threshold=5, max_features=100, sample_size=1000):
    """
    快速计算VIF（方差膨胀因子）

    策略：
    1. 如果特征数 <= max_features，直接计算所有特征的VIF
    2. 如果特征数 > max_features，使用PCA降维后计算，再映射回原始特征

    Parameters:
    -----------
    X_train : DataFrame
        训练数据（仅包含数值特征）
    threshold : float
        VIF阈值
    max_features : int
        最大特征数限制
    sample_size : int
        采样大小

    Returns:
    --------
    dict : 特征名 -> VIF值的字典
    """

    # 处理采样
    if len(X_train) > sample_size:
        X_sample = X_train.sample(n=sample_size, random_state=42)
    else:
        X_sample = X_train.copy()

    # 标准化数据
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X_sample)
    X_scaled = pd.DataFrame(X_scaled, columns=X_sample.columns, index=X_sample.index)

    # 处理常数特征
    constant_cols = []
    for col in X_scaled.columns:
        if X_scaled[col].std() < 1e-10:
            constant_cols.append(col)

    if constant_cols:
        print(f"  ⚠️ 发现 {len(constant_cols)} 个常数特征，计算VIF时会特殊处理")
        X_scaled = X_scaled.drop(columns=constant_cols)

    # 如果特征数太少，直接返回
    if X_scaled.shape[1] < 2:
        vif_dict = {col: 1.0 for col in X_train.columns}
        return vif_dict

    # 根据特征数选择计算策略
    n_features = X_scaled.shape[1]

    if n_features <= max_features:
        # 策略1：直接计算VIF
        vif_dict = calculate_vif_direct(X_scaled)
    else:
        # 策略2：使用PCA加速计算
        vif_dict = calculate_vif_pca_based(X_scaled, n_components=min(max_features, n_features // 2))

    # 处理常数特征（VIF设为无穷大或极大值）
    for col in constant_cols:
        vif_dict[col] = 1e10

    return vif_dict


def calculate_vif_direct(X_scaled):
    """直接计算VIF（适用于特征数较少的情况）"""
    from statsmodels.stats.outliers_influence import variance_inflation_factor

    vif_dict = {}
    n_features = X_scaled.shape[1]

    # 检查是否存在完全共线性
    try:
        # 计算相关矩阵的秩
        corr_matrix = X_scaled.corr().values
        rank = np.linalg.matrix_rank(corr_matrix)

        if rank < n_features:
            print(f"  ⚠️ 存在完全共线性，使用改进方法计算VIF")
            # 使用改进方法
            return calculate_vif_with_ridge(X_scaled)

        # 正常计算VIF
        for i in range(n_features):
            vif = variance_inflation_factor(X_scaled.values, i)
            vif_dict[X_scaled.columns[i]] = vif

    except Exception as e:
        print(f"  ⚠️ 标准VIF计算失败，使用岭回归方法: {str(e)}")
        return calculate_vif_with_ridge(X_scaled)

    return vif_dict


def calculate_vif_with_ridge(X_scaled, alpha=0.1):
    """使用岭回归计算近似VIF（处理共线性问题）"""
    from sklearn.linear_model import Ridge
    from sklearn.metrics import r2_score

    vif_dict = {}
    n_features = X_scaled.shape[1]

    for i in range(n_features):
        y = X_scaled.iloc[:, i]
        X = X_scaled.drop(X_scaled.columns[i], axis=1)

        # 使用岭回归
        ridge = Ridge(alpha=alpha)
        ridge.fit(X, y)

        # 预测并计算R²
        y_pred = ridge.predict(X)
        r2 = r2_score(y, y_pred)

        # 计算VIF
        if r2 < 1:
            vif = 1 / (1 - r2)
        else:
            vif = 1e10  # 完全共线性

        vif_dict[X_scaled.columns[i]] = vif

    return vif_dict


def calculate_vif_pca_based(X_scaled, n_components=None):
    """
    基于PCA的VIF近似计算（适用于特征数很多的情况）

    原理：使用PCA降维后，计算每个原始特征在主要成分上的投影，
    近似估计其VIF值
    """
    from sklearn.decomposition import PCA

    if n_components is None:
        n_components = min(X_scaled.shape[1] // 2, 50)

    # 确保n_components不超过特征数
    n_components = min(n_components, X_scaled.shape[1])

    # 拟合PCA
    pca = PCA(n_components=n_components)
    pca.fit(X_scaled)

    # 计算每个特征在各主成分上的权重
    components = pca.components_
    explained_variance = pca.explained_variance_ratio_

    # 近似计算每个特征的VIF
    vif_dict = {}
    n_features = X_scaled.shape[1]

    for i in range(n_features):
        # 特征在主成分上的投影权重
        weights = components[:, i]

        # 计算该特征解释的方差比例（近似）
        # 使用权重的平方和乘以解释方差
        variance_explained = np.sum((weights ** 2) * explained_variance)

        # 近似VIF = 1 / (1 - variance_explained)
        if variance_explained < 0.999:
            vif = 1 / (1 - variance_explained)
        else:
            vif = 1e10

        vif_dict[X_scaled.columns[i]] = vif

    return vif_dict


def calculate_vif_eigen_based(X_scaled):
    """
    基于特征值的VIF计算方法（更精确但计算量大）
    适用于中等规模数据
    """
    # 计算相关矩阵
    corr_matrix = X_scaled.corr().values

    # 计算特征值
    eigenvalues = np.linalg.eigvalsh(corr_matrix)

    # 条件数检查
    condition_number = np.max(eigenvalues) / np.min(eigenvalues[eigenvalues > 1e-10])

    if condition_number > 1000:
        print(f"  ⚠️ 条件数过大 ({condition_number:.2e})，存在严重共线性")

    # 计算每个特征的VIF（基于特征值分解）
    vif_dict = {}
    n_features = X_scaled.shape[1]

    for i in range(n_features):
        # 计算第i个特征的VIF
        # VIF_i = 1 / (1 - R²_i)
        # 可以用相关矩阵的逆矩阵的对角元素来计算
        try:
            corr_inv = np.linalg.inv(corr_matrix)
            vif = corr_inv[i, i]
            vif_dict[X_scaled.columns[i]] = vif
        except np.linalg.LinAlgError:
            # 矩阵奇异，使用伪逆
            corr_pinv = np.linalg.pinv(corr_matrix)
            vif = corr_pinv[i, i]
            vif_dict[X_scaled.columns[i]] = vif

    return vif_dict


# 完整的 fast_remove_multicollinearity 函数（改进版）
def fast_remove_multicollinearity(X_train, X_test, threshold=5,
                                  max_features=100, sample_size=1000):
    """快速删除共线性特征"""
    print(f"\n  🔍 快速删除共线性特征 (阈值={threshold})...")
    original_features = X_train.shape[1]

    # 确保只处理数值特征
    X_train_num = X_train.select_dtypes(include=[np.number])
    X_test_num = X_test.select_dtypes(include=[np.number])

    if X_train_num.shape[1] < 2:
        print("  ⚠️ 数值特征不足2个，跳过处理")
        return X_train, X_test

    # 处理样本大小
    actual_sample_size = min(sample_size, len(X_train))
    if actual_sample_size < len(X_train):
        print(f"  使用 {actual_sample_size} 个样本进行计算")
        sampled_data = X_train_num.sample(n=actual_sample_size, random_state=42)
    else:
        sampled_data = X_train_num

    # 步骤1：相关系数筛选（改进）
    print(f"  步骤1: 相关系数筛选...")
    corr_matrix = sampled_data.corr().abs()
    upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))

    to_drop_corr = set()
    for col in upper.columns:
        high_corr_cols = []
        for idx in upper.index:
            if upper.loc[idx, col] > 0.95 and idx != col:
                high_corr_cols.append(idx)

        for high_corr_col in high_corr_cols:
            # 保留方差大的特征
            if X_train_num[col].var() < X_train_num[high_corr_col].var():
                to_drop_corr.add(col)
            else:
                to_drop_corr.add(high_corr_col)

    if to_drop_corr:
        print(f"  删除 {len(to_drop_corr)} 个高相关特征")
        X_train = X_train.drop(columns=to_drop_corr)
        X_test = X_test.drop(columns=to_drop_corr)
        X_train_num = X_train.select_dtypes(include=[np.number])
        X_test_num = X_test.select_dtypes(include=[np.number])

    # 步骤2：VIF筛选
    print(f"  步骤2: VIF筛选...")
    iteration = 1
    max_iterations = min(50, X_train_num.shape[1])

    while True:
        if X_train_num.shape[1] < 2:
            break

        vif_dict = fast_vif_calculation(X_train_num, threshold, max_features, actual_sample_size)

        # 过滤出超过阈值的特征
        high_vif_features = {col: vif for col, vif in vif_dict.items() if vif > threshold}

        if not high_vif_features:
            print(f"  ✅ 所有特征VIF <= {threshold}")
            break

        # 删除VIF最高的特征
        drop_col = max(high_vif_features.items(), key=lambda x: x[1])[0]
        X_train = X_train.drop(columns=[drop_col])
        X_test = X_test.drop(columns=[drop_col])
        X_train_num = X_train.select_dtypes(include=[np.number])
        X_test_num = X_test.select_dtypes(include=[np.number])

        print(f"  删除特征: {drop_col} (VIF={high_vif_features[drop_col]:.2f})")

        iteration += 1
        if iteration > max_iterations:
            print(f"  ⚠️ 达到最大迭代次数 {max_iterations}")
            break

    total_removed = original_features - X_train.shape[1]
    print(f"  ✅ 共删除 {total_removed} 个特征，保留 {X_train.shape[1]} 个")

    return X_train, X_test


def remove_multicollinearity(X_train, X_test, threshold=5, method='fast_vif',
                             max_features=100, sample_size=1000):
    """删除共线性特征（主函数）"""

    if method == 'correlation':
        print(f"\n  使用相关系数删除共线性特征...")
        corr_matrix = X_train.corr().abs()
        upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))
        to_drop = [col for col in upper.columns if any(upper[col] > 0.95)]
        if to_drop:
            print(f"  删除 {len(to_drop)} 个高相关特征")
            X_train = X_train.drop(columns=to_drop)
            X_test = X_test.drop(columns=to_drop)
        return X_train, X_test

    elif method == 'vif':
        print(f"\n  🔍 使用标准VIF删除共线性特征...")
        X_train_num = X_train.select_dtypes(include=[np.number]).copy()
        if X_train_num.isnull().any().any():
            print("  ⚠️ 存在缺失值，请先处理")
            return X_train, X_test
        if X_train_num.shape[1] < 2:
            return X_train, X_test

        from statsmodels.stats.outliers_influence import variance_inflation_factor

        dropped_features = []
        while True:
            vif_df = pd.DataFrame({
                'feature': X_train_num.columns,
                'VIF': [variance_inflation_factor(X_train_num.values, i)
                        for i in range(X_train_num.shape[1])]
            })
            max_vif = vif_df['VIF'].max()
            if max_vif <= threshold:
                break
            drop_feature = vif_df.loc[vif_df['VIF'].idxmax(), 'feature']
            dropped_features.append(drop_feature)
            X_train_num = X_train_num.drop(drop_feature, axis=1)
            print(f"  删除 '{drop_feature}' (VIF={max_vif:.2f})")

            if X_train_num.shape[1] < 2:
                break

        X_train_final = X_train.drop(dropped_features, axis=1)
        X_test_final = X_test.drop(dropped_features, axis=1)
        print(f"  保留 {X_train_final.shape[1]} 个特征")
        return X_train_final, X_test_final

    else:  # fast_vif
        return fast_remove_multicollinearity(X_train, X_test, threshold, max_features, sample_size)

调用示例

from sklearn.model_selection import train_test_split



# 生成示例数据（包含多重共线性）
def create_multicollinear_data(n_samples=1000, n_features=50, random_state=42):
    """
    创建包含多重共线性的数据集
    """
    np.random.seed(random_state)

    # 基础特征
    X_base = np.random.randn(n_samples, 20)

    # 创建高度相关的特征（共线性）
    X_corr = []
    for i in range(15):
        # 每个新特征与前面的特征高度相关
        base_idx = i % 10
        noise = np.random.randn(n_samples) * 0.1  # 小噪声
        new_feature = X_base[:, base_idx] * 0.9 + noise
        X_corr.append(new_feature.reshape(-1, 1))

    X_corr = np.hstack(X_corr)

    # 创建一些随机噪声特征
    X_noise = np.random.randn(n_samples, 15)

    # 合并所有特征
    X = np.hstack([X_base, X_corr, X_noise])

    # 创建目标变量（只与前10个特征相关）
    true_coef = np.zeros(X.shape[1])
    true_coef[:10] = np.random.randn(10) * 2
    y = X @ true_coef + np.random.randn(n_samples) * 0.5

    # 创建特征名称
    feature_names = [f'feature_{i:02d}' for i in range(X.shape[1])]

    return pd.DataFrame(X, columns=feature_names), pd.Series(y, name='target')


# 生成数据
print("=" * 60)
print("生成包含多重共线性的数据集...")
X, y = create_multicollinear_data(n_samples=1000, n_features=50)

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

print(f"训练集形状: {X_train.shape}")
print(f"测试集形状: {X_test.shape}")
print("=" * 60)

print("\n" + "=" * 60)
print("示例2: 对比三种方法")
print("=" * 60)

# 方法1: 相关系数法
print("\n【方法1: 相关系数法】")
X_train_corr, X_test_corr = remove_multicollinearity(
    X_train.copy(), X_test.copy(),
    # threshold=0.95,
    method='correlation'
)


print(f"\n最终训练集形状: {X_train_corr.shape}")
print(f"最终测试集形状: {X_test_corr.shape}")

# 方法2: 标准VIF法
print("\n【方法2: 标准VIF法】")
X_train_vif, X_test_vif = remove_multicollinearity(
    X_train.copy(), X_test.copy(),
    threshold=5,
    method='vif'
)

# 方法3: 快速VIF法
print("\n【方法3: 快速VIF法】")
X_train_fast, X_test_fast = remove_multicollinearity(
    X_train.copy(), X_test.copy(),
    threshold=5,
    method='fast_vif',
    max_features=50,
    sample_size=500
)

相关系数删除方法

def comprehensive_correlation_removal(X_train, X_test, threshold=0.95, 
                                      target=None, method='variance'):
    """
    完整的相关系数删除方法
    
    Parameters:
    -----------
    method : str
        'variance' - 保留方差大的
        'target' - 保留与目标变量相关性高的
        'hybrid' - 综合两种策略
    """
    print(f"\n  使用相关系数删除共线性特征 (方法: {method})...")
    original_features = X_train.shape[1]
    
    # 计算相关系数矩阵
    corr_matrix = X_train.corr().abs()
    upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))
    
    # 计算特征与目标的相关性（如果提供）
    if target is not None:
        target_corr = X_train.apply(lambda x: x.corr(target)).abs()
    else:
        target_corr = pd.Series(0, index=X_train.columns)
    
    # 计算特征方差
    variances = X_train.var()
    
    # 找出要删除的特征
    to_drop = set()
    processed_pairs = set()
    
    for col in upper.columns:
        # 找出高相关的特征对
        high_corr_cols = [idx for idx in upper.index 
                         if idx != col and upper.loc[idx, col] > threshold 
                         and (idx, col) not in processed_pairs]
        
        for high_corr_col in high_corr_cols:
            processed_pairs.add((col, high_corr_col))
            
            # 根据策略决定保留哪个
            if method == 'variance':
                # 保留方差大的
                if variances[col] >= variances[high_corr_col]:
                    keep, drop = col, high_corr_col
                else:
                    keep, drop = high_corr_col, col
                    
            elif method == 'target':
                # 保留与目标相关性高的
                if target_corr[col] >= target_corr[high_corr_col]:
                    keep, drop = col, high_corr_col
                else:
                    keep, drop = high_corr_col, col
                    
            else:  # hybrid
                # 综合评分
                score_col = variances[col] * target_corr[col]
                score_high = variances[high_corr_col] * target_corr[high_corr_col]
                if score_col >= score_high:
                    keep, drop = col, high_corr_col
                else:
                    keep, drop = high_corr_col, col
            
            to_drop.add(drop)
            print(f"  删除 '{drop}' (保留 '{keep}', r={upper.loc[high_corr_col, col]:.3f})")
    
    # 执行删除
    if to_drop:
        print(f"\n  删除 {len(to_drop)} 个高相关特征")
        X_train = X_train.drop(columns=to_drop)
        X_test = X_test.drop(columns=to_drop)
    
    total_removed = original_features - X_train.shape[1]
    print(f"  ✅ 保留 {X_train.shape[1]} 个特征")
    
    return X_train, X_test