多项式拟合-线性回归 ·

Machine Learning Exercise 1

实验题目

编写程序:模拟仿真多项式回归参见 textbook p4-12(PRML)。完成以下任务:

1) 生成正弦序列 s(n);
2) 使用噪声函数对正弦序列加噪 x(n)=s(n)+w(n);
3) 使用多项式回归模型对 x(n)进行拟合，并分析过拟合和欠拟合情况

注:参考误差函数式 1-2，带正则项的修正误差函数式 1-4，实验仿真生成图 1- 6、图 1-7，给出模型系数表。

实验过程及代码

实验过程按照题目要求可以拆解为4步完成。首先进行数据采样，添加噪声；接着定义损失函数，以及正则项；
之后采用梯度下降寻找局部最优解；再通过画图可视化实验结果。以下为具体描述：

数据生成

利用numpy在sin(2*pi*x)的基础上生成点，并添加高斯噪声（其中均值设为0，方差设为0.3）。
为了方便后序的实验，我们在生成数据x, y的同时添加阶数，每个x调整为形如
x = [1, x, x^2, ..., x^M]的形式。最前面的1对应目标函数的b。

数据生成代码如下：

# define gaussion noise function
def noise(mu=0, sigma=0.3):
    '''
    dtype:
        mu: default 0
        sigma: default 0.3
    rtype:
        noise: float
    '''
    return np.random.normal(mu, sigma)


def genTargetWithNoise(x, noise=noise):
    ''' target = sin(2*pi*x) + noise
    dtype:
        x: float
        noise: noise function  
    rtype:
        t: float
    '''
#     assert 0 <= x <= 1
    t = np.sin(2*np.pi*x)
    if noise:
        t += noise()
    return t

def regenX(random, poly_degree):
    xs = [[x ** i for i in range(1, poly_degree+1)] for x in random]
    for x in xs:
        x.insert(0, 1.)
    return xs


def genData(numPoints, noise, poly_degree):
    '''Generate data.
    '''
    # random sample floats in the half-open interval [0.0, 1.0).
    random = np.random.random_sample((numPoints,))

    # generate targets for random
    targets = list(map(genTargetWithNoise, random))

    # remake data to ractangle form and add 1. to all x to correspond with b
    xs = regenX(random, poly_degree)

    return np.array(xs), np.array(targets)

损失函数

这里我们按照《PRML》中式1-2, 1-3的定义，使用平方误差与均方误差。同时添加正则项。
如式1-4。

代码如下：

# define a polynomial
def f(xs, theta):
    xs = np.asarray(xs)
#     print(xs.shape)
    if xs.shape[0] > 1:
        return [x.T.dot(theta) for x in xs]
    return xs.T.dot(theta)

# define squares error
def LSE_loss(y_true, y_hat, theta, penalization=0.1):
    y_true = np.asarray(y_true)
    y_hat = np.asarray(y_hat)

    # set the regularizer  
    regularizer = (penalization / 2) * (np.dot(theta.T, theta))

    return (1/2) * np.square(y_hat - y_true).sum() + regularizer

# define mean squares error loss
def MSE_loss(y_true, y_hat, theta, penalization=0.1):
    y_true = np.asarray(y_true)
    y_hat = np.asarray(y_hat)

    # set the regularizer  
    regularizer = (penalization / 2) * (np.dot(theta.T, theta))

    # compute mse
    mse = np.sqrt(np.square(y_hat - y_true).sum()/len(y_hat)) + regularizer

    return mse

梯度下降

上面两步我们分别准备好了数据以及损失函数，这一步我们设置一些常用的参数，利用梯度下降法寻找
局部最优解。定义的多项式函数如下:

这里我们的theta0设置为1，由于我们在准备数据的时候已经处理好了输入x，因此我们的bias项可以直接并入到w中，构成新的theta。

对参数theta求导:

更新theta:

代码如下：

# gd
def gradientDescent(x, y, x_val, y_val, theta, lr, sample_num, numIterations, loss_name='LSE', penalization=0):
    xTrans = x.transpose()
    record = []
    freq = 100 if numIterations < 10000 else 1000
    for i in range(0, numIterations+1):
        y_hat = f(x, theta)

        if loss_name == 'LSE':
            loss = LSE_loss(y, y_hat, theta, penalization)
        elif loss_name == 'MSE':
            loss = MSE_loss(y, y_hat, theta, penalization)

        if i % freq == 0:
            print("Iteration %d | %s loss: %f" % (i, loss_name,  loss))
            y_val_hat = f(x_val, theta)
            if loss_name == 'LSE':
                val_loss = LSE_loss(y_val, y_val_hat, theta)
            elif loss_name == 'MSE':
                val_loss = MSE_loss(y_val, y_val_hat, theta)
            record.append([i, loss, val_loss])
        # avg gradient per example
        gradient = (x.T.dot(y_hat - y) / sample_num) + penalization*theta
        # update
        theta = theta - lr * gradient
    return theta, record

结果绘图

这一步我们通过在训练结果以及在训练中保存的结果分别绘制train loss, val loss对比图，与拟合结果图。
代码如下：

def plotLoss(record):
    x = [it[0] for it in record]
    train_loss = [it[1] for it in record]
    val_loss = [it[2] for it in record]

    plt.figure(figsize=(8, 4))
    plt.plot(x, train_loss, label='$train-loss$', color='green', linewidth=0.5)
    plt.plot(x, val_loss, label='$val-loss$', color='red', linewidth=0.5)

    plt.plot(x, train_loss, 'go', markerfacecolor='none')
    plt.plot(x, val_loss, 'ro', markerfacecolor='none')

    plt.xlabel('Iterations')
    plt.ylabel('Loss')
    plt.title('Train vs Val')
    plt.legend()
    plt.show()

def plotNow(x, y, target_func, cur_func, theta, poly_degree):
    xrange = np.arange(0, 1, 0.01)
    targetfunc = target_func(2*np.pi*xrange)
    re_xrange = regenX(xrange, poly_degree)
    curfunc = cur_func(re_xrange, theta)

    plt.figure(figsize=(8, 4))
    plt.plot(xrange, targetfunc, label='$sin(2πx)$', color='green', linewidth=0.5)
    plt.plot(xrange, curfunc, label='$Polynomial$', color='red', linewidth=0.5)

    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Polynomial')
    # plt.xlim(0,1)
    plt.ylim(-2,2)
    plt.legend()
    plt.plot(x, y, 'bo', markerfacecolor='none')
    plt.show()