def meanFunction(X):
    return np.sinc(2.*X)  # NumPy的sinc函数自动处理了x=0的情况

def generateData(n):
    """
    生成n个数据点
    n: 要生成的数据点数量
    """
    # 生成均匀分布的X，范围是[-1, 1]
    X = np.random.uniform(-1, 1, size=(n, 1))  # 生成n个样本，每个是1维的
    
    # 计算真实的均值函数值
    y_mean = meanFunction(X)  # 计算每个X对应的真实函数值
    
    # 添加高斯噪声，标准差为0.4
    noise = np.random.normal(0, 0.4, size=y_mean.shape)  # 生成与y_mean相同形状的噪声
    
    # 生成最终的y值：真实函数值 + 噪声
    y = y_mean + noise  # 将噪声加到真实值上
    
    return X, y  # 返回特征X和目标值y

def polyReg(degree):
    model = Pipeline([
        ("poly", PolynomialFeatures(degree=degree)),
        ("lr", LinearRegression(fit_intercept=False))
    ])
    return model

degrees = np.arange(1, 12, 2)  # [1, 3, 5, 7, 9, 11]
fig, ax = plt.subplots(2, 3, figsize=(12, 6), constrained_layout=True)

for id in range(6):
    degree = degrees[id]  # 获取当前的多项式次数
    # 创建并训练模型
    model = polyReg(degree).fit(X, y)  # 用20个样本训练
    
    ii = id // 3
    jj = id % 3
    ax[ii,jj].scatter(X, y, s=10, label='Data')
    ax[ii,jj].plot(xplot, yplot, linestyle='dashed', color='grey', label='Mean function')
    # 在密集的点上预测以得到平滑曲线
    ax[ii,jj].plot(xplot, model.predict(xplot.reshape(-1,1)), 
                   linewidth=1, color='black', label='PolyReg')

# 设置候选的多项式次数
degrees = np.arange(1, 8)  # [1, 2, 3, 4, 5, 6, 7]

# 设置查询点（我们要评估的点）
xquery = np.array([[0.]])  # 选择x=0这个点

# 获取查询点的真实值
ref = meanFunction(xquery)  # 计算sin(0)/0 = 1

# Monte Carlo模拟参数
Nmc = 2000  # 重复2000次实验

# 初始化存储预测结果的数组
pred = np.zeros((Nmc, degrees.shape[0]))  # 2000行7列的矩阵

# 执行Monte Carlo模拟
for imc in range(Nmc):  # 每次循环代表一次新实验
    # 生成新的20个样本的数据集
    X_sim, y_sim = generateData(n)  
    
    # 对每个多项式次数进行拟合和预测
    for ideg, degree in enumerate(degrees):
        # 创建并训练模型
        model = polyReg(degree).fit(X_sim, y_sim)  
        
        # 在x=0处进行预测并存储
        pred[imc, ideg] = model.predict(xquery)[0, 0]  # 获取标量预测值

# 计算偏差和方差
# 偏差 = E[预测值] - 真实值
bias = np.mean(pred, axis=0) - ref.flatten()  

# 方差 = E[(预测值 - E[预测值])^2]
variance = np.var(pred, axis=0)

# 设置多项式的次数范围
degrees = np.arange(1, 8, 1)  # 创建数组[1, 2, 3, 4, 5, 6, 7]

# 创建一个列表来存储所有训练好的模型
models = []  # 初始化空列表

# 训练每个次数的多项式模型
for degree in degrees:
    # 为当前次数创建多项式回归pipeline
    model = polyReg(degree)  # 调用之前定义的polyReg函数
    
    # 使用练习1中的数据(X, y)训练模型
    # X和y是20个样本的训练数据
    model.fit(X, y)  # 拟合模型
    
    # 将训练好的模型添加到列表中
    models.append(model)  # 保存模型供后续使用

# 生成大规模测试集
ntest = 10000  # 测试集大小
Xtest, ytest = generateData(ntest)  # 使用相同的数据生成函数

# 初始化存储误差的列表
erTrain = []  # 训练误差列表
erTest = []   # 测试误差列表

