<reponame>Rouzip/Stanford_Machine_Learning<filename>python/ex5/ex5.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.optimize import minimize
from time import sleep


def load_data(filename):
    data = loadmat(filename)
    return data


def plot_data_orginal(X, y):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(X, y, 'rx', 10, 1.5)
    plt.xlabel('Change in water level (x)')
    plt.ylabel('Water flowing out of the dam (y)')
    plt.show()
    plt.close()


def linear_regcost_function(X, y, theta, lambda_, grad_value=False):
    m = y.shape[0]
    y = y.reshape(-1, 1)
    theta = np.reshape(theta, (-1, y.shape[1]))
    h = X@theta
    reg = lambda_ * (theta[1:]**2).sum() / (2 * m)
    J = 1 / (2 * m) * ((h - y)**2).sum() + reg
    if grad_value:
        tmp_theta = theta.copy()
        tmp_theta[0] = 0
        grad = X.T@(h - y) / m + lambda_ * tmp_theta / m
        return J, grad.flatten()
    else:
        return J


def train_linear_reg(X, y, lambda_):
    def cost_func(p):
        # 简化训练参数输入
        return linear_regcost_function(X, y, p, lambda_, True)
    initital_theta = np.zeros(X.shape[1])
    myoptins = {'maxiter': 200, 'disp': False}
    result = minimize(cost_func, x0=initital_theta, options=myoptins,
                      method='L-BFGS-B', jac=True)
    theta = result['x']
    return theta


def learn_curve(X, y, Xval, yval, lambda_):
    '''
    学习曲线图，将训练集和验证集分别的损失计算并返回
    '''
    m = X.shape[0]
    m_val = Xval.shape[0]
    error_train = np.zeros((m, 1))
    error_val = np.zeros((m, 1))
    for i in range(m):
        X_train = X[:i + 1]
        y_train = y[:i + 1]
        theta = train_linear_reg(X_train, y_train, lambda_)
        error_train[i] = linear_regcost_function(
            X_train, y_train, theta, lambda_)
        error_val[i] = linear_regcost_function(Xval, yval, theta, lambda_)
    return error_train, error_val[:m]


def ploy_feature(X, p):
    '''
    将X进行多项式展开
    '''
    m = X.shape[0]
    X_poly = np.matrix(np.zeros((m, p)))

    for i in range(p):
        X_poly[:, i] = X**(i + 1)
    return np.array(X_poly)


def feature_normalize(X):
    mu = np.mean(X, axis=0)
    X_norm = X - mu

    # python默认使用n，而matlab则默认使用n-1
    sigma = np.std(X_norm, axis=0, ddof=1)
    X_norm = X_norm / sigma
    return X_norm, mu, sigma


def plotfit(min_x, max_x, mu, sigma, theta, p, plt):
    X = np.arange(min_x - 15, max_x + 25, 0.05)
    X_poly = ploy_feature(X.reshape(-1, 1), p)
    X_poly = X_poly - mu
    X_poly = X_poly / sigma
    X_poly = np.hstack((np.ones((X.shape[0], 1)), X_poly))
    plt.plot(X, X_poly@theta, '--', linewidth=2)


def validation_curve(X, y, Xval, yval):
    lambda_vec = [0, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10]
    length = len(lambda_vec)
    error_train = np.zeros((length, 1))
    error_val = np.zeros((length, 1))
    for i in range(length):
        lambda_ = lambda_vec[i]
        theta = train_linear_reg(X, y, lambda_)
        error_train[i] = linear_regcost_function(X, y, theta, 0)
        error_val[i] = linear_regcost_function(Xval, yval, theta, 0)
    return lambda_vec, error_train, error_val


