import numpy as np
from scipy import stats
from scipy.linalg import inv
import pandas as pd

# 合并所有数据
all_provinces = ['山西', '内蒙古', '吉林', '黑龙江', '河南', '甘肃', '青海', '河北', '陕西', '宁夏',
                '新疆', '湖北', '云南', '湖南', '安徽', '贵州', '辽宁', '四川', '山东', '江西',
                '福建', '广西', '海南', '天津', '江苏', '浙江', '北京', '西藏', '上海', '广东']

all_data = [
    # 前20个（第1类G1）
    [8.35, 23.53, 7.51, 8.62, 17.42, 10.00, 1.04, 11.21],
    [9.25, 23.75, 6.61, 9.19, 17.77, 10.48, 1.72, 10.51],
    [8.19, 30.50, 4.72, 9.78, 16.28, 7.60, 2.52, 10.32],
    [7.73, 29.20, 5.42, 9.43, 19.29, 8.49, 2.52, 10.00],
    [9.42, 27.93, 8.20, 8.14, 16.17, 9.42, 1.55, 9.76],
    [9.16, 27.98, 9.01, 9.32, 15.99, 9.10, 1.82, 11.35],
    [10.06, 28.64, 10.52, 10.05, 16.18, 8.39, 1.96, 10.81],
    [9.09, 28.12, 7.40, 9.62, 17.26, 11.12, 2.49, 12.65],
    [9.41, 28.20, 5.77, 10.80, 16.36, 11.56, 1.53, 12.17],
    [8.70, 28.12, 7.21, 10.53, 19.45, 13.30, 1.66, 11.96],
    [6.93, 29.85, 4.54, 9.49, 16.62, 10.65, 1.88, 13.61],
    [8.67, 36.05, 7.31, 7.75, 16.67, 11.68, 2.38, 12.88],
    [9.98, 37.69, 7.01, 8.94, 16.15, 11.08, 0.83, 11.67],
    [6.77, 38.69, 6.01, 8.82, 14.79, 11.44, 1.74, 13.23],
    [8.14, 37.75, 9.61, 8.49, 13.15, 9.76, 1.28, 11.28],
    [7.67, 35.71, 8.04, 8.31, 15.13, 7.76, 1.41, 13.25],
    [7.90, 39.77, 8.49, 12.94, 19.27, 11.05, 2.04, 13.29],
    [7.18, 40.91, 7.32, 8.94, 17.60, 12.75, 1.14, 14.80],
    [8.82, 33.70, 7.59, 10.98, 18.82, 14.73, 1.78, 10.10],
    [6.25, 35.02, 4.72, 6.28, 10.03, 7.15, 1.93, 10.39],
    # 21-27号（第2类G2）
    [10.60, 52.41, 7.70, 9.98, 12.53, 11.70, 2.31, 14.69],  # 福建
    [7.27, 52.65, 3.84, 9.16, 13.03, 15.26, 1.98, 14.57],   # 广西
    [13.45, 55.85, 5.50, 7.45, 9.55, 9.52, 2.21, 16.30],    # 海南
    [10.85, 44.68, 7.32, 14.51, 17.13, 12.08, 1.26, 11.57], # 天津
    [7.21, 45.79, 7.66, 10.36, 16.56, 12.86, 2.25, 11.69],  # 江苏
    [7.68, 50.37, 11.35, 13.30, 19.25, 14.59, 2.75, 14.87], # 浙江
    [7.78, 48.44, 8.00, 20.51, 22.12, 15.73, 1.15, 16.61],  # 北京
    # 最后三个待判别
    [7.94, 39.65, 20.97, 20.82, 22.52, 12.41, 1.75, 7.90],  # 西藏
    [8.28, 64.34, 8.00, 22.22, 20.06, 15.12, 0.72, 22.89],  # 上海
    [12.47, 76.39, 5.52, 11.24, 14.52, 22.00, 5.46, 25.50]  # 广东
]

# 数据分组
X1 = np.array(all_data[:20])     # 第1类G1：前20个省、自治区
X2 = np.array(all_data[20:27])   # 第2类G2：21-27号省、自治区、直辖市
X_test = np.array(all_data[27:]) # 待判别：西藏、上海、广东

print("=" * 60)
print("消费指标数据判别分析")
print("=" * 60)

print(f"\n第1类（G1）样本数量: {len(X1)} (前20个省、自治区)")
print(f"第2类（G2）样本数量: {len(X2)} (21-27号省、自治区、直辖市)")
print(f"待判别样本数量: {len(X_test)} (西藏、上海、广东)")