# 计算每个模型的训练误差和测试误差
for i, model in enumerate(models):
    # 步骤1：计算训练误差
    # 在训练数据上进行预测
    y_train_pred = model.predict(X)  # 用训练数据X预测
    # 计算均方误差
    train_mse = mean_squared_error(y, y_train_pred)  # 真实y vs 预测y
    erTrain.append(train_mse)  # 存储训练误差
    
    # 步骤2：计算测试误差（泛化误差）
    # 在测试数据上进行预测
    y_test_pred = model.predict(Xtest)  # 用测试数据Xtest预测
    # 计算均方误差
    test_mse = mean_squared_error(ytest, y_test_pred)  # 真实ytest vs 预测值
    erTest.append(test_mse)  # 存储测试误差

# 设置留一法交叉验证参数
n_splits = n  # 留一法：折数等于样本数(n=20)

# 创建KFold交叉验证对象
kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)
# 参数解释：
# n_splits=20: 将数据分成20份，每次留出1个样本
# shuffle=True: 打乱数据顺序，避免数据顺序带来的偏差
# random_state=42: 固定随机种子，保证结果可重现

# 初始化CV误差存储矩阵
erCV = np.zeros((n, len(degrees)))  # 20行（每个fold）x 7列（每个degree）

# 执行交叉验证循环
for icv, (itrain, itest) in enumerate(kf.split(X)):
    # icv: 当前折的索引（0到19）
    # itrain: 训练集索引数组（19个索引）
    # itest: 测试集索引数组（1个索引，留一法）
    
    # 提取当前折的训练数据和测试数据
    X_train_fold = X[itrain]  # 19个训练样本的特征
    y_train_fold = y[itrain]  # 19个训练样本的目标值
    X_test_fold = X[itest]    # 1个测试样本的特征
    y_test_fold = y[itest]    # 1个测试样本的目标值
    
    # 对每个多项式次数进行训练和测试
    for ideg, degree in enumerate(degrees):
        # 创建新的模型实例（每折都要重新训练）
        model_cv = polyReg(degree)
        
        # 在当前折的训练数据上训练模型
        model_cv.fit(X_train_fold, y_train_fold)  # 用19个样本训练
        
        # 在留出的测试样本上进行预测
        y_pred_fold = model_cv.predict(X_test_fold)  # 预测1个样本
        
        # 计算预测误差（平方误差）
        error = (y_test_fold - y_pred_fold) ** 2  # 计算平方误差
        erCV[icv, ideg] = error[0, 0]  # 存储误差值，注意取标量

# 计算每个degree的平均CV误差
mean_cv_error = np.mean(erCV, axis=0)  # 对20个fold求平均

from sklearn.model_selection import GridSearchCV

# Define the model
model = polyReg(None)

# define the hyperparameters (in this case, the degrees)
param_grid = {
    "poly__degree": degrees  # 双下划线访问Pipeline内部参数
}

polyCV = GridSearchCV(
    estimator=model,
    param_grid=param_grid,
    cv=n,  # 留一法
    scoring="neg_mean_squared_error",  # 负MSE因为要最大化
)

polyCV.fit(X, y)
print('The CV scores are:', -polyCV.cv_results_['mean_test_score'])
print('The best parameter index is:', polyCV.best_index_)
print('The best model is:', polyCV.best_estimator_)

def regularizedPolyReg(degree: int, regType: str, hyper: float):
    """
    构建正则化的多项式回归模型
    
    参数：
    degree: 多项式的次数
    regType: 正则化类型，'ridge'或'lasso'
    hyper: 正则化超参数λ
    """
    
    # 步骤1：根据regType选择回归方法
    if regType.lower() == 'ridge':
        # 创建Ridge回归器，alpha是正则化强度
        regressor = Ridge(alpha=hyper, fit_intercept=False)
    elif regType.lower() == 'lasso':
        # 创建Lasso回归器，alpha是正则化强度
        # max_iter增加到1000000以确保收敛
        regressor = Lasso(alpha=hyper, fit_intercept=False, max_iter=1000000)
    
    # 步骤2：构建Pipeline
    model = Pipeline([
        # 第一步：创建多项式特征
        ("poly", PolynomialFeatures(degree=degree)),
        
        # 第二步：标准化特征（中心化并缩放）
        # 这很重要！确保所有特征在相同的尺度上
        ("scaler", StandardScaler()),
        
        # 第三步：应用正则化回归
        ("regression", regressor)
    ])
    
    return model  # 返回构建好的模型

# 使用练习1中的数据（20个样本）
# 设置10次多项式
degree = 10  # 高次多项式，容易过拟合

# 选择正则化参数
# 注意：由于缩放因子不同，我们为两种方法选择不同的λ
lambda_ridge = 1.0    # Ridge的λ值
lambda_lasso = 0.01   # Lasso的λ值（更小，因为损失函数已经除以n）

