'''
在本练习中，您将实现正则化的线性回归和多项式回归，并使用它来研究具有不同偏差-方差属性的模型.
在前半部分的练习中，你将实现正则化线性回归，以预测水库中的水位变化，从而预测大坝流出的水量。
在下半部分中，您将通过一些调试学习算法的诊断，并检查偏差 v.s. 方差的影响。
我们将从可视化数据集开始，其中包含水位变化的历史记录，x，以及从大坝流出的水量，y。
这个数据集分为了三个部分：
training set 训练集：训练模型
cross validation set 交叉验证集：选择正则化参数
test set 测试集：评估性能，模型训练中不曾用过的样本
'''
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
import scipy.optimize as opt

# 读取数据
path = 'data/ex5data1.mat'  # 数据路径
data = loadmat(path)  # 加载数据
X, y = data['X'], data['y']  # 获取 训练集
Xval, yval = data['Xval'], data['yval']  # 获取 交叉验证集
Xtest, ytest = data['Xtest'], data['ytest']  # 训练集
X = np.insert(X, 0, 1, axis=1)  # X的第一列插入1
Xval = np.insert(Xval, 0, 1, axis=1)  # Xval的第一列插入1
Xtest = np.insert(Xtest, 0, 1, axis=1)  # Xtest的第一列插入1
# 打印X，Xval，Xtest的大小
# print('X={},y={}'.format(X.shape, y.shape))
# print('Xval={},yval={}'.format(Xval.shape, yval.shape))
# print('Xtest={},ytest={}'.format(Xtest.shape, ytest.shape))
'''
结果：
X=(12, 2),y=(12, 1)
Xval=(21, 2),yval=(21, 1)
Xtest=(21, 2),ytest=(21, 1)
'''


def plotData():  # 将训练集的数据可视化为散点图
    plt.figure(figsize=(8, 5))  # 图大小
    plt.scatter(X[:, 1:], y, c='r', marker='x')  # X轴，y轴，颜色，点样式
    plt.xlabel('Change in water level (x)')
    plt.ylabel('Water flowing out of the dam (y)')
    plt.grid(True)  # 生成网格


# plotData() #运行画图函数
# plt.show()

# 计算带正则化项的线性代价值
def costReg(theta, X, y, l):
    """
        不正则化 theta0
        theta:(n+1,)
        X:(m, n+1)
        y:(m, 1)
    """
    cost = ((X @ theta - y.flatten()) ** 2).sum()
    regterm = l * (theta[1:] @ theta[1:])  # 正则化项
    return (cost + regterm) / (2 * len(X))  # 返回带正则化项的线性代价值


# 使用theta=[1，1]，lambda = 1运行代价计算
# theta = np.ones(X.shape[1])  # [1,1]
# print(costReg(theta, X, y, l=1))  # 结果：303.9931922202643


# 计算带正则化项的线性梯度值
def gradientReg(theta, X, y, l):
    """
        不正则化 theta0
        theta:(n+1,) theta:(2,)
        X:(m, n+1) X: (12, 2)
        y:(m, 1)  y: (12, 1)
        返回结果：(2,)
    """
    grad = (X @ theta - y.flatten()) @ X
    regterm = l * theta  # 正则化项
    regterm[0] = 0
    return (grad + regterm) / len(X)  # 返回带正则化项的线性梯度值


# 使用theta=[1，1]，lambda = 1运行梯度计算
# theta = np.ones(X.shape[1])
# print(gradientReg(theta, X, y, 1))  # 结果：[-15.30301567 598.25074417]

# 拟合线性回归
def trainLinearReg(X, y, l):
    theta = np.ones(X.shape[1])  # 初始化theta为[1，1]
    res = opt.minimize(fun=costReg, x0=theta, args=(X, y, l), method='TNC', jac=gradientReg)
    return res.x  # 结果返回优化后的theta:(2,)


fit_theta = trainLinearReg(X, y, 0)
print('优化后的theta:', fit_theta)  # 结果：[13.08790367  0.36777923]
plotData()
plt.plot(X[:, 1:], X @ fit_theta)  # 画线性拟合直线

