第23章 马氏距离

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.covariance import EmpiricalCovariance
from scipy.stats import chi2
import matplotlib.pyplot as plt

# 步骤1：数据加载与预处理
def preprocess_data(df):
    """预处理信用卡交易数据"""
    # 移除缺失值
    df = df.dropna()
    
    # 分离正常交易（用于训练）
    normal_transactions = df[df['fraud'] == 0].copy()
    X_normal = normal_transactions.drop(['fraud'], axis=1)
    
    # 特征选择：移除高度相关的特征
    corr_matrix = X_normal.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)]
    X_normal = X_normal.drop(to_drop, axis=1)
    
    return X_normal, df

# 步骤2：构建马氏距离模型
class FraudDetector:
    def __init__(self, contamination=0.01, use_pca=True):
        self.contamination = contamination
        self.use_pca = use_pca
        self.pca = None
        self.cov_estimator = None
        self.threshold = None
        
    def fit(self, X_normal):
        """在正常交易上训练模型"""
        n, d = X_normal.shape
        
        # 降维（如果维度过高）
        if self.use_pca and d > 20:
            # 保留95%方差
            self.pca = PCA(n_components=0.95)
            X_train = self.pca.fit_transform(X_normal)
            print(f"PCA降维：{d} → {X_train.shape[1]}")
        else:
            X_train = X_normal
        
        # 估计协方差（根据样本量选择方法）
        n_train, d_train = X_train.shape
        if n_train > 5 * d_train:
            # 样本充足，使用经典估计
            self.cov_estimator = EmpiricalCovariance()
            print("使用经典协方差估计")
        else:
            # 样本不足，使用收缩估计
            from sklearn.covariance import LedoitWolf
            self.cov_estimator = LedoitWolf()
            print("使用Ledoit-Wolf收缩估计")
        
        self.cov_estimator.fit(X_train)
        
        # 计算训练数据的马氏距离
        dist_sq = self.cov_estimator.mahalanobis(X_train)
        
        # 设置阈值（卡方分布）
        self.threshold = chi2.ppf(
            1 - self.contamination, 
            df=d_train
        )
        
        print(f"阈值设置为：{np.sqrt(self.threshold):.2f}")
        print(f"训练集异常率：{(dist_sq > self.threshold).mean():.2%}")
        
        return self
    
    def predict(self, X):
        """预测新交易是否为欺诈"""
        # 应用相同的预处理
        if self.pca is not None:
            X = self.pca.transform(X)
        
        # 计算马氏距离
        dist_sq = self.cov_estimator.mahalanobis(X)
        
        # 判断异常
        is_fraud = dist_sq > self.threshold
        
        return is_fraud, np.sqrt(dist_sq)
    
    def score(self, X):
        """返回异常分数（距离）"""
        if self.pca is not None:
            X = self.pca.transform(X)
        dist_sq = self.cov_estimator.mahalanobis(X)
        return np.sqrt(dist_sq)

# 步骤3：模型评估
def evaluate_model(detector, X_test, y_test):
    """评估模型性能"""
    predictions, distances = detector.predict(X_test)
    
    from sklearn.metrics import (
        classification_report,
        confusion_matrix,
        roc_auc_score,
        roc_curve
    )
    
    # 混淆矩阵
    cm = confusion_matrix(y_test, predictions)
    print("\n混淆矩阵：")
    print(cm)
    
    # 分类报告
    print("\n分类报告：")
    print(classification_report(y_test, predictions))
    
    # ROC-AUC
    auc = roc_auc_score(y_test, distances)
    print(f"\nROC-AUC: {auc:.4f}")
    
    # 可视化
    plt.figure(figsize=(12, 4))
    
    # 距离分布
    plt.subplot(1, 3, 1)
    plt.hist(distances[y_test==0], bins=50, 
            alpha=0.6, label='Normal', density=True)
    plt.hist(distances[y_test==1], bins=30, 
            alpha=0.6, label='Fraud', density=True)
    plt.axvline(np.sqrt(detector.threshold), 
               color='red', linestyle='--', label='Threshold')
    plt.xlabel('Mahalanobis Distance')
    plt.legend()
    plt.title('Distance Distribution')
    
    # ROC曲线
    plt.subplot(1, 3, 2)
    fpr, tpr, _ = roc_curve(y_test, distances)
    plt.plot(fpr, tpr, label=f'AUC={auc:.3f}')
    plt.plot([0, 1], [0, 1], 'k--')
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.legend()
    plt.title('ROC Curve')
    
    # 混淆矩阵热图
    plt.subplot(1, 3, 3)
    plt.imshow(cm, cmap='Blues')
    plt.colorbar()
    plt.title('Confusion Matrix')
    plt.tight_layout()
    plt.show()