if __name__ == '__main__':
    # part0 加载并可视化数据
    data = load_data('./ex5data1.mat')
    X = data['X']
    y = data['y']
    m = X.shape[0]
    train = np.concatenate((np.ones((m, 1)), X), axis=1)
    plot_data_orginal(X, y)
    input('next step')

    # part1 正则化线性回归损失
    theta = np.array([1, 1]).reshape(-1, 1)
    print('预期值为303.993192')
    J = linear_regcost_function(train, y, theta, 1)
    print(J)
    input('next step')

    # part2 正则化线性回归梯度
    print('预期梯度值为-15.303016; 598.250744')
    J, grad = linear_regcost_function(train, y, theta, 1, True)
    print(grad)
    input('next step')

    # part3 训练线性回归
    lambda_ = 0
    theta = train_linear_reg(train, y, lambda_)
    plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
    plt.xlabel('Change in water level (x)')
    plt.ylabel('Water flowing out of the dam (y)')
    plt.plot(X, np.hstack((np.ones((m, 1)), X))@theta, '--', linewidth=2)
    plt.show()
    input('next step')

    # part4 线性回归学习曲线
    Xval = data['Xval']
    Xval_train = np.hstack((np.ones((Xval.shape[0], 1)), Xval))
    yval = data['yval'].flatten()
    error_train, error_val = learn_curve(train, y,
                                         Xval_train, yval, lambda_)
    xaxis = np.array(range(m))
    plt.close()
    plt.plot(xaxis, error_train, xaxis, error_val)
    plt.title('Learning curve for linear regression')
    plt.legend(['Train', 'Cross Validation'])
    plt.xlabel('Number of training examples')
    plt.ylabel('Error')
    plt.axis([0, 13, 0, 150])
    plt.show()

    # part5 多项式回归的特征映射
    p = 8
    X_poly = ploy_feature(X, p)
    X_poly, mu, sigma = feature_normalize(X_poly)
    X_poly = np.hstack((np.ones((m, 1)), X_poly))

    Xtest = data['Xtest']
    ytest = data['ytest']
    X_poly_test = ploy_feature(Xtest, p)
    X_poly_test = X_poly_test - mu
    X_poly_test = X_poly_test / sigma
    X_poly_test = np.hstack((np.ones((X_poly_test.shape[0], 1)), X_poly_test))

    X_poly_val = ploy_feature(Xval, p)
    X_poly_val = X_poly_val - mu
    X_poly_val = X_poly_val / sigma
    X_poly_val = np.hstack((np.ones((X_poly_val.shape[0], 1)), X_poly_val))
    print('归一化后第一行为：{}'.format(X_poly[0, :]))
    input('next step')

    # part6 多项式回归的学习曲线
    lambda_ = 0
    theta = train_linear_reg(X_poly, y, lambda_)
    print('theta is: {}'.format(theta))
    fig1 = plt.figure(1)
    plt.close()
    plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
    plotfit(min(X), max(X), mu, sigma, theta, p, plt)
    plt.show()

    fig2 = plt.figure(2)
    error_train, error_val = learn_curve(X_poly, y, X_poly_val, yval, lambda_)
    plt.close()
    plt.plot(xaxis, error_train, xaxis, error_val)
    plt.title('Polynomial Regression Learning Curve (lambda = {})'.format(lambda_))
    plt.xlabel('Number of training examples')
    plt.ylabel('Error')
    plt.legend(['Train', 'Cross Validation'])
    plt.axis([0, 13, 0, 100])
    plt.show()

    print('Polynomial Regression (lambda = {})\n\n'.format(lambda_))
    print('#Training Examples\tTrain Error\tCross Validation Error\n')
    for i in range(m):
        print('\t {}\t\t{}\t{}\n'.format(i + 1, error_train[i], error_val[i]))

    input('next step')

    # part7 对于lambda进行挑选
    lambda_vec, error_train, error_val = validation_curve(
        X_poly, y, X_poly_val, yval)
    plt.close('all')
    p1, p2 = plt.plot(lambda_vec, error_train, lambda_vec, error_val)
    plt.legend(['Train', 'Cross Validation'])
    plt.xlabel('lambda')
    plt.ylabel('Error')
    plt.show()

    print('lambda\t\tTrain Error\tValidation Error\n')
    for i in range(len(lambda_vec)):
        print(' {}\t{}\t{}\n'.format(
            lambda_vec[i], error_train[i], error_val[i]))
    input('next step')

    # part8 将最优lambda在测试集中计算损失
    lambda_val = 3
    theta = train_linear_reg(X_poly, y, lambda_val)
    error_test = linear_regcost_function(X_poly_test, ytest, theta, 0)
    # 预期值为3.859
    print('Test set error :{}\n'.format(error_test))
    input('next step')

    # part9 随机挑选几个例子进行学习曲线可视化
    lambda_val = 0.01
    times = 50

    error_train_rand = np.zeros((m, times))
    error_val_rand = np.zeros((m, times))

    for i in range(1, m + 1):
        for j in range(times):
            # 随机从X_ploy和X_val中选取一些索引训练theta
            rand_sample = np.random.permutation(m)
            rand_train_index = rand_sample[:i]
            X_train_rand = X_poly[rand_train_index, :]
            y_train_rand = y[rand_train_index]

            rand_sample = np.random.permutation(X_poly_test.shape[0])
            rand_val_index = rand_sample[:i]
            X_val_rand = X_poly_val[rand_val_index, :]
            y_val_rand = yval[rand_val_index]

            theta = train_linear_reg(X_train_rand, y_train_rand, lambda_val)
            error_train_rand[i - 1, j] = \
                linear_regcost_function(X_train_rand, y_train_rand, theta, 0)
            error_val_rand[i - 1, j] = \
                linear_regcost_function(X_val_rand, y_val_rand, theta, 0)
    error_train = np.mean(error_train_rand, axis=1)
    error_val = np.mean(error_val_rand, axis=1)

    # 可视化曲线
    plt.close('all')
    p1, p2 = plt.plot(xaxis, error_train, xaxis, error_val)
    plt.title(
        'Polynomial Regression Learning Curve (lambda = {})'.format(lambda_val))
    plt.xlabel('Number of training examples')
    plt.ylabel('Error')
    plt.legend((p1, p2), ('Train', 'Cross Validation'))
    plt.axis([0, 13, 0, 150])
    plt.show()

    print('# Training Examples\tTrain Error\tCross Validation Error\n')
    for i in range(m):
        print('  \t{:d}\t\t{:f}\t{:f}\n'.format(
            i + 1, error_train[i], error_val[i]))

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.optimize import minimize
from time import sleep