# 计算各类的均值和协方差矩阵
mean1 = np.mean(X1, axis=0)
mean2 = np.mean(X2, axis=0)

cov1 = np.cov(X1.T)
cov2 = np.cov(X2.T)

print(f"\n第1类均值向量:")
print(mean1)
print(f"\n第2类均值向量:")
print(mean2)

# (1) 进行距离判别、线性及二次判别函数计算，并计算误判率的回代估计与交叉确认估计

print("\n" + "=" * 40)
print("(1) 距离判别、线性及二次判别分析")
print("=" * 40)

# 距离判别函数（马氏距离）
def distance_discriminant(x, mean1, mean2, cov1, cov2):
    """距离判别函数：基于马氏距离"""
    try:
        inv_cov1 = inv(cov1)
        inv_cov2 = inv(cov2)
    except:
        # 如果矩阵奇异，使用伪逆
        inv_cov1 = np.linalg.pinv(cov1)
        inv_cov2 = np.linalg.pinv(cov2)
    
    diff1 = x - mean1
    diff2 = x - mean2
    
    # 计算马氏距离的平方
    d1_squared = diff1 @ inv_cov1 @ diff1
    d2_squared = diff2 @ inv_cov2 @ diff2
    
    return 1 if d1_squared < d2_squared else 2

# 线性判别函数（假设协方差矩阵相等）
pooled_cov = ((len(X1)-1) * cov1 + (len(X2)-1) * cov2) / (len(X1) + len(X2) - 2)

def linear_discriminant(x, mean1, mean2, pooled_cov):
    """线性判别函数"""
    try:
        inv_cov = inv(pooled_cov)
        w = inv_cov @ (mean1 - mean2)
        threshold = 0.5 * (mean1 + mean2) @ inv_cov @ (mean1 - mean2)
        score = w @ x - threshold
        return 1 if score > 0 else 2
    except:
        # 如果矩阵奇异，使用伪逆
        pinv_cov = np.linalg.pinv(pooled_cov)
        w = pinv_cov @ (mean1 - mean2)
        threshold = 0.5 * (mean1 + mean2) @ pinv_cov @ (mean1 - mean2)
        score = w @ x - threshold
        return 1 if score > 0 else 2

# 二次判别函数
def quadratic_discriminant(x, mean1, mean2, cov1, cov2):
    """二次判别函数"""
    try:
        inv_cov1 = inv(cov1)
        inv_cov2 = inv(cov2)
    except:
        # 如果矩阵奇异，使用伪逆
        inv_cov1 = np.linalg.pinv(cov1)
        inv_cov2 = np.linalg.pinv(cov2)
    
    # 计算判别得分
    diff1 = x - mean1
    diff2 = x - mean2
    
    score1 = -0.5 * diff1 @ inv_cov1 @ diff1 - 0.5 * np.log(np.linalg.det(cov1) + 1e-10)
    score2 = -0.5 * diff2 @ inv_cov2 @ diff2 - 0.5 * np.log(np.linalg.det(cov2) + 1e-10)
    
    return 1 if score1 > score2 else 2

# 回代估计（训练样本误判率）
def calculate_resubstitution_error():
    """计算回代估计误判率"""
    # 距离判别
    dist_errors = 0
    for i, x in enumerate(X1):
        if distance_discriminant(x, mean1, mean2, cov1, cov2) != 1:
            dist_errors += 1
    for i, x in enumerate(X2):
        if distance_discriminant(x, mean1, mean2, cov1, cov2) != 2:
            dist_errors += 1
    dist_error_rate = dist_errors / (len(X1) + len(X2))
    
    # 线性判别
    linear_errors = 0
    for i, x in enumerate(X1):
        if linear_discriminant(x, mean1, mean2, pooled_cov) != 1:
            linear_errors += 1
    for i, x in enumerate(X2):
        if linear_discriminant(x, mean1, mean2, pooled_cov) != 2:
            linear_errors += 1
    linear_error_rate = linear_errors / (len(X1) + len(X2))
    
    # 二次判别
    quad_errors = 0
    for i, x in enumerate(X1):
        if quadratic_discriminant(x, mean1, mean2, cov1, cov2) != 1:
            quad_errors += 1
    for i, x in enumerate(X2):
        if quadratic_discriminant(x, mean1, mean2, cov1, cov2) != 2:
            quad_errors += 1
    quad_error_rate = quad_errors / (len(X1) + len(X2))
    
    return dist_error_rate, linear_error_rate, quad_error_rate