'''
这里我们把λ = 0，因为线性回归只有两个参数，这么低的维度，正则化并没有用。
从图中可以看到，拟合最好的这条直线告诉我们这个模型并不适合这个数据。
接下来，您将实现一个函数来生成学习曲线，它可以帮助您调试学习算法。
'''
'''
机器学习中一个重要的概念是偏差（bias）和方差（variance）的权衡。
高偏差意味着欠拟合，高方差意味着过拟合。
在这部分练习中，您将在学习曲线上绘制训练误差和验证误差，以诊断bias-variance问题。
'''
'''
训练样本数量从1开始逐渐增加，训练出不同的参数向量θ。
接着通过交叉验证样本Xval计算验证误差。
切记此时不要使用正则化，将λ=0。
计算交叉验证代价时记得整个交叉验证集来计算。
'''


# 画出学习曲线，即交叉验证误差和训练误差随样本数量的变化的变化
def plot_learning_curve(X, y, Xval, yval, l):
    xx = range(1, len(X) + 1)  # [1,2,3,4,...,12]
    training_cost, cv_cost = [], []
    for i in xx:
        res = trainLinearReg(X[:i], y[:i], l)  # 样本数量逐层增加,res存储每次的theta值
        training_cost_i = costReg(res, X[:i], y[:i], 0)  # 训练集代价
        cv_cost_i = costReg(res, Xval, yval, 0)  # 交叉验证集代价
        training_cost.append(training_cost_i)
        cv_cost.append(cv_cost_i)
    # 画学习曲线图
    plt.figure(figsize=(8, 5))
    plt.plot(xx, training_cost, label='training cost')
    plt.plot(xx, cv_cost, label='cv cost')
    plt.legend()
    plt.xlabel('Number of training examples')
    plt.ylabel('Cost')
    plt.title('Learning curve for linear regression')
    plt.grid(True)


plot_learning_curve(X, y, Xval, yval, l=0)
# 从图中看出来，随着样本数量的增加，训练误差和交叉验证误差都很高，这属于高偏差，欠拟合。


'''
我们的线性模型对于数据来说太简单了，导致了欠拟合(高偏差)。
接下来，您将通过添加更多的特性来解决这个问题。

数据预处理:
X，Xval，Xtest都需要添加多项式特征，这里我们选择增加到6次方(不要忘了标准化)。
'''


def genPolyFeatures(X, power):
    """ 添加多项式特征:
        每次在最后一列插入第二列的 2,3,...,power 次方（第一列为偏置）
        从二次方开始插入（因为本身含有一列一次方）
    """
    Xpoly = X.copy()  # 浅拷贝
    for i in range(2, power + 1):
        Xpoly = np.insert(Xpoly, Xpoly.shape[1], np.power(Xpoly[:, 1], i), axis=1)  # 按列插入
    return Xpoly  # 返回添加多项式特征后的数据集


def get_means_std(X):
    """
    获取训练集的均值和误差，用来标准化所有数据。
    """
    means = np.mean(X, axis=0)  # 对各列求均值，返回 1* n 矩阵
    """
    numpy.mean(a, axis, dtype, out，keepdims )
    函数功能：求取均值
    经常操作的参数为axis，以m * n矩阵举例：
    axis 不设置值，对 m*n 个数求均值，返回一个实数
    axis = 0：压缩行，对各列求均值，返回 1* n 矩阵
    axis =1 ：压缩列，对各行求均值，返回 m *1 矩阵
    """
    stds = np.std(X, axis=0, ddof=1)  # 计算每一列的无偏样本标准差
    """
    numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=<no value>)
    函数功能：求标准差
    axis=0时，表示求每一列标准差
    axis=1时，表示求每一行标准差
    axis=None时，表示求全局标准差。
    
    ddof=0时，计算有偏样本标准差
    (一般在拥有所有数据的情况下，计算所有数据的标准差时使用，即最终除以n)
    ddof=1时，表示计算无偏样本标准差
    (最终除以n-1)
    """
    return means, stds


def featureNormalize(myX, means, stds):
    """标准化
    所有数据集应该都用训练集的均值和样本标准差处理。
    """
    X_norm = myX.copy()
    X_norm[:, 1:] = X_norm[:, 1:] - means[1:]
    X_norm[:, 1:] = X_norm[:, 1:] / stds[1:]
    return X_norm