# 步骤4：使用流程
# df = pd.read_csv('credit_card_transactions.csv')
# X_normal, df = preprocess_data(df)
# 
# detector = FraudDetector(contamination=0.01, use_pca=True)
# detector.fit(X_normal)
# 
# X_test = df[df['is_test']].drop(['fraud', 'is_test'], axis=1)
# y_test = df[df['is_test']]['fraud']
# evaluate_model(detector, X_test, y_test)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
from matplotlib.patches import Ellipse

def plot_mahalanobis_ellipses(X, levels=[0.95, 0.99]):
    """
    绘制马氏距离等距离椭圆
    
    Parameters:
    -----------
    X : array-like (n, 2)
        二维数据
    levels : list
        置信水平列表
    """
    # 计算均值和协方差
    mean = np.mean(X, axis=0)
    cov = np.cov(X.T)
    
    # 创建图形
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # 子图1：原始数据 + 马氏距离椭圆
    ax1 = axes[0]
    ax1.scatter(X[:, 0], X[:, 1], alpha=0.5, s=20)
    ax1.scatter(mean[0], mean[1], c='red', s=100, 
               marker='x', linewidths=3, label='Mean')
    
    # 绘制马氏距离椭圆
    colors = ['blue', 'green', 'orange']
    for level, color in zip(levels, colors):
        # 卡方分布临界值
        chi2_val = chi2.ppf(level, df=2)
        
        # 计算椭圆参数
        eigenvalues, eigenvectors = np.linalg.eigh(cov)
        angle = np.degrees(np.arctan2(eigenvectors[1, 1], 
                                      eigenvectors[0, 1]))
        width = 2 * np.sqrt(chi2_val * eigenvalues[0])
        height = 2 * np.sqrt(chi2_val * eigenvalues[1])
        
        # 绘制椭圆
        ellipse = Ellipse(mean, width, height, angle=angle,
                         facecolor='none', edgecolor=color, 
                         linewidth=2, label=f'{level*100:.0f}%')
        ax1.add_patch(ellipse)
    
    ax1.set_title('Mahalanobis Distance Ellipses', fontsize=12)
    ax1.set_xlabel('X1')
    ax1.set_ylabel('X2')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.axis('equal')
    
    # 子图2：欧氏距离圆
    ax2 = axes[1]
    ax2.scatter(X[:, 0], X[:, 1], alpha=0.5, s=20)
    ax2.scatter(mean[0], mean[1], c='red', s=100, 
               marker='x', linewidths=3)
    
    # 绘制欧氏距离圆
    for level, color in zip(levels, colors):
        # 计算半径（使用标准差的平均值）
        radius = np.sqrt(chi2.ppf(level, df=2)) * np.mean(np.sqrt(np.diag(cov)))
        circle = plt.Circle(mean, radius, facecolor='none', 
                          edgecolor=color, linewidth=2, 
                          label=f'{level*100:.0f}%')
        ax2.add_patch(circle)
    
    ax2.set_title('Euclidean Distance Circles', fontsize=12)
    ax2.set_xlabel('X1')
    ax2.set_ylabel('X2')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.axis('equal')
    
    # 子图3：马氏距离热图
    ax3 = axes[2]
    
    # 创建网格
    x_min, x_max = X[:, 0].min() - 2, X[:, 0].max() + 2
    y_min, y_max = X[:, 1].min() - 2, X[:, 1].max() + 2
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
                         np.linspace(y_min, y_max, 100))
    
    # 计算每个点的马氏距离
    grid_points = np.c_[xx.ravel(), yy.ravel()]
    diff = grid_points - mean
    inv_cov = np.linalg.inv(cov)
    dist_sq = np.sum(np.dot(diff, inv_cov) * diff, axis=1)
    dist = np.sqrt(dist_sq).reshape(xx.shape)
    
    # 绘制热图
    contour = ax3.contourf(xx, yy, dist, levels=20, cmap='viridis')
    ax3.scatter(X[:, 0], X[:, 1], c='red', s=20, alpha=0.6)
    ax3.scatter(mean[0], mean[1], c='white', s=100, 
               marker='x', linewidths=3)
    
    # 添加等距离线
    contour_lines = ax3.contour(xx, yy, dist, 
                                levels=[np.sqrt(chi2.ppf(0.95, df=2)),
                                       np.sqrt(chi2.ppf(0.99, df=2))],
                                colors=['white', 'yellow'], 
                                linewidths=2, linestyles='--')
    ax3.clabel(contour_lines, inline=True, fontsize=8)
    
    plt.colorbar(contour, ax=ax3, label='Mahalanobis Distance')
    ax3.set_title('Mahalanobis Distance Heatmap', fontsize=12)
    ax3.set_xlabel('X1')
    ax3.set_ylabel('X2')
    
    plt.tight_layout()
    plt.show()

# 示例使用
# 生成相关数据
np.random.seed(42)
mean = [0, 0]
cov = [[2, 1.5], [1.5, 1]]  # 相关性 = 1.5/√(2*1) ≈ 0.75
X = np.random.multivariate_normal(mean, cov, 200)

# 绘制
plot_mahalanobis_ellipses(X, levels=[0.68, 0.95, 0.997])