# 交叉验证估计（留一法）
def calculate_cross_validation_error():
    """计算交叉验证误判率（留一法）"""
    dist_errors = 0
    linear_errors = 0
    quad_errors = 0
    total_samples = len(X1) + len(X2)
    
    # 对第1类进行交叉验证
    for i in range(len(X1)):
        # 留出第i个样本
        X1_cv = np.delete(X1, i, axis=0)
        X2_cv = X2.copy()
        
        # 重新计算统计量
        mean1_cv = np.mean(X1_cv, axis=0)
        mean2_cv = np.mean(X2_cv, axis=0)
        cov1_cv = np.cov(X1_cv.T)
        cov2_cv = np.cov(X2_cv.T)
        pooled_cov_cv = ((len(X1_cv)-1) * cov1_cv + (len(X2_cv)-1) * cov2_cv) / (len(X1_cv) + len(X2_cv) - 2)
        
        # 预测留出的样本
        x_test = X1[i]
        if distance_discriminant(x_test, mean1_cv, mean2_cv, cov1_cv, cov2_cv) != 1:
            dist_errors += 1
        if linear_discriminant(x_test, mean1_cv, mean2_cv, pooled_cov_cv) != 1:
            linear_errors += 1
        if quadratic_discriminant(x_test, mean1_cv, mean2_cv, cov1_cv, cov2_cv) != 1:
            quad_errors += 1
    
    # 对第2类进行交叉验证
    for i in range(len(X2)):
        # 留出第i个样本
        X1_cv = X1.copy()
        X2_cv = np.delete(X2, i, axis=0)
        
        # 重新计算统计量
        mean1_cv = np.mean(X1_cv, axis=0)
        mean2_cv = np.mean(X2_cv, axis=0)
        cov1_cv = np.cov(X1_cv.T)
        cov2_cv = np.cov(X2_cv.T)
        pooled_cov_cv = ((len(X1_cv)-1) * cov1_cv + (len(X2_cv)-1) * cov2_cv) / (len(X1_cv) + len(X2_cv) - 2)
        
        # 预测留出的样本
        x_test = X2[i]
        if distance_discriminant(x_test, mean1_cv, mean2_cv, cov1_cv, cov2_cv) != 2:
            dist_errors += 1
        if linear_discriminant(x_test, mean1_cv, mean2_cv, pooled_cov_cv) != 2:
            linear_errors += 1
        if quadratic_discriminant(x_test, mean1_cv, mean2_cv, cov1_cv, cov2_cv) != 2:
            quad_errors += 1
    
    return dist_errors/total_samples, linear_errors/total_samples, quad_errors/total_samples

# 计算误判率
print("\n误判率估计:")
resub_dist, resub_linear, resub_quad = calculate_resubstitution_error()
cv_dist, cv_linear, cv_quad = calculate_cross_validation_error()

print(f"回代估计误判率:")
print(f"  距离判别(马氏距离): {resub_dist:.4f} ({resub_dist*100:.2f}%)")
print(f"  线性判别: {resub_linear:.4f} ({resub_linear*100:.2f}%)")
print(f"  二次判别: {resub_quad:.4f} ({resub_quad*100:.2f}%)")

print(f"\n交叉验证误判率:")
print(f"  距离判别(马氏距离): {cv_dist:.4f} ({cv_dist*100:.2f}%)")
print(f"  线性判别: {cv_linear:.4f} ({cv_linear*100:.2f}%)")
print(f"  二次判别: {cv_quad:.4f} ({cv_quad*100:.2f}%)")

# (2) Bayes判别分析
print("\n" + "=" * 40)
print("(2) Bayes判别分析")
print("=" * 40)

# 假设两总体服从正态分布，先验概率按比例分配且误判损失相同
n1, n2 = len(X1), len(X2)
prior1 = n1 / (n1 + n2)  # 先验概率
prior2 = n2 / (n1 + n2)

print(f"先验概率:")
print(f"  P(G1) = {prior1:.4f}")
print(f"  P(G2) = {prior2:.4f}")