# 获取添加多项式特征以及标准化之后的数据集。
power = 6  # 扩展到x的6次方
train_means, train_stds = get_means_std(genPolyFeatures(X, power))  # 获得训练集的均值和标准差
X_norm = featureNormalize(genPolyFeatures(X, power), train_means, train_stds)  # 扩展特征后的训练集
Xval_norm = featureNormalize(genPolyFeatures(Xval, power), train_means, train_stds)  # 扩展特征后的交叉验证集
Xtest_norm = featureNormalize(genPolyFeatures(Xtest, power), train_means, train_stds)  # 扩展特征后的测试集


# 画出多项式回归的拟合曲线
def plot_fit(means, stds, l):
    theta = trainLinearReg(X_norm, y, l)  # 得到扩大特征后的训练集的最优theta值 (7, )
    print('此时的theta:', theta)
    x = np.linspace(-75, 55, 50)  # 50表示要生成的样本数
    xmat = x.reshape(-1, 1)  # 重塑成一列 (50,1)
    xmat = np.insert(xmat, 0, 1, axis=1)  # 在第一列插入1 (50,2)
    Xmat = genPolyFeatures(xmat, power)  # 扩大特征 (50,7)
    Xmat_norm = featureNormalize(Xmat, means, stds)  # 标准化

    plotData()  # 画原数据散点图
    plt.plot(x, Xmat_norm @ theta, 'b--')  # 画拟合曲线


plot_fit(train_means, train_stds, 0)
# 此时的theta: [ 11.21817067  11.36845531  13.44911626  10.73084633  -4.44505759  -11.89800257  -5.06836084]
plot_learning_curve(X_norm, y, Xval_norm, yval, 0)

plot_fit(train_means, train_stds, 100)
# 此时的theta: [11.21758894  0.98876228  0.30025154  0.77917857  0.11521526  0.5903675  -0.0118514 ]
plot_learning_curve(X_norm, y, Xval_norm, yval, 100)
'''
l=0时，不进行正则化，图像显示过拟合特征。
l=100时，惩罚过多，图像显示欠拟合特征。
'''

# 使用交叉验证集来选择最合适的正则化lambda值
lambdas = [0., 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1., 3., 10.]  # lambda候选值
error_train, error_val = [], []
for l in lambdas:
    theta = trainLinearReg(X_norm, y, l)
    error_train.append(costReg(theta, X_norm, y, 0))  # 这里计算error时，令l=0
    error_val.append(costReg(theta, Xval_norm, yval, 0))
    '''
    正则化是用来进行梯度下降的，惩罚某些特征。
    当梯度下降完成，用得到的最优参数计算损失函数，就不用正则化那一项了。
    因此这里计算error时，令l=0。
    '''

l = lambdas[np.argmin(error_val)]  # np.argmin(error_val)：获得error_val最小值处的索引号
print('最优lambda值为:', l)
theta = trainLinearReg(X_norm, y, l)
print('此时theta值为：', theta)
print('test cost=', costReg(theta, Xtest_norm, ytest, 0))
'''
运行结果：
最优lambda值为: 3.0
此时theta值为： [11.21758982  6.77857869  3.98943697  3.99183811  2.24949203  2.45682122  1.25070152]
test cost= 4.755261169964735

若power=8，则结果为：
最优lambda值为: 3.0
此时theta值为： [11.2175881   6.65682303  3.89421607  3.75899131  2.25315894  2.20639341  1.33198932  1.35014165  0.76188972]
test cost= 3.859893821669393
'''

# 画lambda-cost图
plt.figure(figsize=(8, 5))
plt.plot(lambdas, error_train, label='Train')
plt.plot(lambdas, error_val, label='Cross Validation')
plt.legend()
plt.xlabel('lambda')
plt.ylabel('Error')
plt.grid(True)
plt.show()

运行结果：

线性回归及其学习曲线：



多项式回归及其学习曲线：

过拟合：


欠拟合：


lambda的选取：

正则化代价函数：

正则化线性梯度值:

训练集代价函数:

多项式回归:

linspace函数:

reshape函数:
在创建DataFrame的时候常常使用reshape来更改数据的列数和行数。
在使用了reshape（-1，1）之后，数据集变成了一列。reshape(1,-1)变成了一行。