# 训练Ridge回归模型
model_ridge = regularizedPolyReg(degree, 'ridge', lambda_ridge)
model_ridge.fit(X, y)  # 用20个样本训练

# 训练Lasso回归模型
model_lasso = regularizedPolyReg(degree, 'lasso', lambda_lasso)
model_lasso.fit(X, y)  # 用20个样本训练

# 训练普通的10次多项式（作为对比）
model_normal = polyReg(degree)  # 使用练习1的函数
model_normal.fit(X, y)

# 从Pipeline中提取系数
coef_ridge = model_ridge.named_steps['regression'].coef_
coef_lasso = model_lasso.named_steps['regression'].coef_

# Ridge：所有系数都非零但被压缩
print(f"Ridge系数的L2范数: {np.linalg.norm(coef_ridge):.4f}")
print(f"Ridge非零系数个数: {np.sum(np.abs(coef_ridge) > 1e-10)}/11")

# Lasso：产生稀疏解
print(f"Lasso系数的L1范数: {np.linalg.norm(coef_lasso, 1):.4f}")
print(f"Lasso非零系数个数: {np.sum(np.abs(coef_lasso) > 1e-10)}/11")

# 显示Lasso的稀疏性
for i, coef in enumerate(coef_lasso[:11]):
    if abs(coef) < 1e-10:
        print(f"  x^{i}: 系数被压缩到0")
    else:
        print(f"  x^{i}: 系数 = {coef:.4f}")

# 生成更大的数据集
np.random.seed(120)  # 设置随机种子
n = 75  # 总样本数
X_full, y_full = generateData(n)  # 生成75个样本

# 步骤1：分割数据集
X_train, X_test, y_train, y_test = train_test_split(
    X_full, y_full,  # 输入数据
    test_size=25,    # 测试集大小（25个样本）
    random_state=42  # 随机种子，保证结果可重现
)

print(f"训练集大小: {X_train.shape[0]} 样本")
print(f"测试集大小: {X_test.shape[0]} 样本")

# 步骤2：为Ridge和Lasso设置超参数候选值
# Ridge的λ范围
ridge_lambdas = np.logspace(-1, 2, 15)  # 从0.1到100，15个值
# Lasso的λ范围
lasso_lambdas = np.logspace(-3, 0, 15)  # 从0.001到1，15个值

# 步骤3：使用GridSearchCV进行超参数选择
from sklearn.model_selection import GridSearchCV

# Ridge的交叉验证
ridge_model = regularizedPolyReg(10, 'ridge', None)  # λ稍后由GridSearch设置
ridge_param_grid = {'regression__alpha': ridge_lambdas}

ridge_cv = GridSearchCV(
    ridge_model,
    ridge_param_grid,
    cv=10,  # 10折交叉验证
    scoring='neg_mean_squared_error',
    n_jobs=-1  # 使用所有CPU核心加速
)
ridge_cv.fit(X_train, y_train)

# Lasso的交叉验证
lasso_model = regularizedPolyReg(10, 'lasso', None)
lasso_param_grid = {'regression__alpha': lasso_lambdas}

lasso_cv = GridSearchCV(
    lasso_model,
    lasso_param_grid,
    cv=10,  # 10折交叉验证
    scoring='neg_mean_squared_error',
    n_jobs=-1
)
lasso_cv.fit(X_train, y_train)

# 步骤4：评估最佳模型在测试集上的表现
# Ridge的最佳模型
best_ridge_lambda = ridge_cv.best_params_['regression__alpha']
ridge_train_score = -ridge_cv.best_score_  # 转换回正的MSE
ridge_test_pred = ridge_cv.predict(X_test)
ridge_test_score = mean_squared_error(y_test, ridge_test_pred)

# Lasso的最佳模型
best_lasso_lambda = lasso_cv.best_params_['regression__alpha']
lasso_train_score = -lasso_cv.best_score_
lasso_test_pred = lasso_cv.predict(X_test)
lasso_test_score = mean_squared_error(y_test, lasso_test_pred)

# 比较分析
if ridge_test_score < lasso_test_score:
    print(f"Ridge表现更好，测试误差低 {(lasso_test_score - ridge_test_score):.4f}")
else:
    print(f"Lasso表现更好，测试误差低 {(ridge_test_score - lasso_test_score):.4f}")