def bayes_discriminant_equal_cov(x, mean1, mean2, pooled_cov, prior1, prior2):
    """Bayes判别（协方差矩阵相等）"""
    try:
        inv_cov = inv(pooled_cov)
        # 线性判别函数
        w = inv_cov @ (mean1 - mean2)
        threshold = 0.5 * (mean1 + mean2) @ inv_cov @ (mean1 - mean2) + np.log(prior2/prior1)
        score = w @ x - threshold
        return 1 if score > 0 else 2
    except:
        pinv_cov = np.linalg.pinv(pooled_cov)
        w = pinv_cov @ (mean1 - mean2)
        threshold = 0.5 * (mean1 + mean2) @ pinv_cov @ (mean1 - mean2) + np.log(prior2/prior1)
        score = w @ x - threshold
        return 1 if score > 0 else 2

def bayes_discriminant_unequal_cov(x, mean1, mean2, cov1, cov2, prior1, prior2):
    """Bayes判别（协方差矩阵不等）"""
    try:
        inv_cov1 = inv(cov1)
        inv_cov2 = inv(cov2)
    except:
        inv_cov1 = np.linalg.pinv(cov1)
        inv_cov2 = np.linalg.pinv(cov2)
    
    diff1 = x - mean1
    diff2 = x - mean2
    
    score1 = -0.5 * diff1 @ inv_cov1 @ diff1 - 0.5 * np.log(np.linalg.det(cov1) + 1e-10) + np.log(prior1)
    score2 = -0.5 * diff2 @ inv_cov2 @ diff2 - 0.5 * np.log(np.linalg.det(cov2) + 1e-10) + np.log(prior2)
    
    return 1 if score1 > score2 else 2

# 计算Bayes判别的误判率
def calculate_bayes_error():
    """计算Bayes判别误判率"""
    # 协方差相等情况
    equal_errors = 0
    for x in X1:
        if bayes_discriminant_equal_cov(x, mean1, mean2, pooled_cov, prior1, prior2) != 1:
            equal_errors += 1
    for x in X2:
        if bayes_discriminant_equal_cov(x, mean1, mean2, pooled_cov, prior1, prior2) != 2:
            equal_errors += 1
    equal_error_rate = equal_errors / (len(X1) + len(X2))
    
    # 协方差不等情况
    unequal_errors = 0
    for x in X1:
        if bayes_discriminant_unequal_cov(x, mean1, mean2, cov1, cov2, prior1, prior2) != 1:
            unequal_errors += 1
    for x in X2:
        if bayes_discriminant_unequal_cov(x, mean1, mean2, cov1, cov2, prior1, prior2) != 2:
            unequal_errors += 1
    unequal_error_rate = unequal_errors / (len(X1) + len(X2))
    
    return equal_error_rate, unequal_error_rate

bayes_equal_error, bayes_unequal_error = calculate_bayes_error()

print(f"\nBayes判别误判率:")
print(f"  协方差矩阵相等: {bayes_equal_error:.4f} ({bayes_equal_error*100:.2f}%)")
print(f"  协方差矩阵不等: {bayes_unequal_error:.4f} ({bayes_unequal_error*100:.2f}%)")

# (3) 判别西藏、上海、广东各属哪一类
print("\n" + "=" * 40)
print("(3) 判别西藏、上海、广东的归属")
print("=" * 40)

test_regions = [('西藏', X_test[0]), ('上海', X_test[1]), ('广东', X_test[2])]

for region_name, region_data in test_regions:
    print(f"\n{region_name}的判别结果:")
    
    # 距离判别
    dist_result = distance_discriminant(region_data, mean1, mean2, cov1, cov2)
    print(f"  距离判别(马氏距离): 第{dist_result}类")
    
    # 线性判别
    linear_result = linear_discriminant(region_data, mean1, mean2, pooled_cov)
    print(f"  线性判别: 第{linear_result}类")
    
    # 二次判别
    quad_result = quadratic_discriminant(region_data, mean1, mean2, cov1, cov2)
    print(f"  二次判别: 第{quad_result}类")
    
    # Bayes判别（协方差相等）
    bayes_equal_result = bayes_discriminant_equal_cov(region_data, mean1, mean2, pooled_cov, prior1, prior2)
    print(f"  Bayes判别(协方差相等): 第{bayes_equal_result}类")
    
    # Bayes判别（协方差不等）
    bayes_unequal_result = bayes_discriminant_unequal_cov(region_data, mean1, mean2, cov1, cov2, prior1, prior2)
    print(f"  Bayes判别(协方差不等): 第{bayes_unequal_result}类")

print("\n" + "=" * 60)
print("分析完成")
print("说明:")
print("第1类(G1): 前20个省、自治区")
print("第2类(G2): 21-27号省、自治区、直辖市(福建、广西、海南、天津、江苏、浙江、北京)")
print("待判别: 西藏、上海、广东")
print("=" * 60)