def load_data(filename):
    data = loadmat(filename)
    return data


def plot_data_orginal(X, y):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(X, y, 'rx', 10, 1.5)
    plt.xlabel('Change in water level (x)')
    plt.ylabel('Water flowing out of the dam (y)')
    plt.show()
    plt.close()


def linear_regcost_function(X, y, theta, lambda_, grad_value=False):
    m = y.shape[0]
    y = y.reshape(-1, 1)
    theta = np.reshape(theta, (-1, y.shape[1]))
    h = X@theta
    reg = lambda_ * (theta[1:]**2).sum() / (2 * m)
    J = 1 / (2 * m) * ((h - y)**2).sum() + reg
    if grad_value:
        tmp_theta = theta.copy()
        tmp_theta[0] = 0
        grad = X.T@(h - y) / m + lambda_ * tmp_theta / m
        return J, grad.flatten()
    else:
        return J


def train_linear_reg(X, y, lambda_):
    def cost_func(p):
        # 简化训练参数输入
        return linear_regcost_function(X, y, p, lambda_, True)
    initital_theta = np.zeros(X.shape[1])
    myoptins = {'maxiter': 200, 'disp': False}
    result = minimize(cost_func, x0=initital_theta, options=myoptins,
                      method='L-BFGS-B', jac=True)
    theta = result['x']
    return theta


def learn_curve(X, y, Xval, yval, lambda_):
    '''
    学习曲线图，将训练集和验证集分别的损失计算并返回
    '''
    m = X.shape[0]
    m_val = Xval.shape[0]
    error_train = np.zeros((m, 1))
    error_val = np.zeros((m, 1))
    for i in range(m):
        X_train = X[:i + 1]
        y_train = y[:i + 1]
        theta = train_linear_reg(X_train, y_train, lambda_)
        error_train[i] = linear_regcost_function(
            X_train, y_train, theta, lambda_)
        error_val[i] = linear_regcost_function(Xval, yval, theta, lambda_)
    return error_train, error_val[:m]


def ploy_feature(X, p):
    '''
    将X进行多项式展开
    '''
    m = X.shape[0]
    X_poly = np.matrix(np.zeros((m, p)))

    for i in range(p):
        X_poly[:, i] = X**(i + 1)
    return np.array(X_poly)


def feature_normalize(X):
    mu = np.mean(X, axis=0)
    X_norm = X - mu

    # python默认使用n，而matlab则默认使用n-1
    sigma = np.std(X_norm, axis=0, ddof=1)
    X_norm = X_norm / sigma
    return X_norm, mu, sigma


def plotfit(min_x, max_x, mu, sigma, theta, p, plt):
    X = np.arange(min_x - 15, max_x + 25, 0.05)
    X_poly = ploy_feature(X.reshape(-1, 1), p)
    X_poly = X_poly - mu
    X_poly = X_poly / sigma
    X_poly = np.hstack((np.ones((X.shape[0], 1)), X_poly))
    plt.plot(X, X_poly@theta, '--', linewidth=2)


def validation_curve(X, y, Xval, yval):
    lambda_vec = [0, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10]
    length = len(lambda_vec)
    error_train = np.zeros((length, 1))
    error_val = np.zeros((length, 1))
    for i in range(length):
        lambda_ = lambda_vec[i]
        theta = train_linear_reg(X, y, lambda_)
        error_train[i] = linear_regcost_function(X, y, theta, 0)
        error_val[i] = linear_regcost_function(Xval, yval, theta, 0)
    return lambda_vec, error_train, error_val


if __name__ == '__main__':
    # part0 加载并可视化数据
    data = load_data('./ex5data1.mat')
    X = data['X']
    y = data['y']
    m = X.shape[0]
    train = np.concatenate((np.ones((m, 1)), X), axis=1)
    plot_data_orginal(X, y)
    input('next step')

    # part1 正则化线性回归损失
    theta = np.array([1, 1]).reshape(-1, 1)
    print('预期值为303.993192')
    J = linear_regcost_function(train, y, theta, 1)
    print(J)
    input('next step')

    # part2 正则化线性回归梯度
    print('预期梯度值为-15.303016; 598.250744')
    J, grad = linear_regcost_function(train, y, theta, 1, True)
    print(grad)
    input('next step')

    # part3 训练线性回归
    lambda_ = 0
    theta = train_linear_reg(train, y, lambda_)
    plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
    plt.xlabel('Change in water level (x)')
    plt.ylabel('Water flowing out of the dam (y)')
    plt.plot(X, np.hstack((np.ones((m, 1)), X))@theta, '--', linewidth=2)
    plt.show()
    input('next step')

    # part4 线性回归学习曲线
    Xval = data['Xval']
    Xval_train = np.hstack((np.ones((Xval.shape[0], 1)), Xval))
    yval = data['yval'].flatten()
    error_train, error_val = learn_curve(train, y,
                                         Xval_train, yval, lambda_)
    xaxis = np.array(range(m))
    plt.close()
    plt.plot(xaxis, error_train, xaxis, error_val)
    plt.title('Learning curve for linear regression')
    plt.legend(['Train', 'Cross Validation'])
    plt.xlabel('Number of training examples')
    plt.ylabel('Error')
    plt.axis([0, 13, 0, 150])
    plt.show()

    # part5 多项式回归的特征映射
    p = 8
    X_poly = ploy_feature(X, p)
    X_poly, mu, sigma = feature_normalize(X_poly)
    X_poly = np.hstack((np.ones((m, 1)), X_poly))

    Xtest = data['Xtest']
    ytest = data['ytest']
    X_poly_test = ploy_feature(Xtest, p)
    X_poly_test = X_poly_test - mu
    X_poly_test = X_poly_test / sigma
    X_poly_test = np.hstack((np.ones((X_poly_test.shape[0], 1)), X_poly_test))

    X_poly_val = ploy_feature(Xval, p)
    X_poly_val = X_poly_val - mu
    X_poly_val = X_poly_val / sigma
    X_poly_val = np.hstack((np.ones((X_poly_val.shape[0], 1)), X_poly_val))
    print('归一化后第一行为：{}'.format(X_poly[0, :]))
    input('next step')

    # part6 多项式回归的学习曲线
    lambda_ = 0
    theta = train_linear_reg(X_poly, y, lambda_)
    print('theta is: {}'.format(theta))
    fig1 = plt.figure(1)
    plt.close()
    plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
    plotfit(min(X), max(X), mu, sigma, theta, p, plt)
    plt.show()

    fig2 = plt.figure(2)
    error_train, error_val = learn_curve(X_poly, y, X_poly_val, yval, lambda_)
    plt.close()
    plt.plot(xaxis, error_train, xaxis, error_val)
    plt.title('Polynomial Regression Learning Curve (lambda = {})'.format(lambda_))
    plt.xlabel('Number of training examples')
    plt.ylabel('Error')
    plt.legend(['Train', 'Cross Validation'])
    plt.axis([0, 13, 0, 100])
    plt.show()

    print('Polynomial Regression (lambda = {})\n\n'.format(lambda_))
    print('#Training Examples\tTrain Error\tCross Validation Error\n')
    for i in range(m):
        print('\t {}\t\t{}\t{}\n'.format(i + 1, error_train[i], error_val[i]))

    input('next step')

    # part7 对于lambda进行挑选
    lambda_vec, error_train, error_val = validation_curve(
        X_poly, y, X_poly_val, yval)
    plt.close('all')
    p1, p2 = plt.plot(lambda_vec, error_train, lambda_vec, error_val)
    plt.legend(['Train', 'Cross Validation'])
    plt.xlabel('lambda')
    plt.ylabel('Error')
    plt.show()

    print('lambda\t\tTrain Error\tValidation Error\n')
    for i in range(len(lambda_vec)):
        print(' {}\t{}\t{}\n'.format(
            lambda_vec[i], error_train[i], error_val[i]))
    input('next step')

    # part8 将最优lambda在测试集中计算损失
    lambda_val = 3
    theta = train_linear_reg(X_poly, y, lambda_val)
    error_test = linear_regcost_function(X_poly_test, ytest, theta, 0)
    # 预期值为3.859
    print('Test set error :{}\n'.format(error_test))
    input('next step')

    # part9 随机挑选几个例子进行学习曲线可视化
    lambda_val = 0.01
    times = 50

    error_train_rand = np.zeros((m, times))
    error_val_rand = np.zeros((m, times))

    for i in range(1, m + 1):
        for j in range(times):
            # 随机从X_ploy和X_val中选取一些索引训练theta
            rand_sample = np.random.permutation(m)
            rand_train_index = rand_sample[:i]
            X_train_rand = X_poly[rand_train_index, :]
            y_train_rand = y[rand_train_index]

            rand_sample = np.random.permutation(X_poly_test.shape[0])
            rand_val_index = rand_sample[:i]
            X_val_rand = X_poly_val[rand_val_index, :]
            y_val_rand = yval[rand_val_index]

            theta = train_linear_reg(X_train_rand, y_train_rand, lambda_val)
            error_train_rand[i - 1, j] = \
                linear_regcost_function(X_train_rand, y_train_rand, theta, 0)
            error_val_rand[i - 1, j] = \
                linear_regcost_function(X_val_rand, y_val_rand, theta, 0)
    error_train = np.mean(error_train_rand, axis=1)
    error_val = np.mean(error_val_rand, axis=1)

    # 可视化曲线
    plt.close('all')
    p1, p2 = plt.plot(xaxis, error_train, xaxis, error_val)
    plt.title(
        'Polynomial Regression Learning Curve (lambda = {})'.format(lambda_val))
    plt.xlabel('Number of training examples')
    plt.ylabel('Error')
    plt.legend((p1, p2), ('Train', 'Cross Validation'))
    plt.axis([0, 13, 0, 150])
    plt.show()

    print('# Training Examples\tTrain Error\tCross Validation Error\n')
    for i in range(m):
        print('  \t{:d}\t\t{:f}\t{:f}\n'.format(
            i + 1, error_train[i], error_val[i]))
# 简化训练参数输入 学习曲线图，将训练集和验证集分别的损失计算并返回 将X进行多项式展开 # python默认使用n，而matlab则默认使用n-1 # part0 加载并可视化数据 # part1 正则化线性回归损失 # part2 正则化线性回归梯度 # part3 训练线性回归 # part4 线性回归学习曲线 # part5 多项式回归的特征映射 # part6 多项式回归的学习曲线 # part7 对于lambda进行挑选 # part8 将最优lambda在测试集中计算损失 # 预期值为3.859 # part9 随机挑选几个例子进行学习曲线可视化 # 随机从X_ploy和X_val中选取一些索引训练theta # 可视化曲线