最小二乘法拟合直线（问题i）

问题描述

有一个函数，现在不知道函数的的具体形式，给定满足函数关系的一组训练样本，请使用线性回归模型拟合出函数y=f(x)。
(可尝试一种或几种不同的基函数，如多项式、高斯或sigmoid基函数）

要求：

i. 先完成最小二乘法的优化 (参考书中第二章 2.3节中的公式)

ii. 附加题：实现“多项式基函数”以及“高斯基函数”（可参考PRML）

iii. 附加题：完成梯度下降的优化 (参考书中第二章 2.3节中的公式)
使用训练集train.txt 进行训练，使用测试集test.txt 进行评估（标准差），训练模型时请不要使用测试集。

读取数据

我们先将采集到的数据个读出来：这里以txt的格式保存，利用numpy库从文本中读取数据（注意jupyter notebook中数据必须保存到所建立）：

import numpy as np #导入numpy库
import matplotlib.pyplot as plt 
train_dataset = np.loadtxt('train.txt')
test_dataset = np.loadtxt('test.txt')
#print(len(train_dataset))
#print(len(test_dataset))
#print(train_dataset)
#print(test_dataset)

在读出所给的数据之后，我们按照题目所给的要求完成最小二乘法的优化（参照书上2.3节公式）这里单单我们看书是很难理解的，因为合理给出了优化参数的方法，但是我拿到这道题目的时候被难住的是最开始的基本线性回归的知识。题目中提到了用到基函数（内心很是mengbi）下面将一些基础知识整理一下，然后自己再试着编程看一看。

线性回归（线性基函数）

最简单的线性回归模型的形式如下所示：
$$y\left(\mathbf{x},\mathbf{w}\right)=w_0+w_{1}x+\cdot \cdot \cdot +w_n{x}_n$$
其中$\mathbf{x}={\left(x_1,x_2,\cdots ,x_n\right)}^T$,$y$是输入变量$x_1$的线性组合，这种模型具有很大的局限性。
一般用向量形式写成：$$f\left(x\right)={\mathbf{w}}^{T}x+b$$
其中$\mathbf{w}=\left( w_1,w_2,\cdots ,w_n \right) $,$w$和$b$学得之后，模型就可以确定。

我们可以将其扩展为输入变量通过有限个非线性函数转换后的线性组合，如下式表示：$$y\left(\mathbf{x},\mathbf{w}\right)=w_0+\sum _{j=1}^{n-1}{w}_j{\varphi }_j\left(\mathbf{x}\right)$$
这里$\phi \left( \mathbf{x} \right) $就是基函数。我们可以添加一个基函数$\phi 0=1$，则上式就可以改为$$y\left( \mathbf{x,w} \right) =\sum{j=0}^{M-1}{w_j\phi _j\left( \mathbf{x} \right)}=\mathbf{w}^T\phi \left( \mathbf{x} \right) $$
这里$\mathbf{w}={\left(w_0,w_1,\cdots w_{n-1}\right)}^T$
,$\varphi ={\left({\varphi }_0,{\varphi }_1,\cdots {\varphi }_{n-1}\right)}^T$基函数$\phi _j\left( x \right) $可以有多种不同类型，常见的基函数有：

·多项式基函数：$${\varphi }_j\left(x\right)=x^j$$

·高斯基函数：$${\varphi }_j\left(x\right)=e^{\left\{-\frac{{\left(x-u_j\right)}^2}{2s^2}\right\}}$$
这里$u$控制位置$s$控制跨度

·$sigmoid$基函数$${\varphi }_j\left(x\right)=\sigma \left(\frac{x-u_j}{s}\right)$$
$$\sigma \left(a\right)=\frac{1}{1+e^{\left(-a\right)}}$$

参数估计

一般的情况是数据集$D_i$样本是由$d$个属性描述，此时我们试图学得$$f\left(x_i\right)={\mathbf{w}}^T{x}_i+b$$使得$f\left(x_i\right)\backsimeq y_i$,这称为“多元线性回归”(multivariate linear reegression)

类似的我们可以用最小二乘法对$w$和$b$进行估计，为了便于讨论我们将$w$和$b$吸入向量形式$\mathbf{\hat{w}}=\left( \mathbf{w:b} \right) $，相应的，把数据集$D$表示为一个$m×(d+1)$大小的矩阵$X$,其中每行对应于一个示例表示为下述矩阵形式：
$$\mathbf{X}=\left[\begin{array}{ccccc}l{x}_{11}& x_{12}& \cdots & x_{1d}& 1\\ x_{21}& x_{22}& \cdots & x_{2d}& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ x_{m1}& x_{m2}& \cdots & x_{md}& 1\end{array}\right]=\left[\begin{array}{cc}l{x}_1^T& 1\\ x_2^T& 1\\ ⋮& ⋮\\ x_m^T& 1\end{array}\right]$$
当$X^{T}X$为满秩矩阵的时候或者为正定矩阵的时候令上式为零我们可以很快速的求得参数。$${\hat{w}}^{\ast }={\left(X^{T}X\right)}^{-1}{X}^{T}y$$
这里我们可以很快的求得多元线性回归模型为：$$f\left({\hat{x}}_i\right)={\hat{x}}_i^T{\left(X^{T}X\right)}^{-1}{X}^{T}y$$然而现实中$X^{T}X$往往不是满秩矩阵，例如很多任务中我们会遇到大量的变量，其数目超过了样例数，导致$X$的列数多于行数，$X^{T}X$显然不是满秩矩阵此时有多个$\hat{w}$，它们都能够是均方差最小化，这里将引入正则化(regularization)项。这里将先不论述正则化，待以后再写一篇正则化的博客。

下面我们根据上述的要求来进行编程实现。

按照最小二乘法进行处理

from numpy import mat
#先将训练集x数据导出来做成矩阵的形式
x_train_dataset=train_dataset[:,0]
#将一维数组转换为矩阵
x_train_dataset=mat(x_train_dataset).reshape(300,1)
a = np.ones(300)
x_train_dataset=np.column_stack((x_train_dataset,a))
#print(x_train_dataset)
print(x_train_dataset.shape)

#将训练集y数据导出来
y_train_dataset=train_dataset[:,1]
y_train_dataset=mat(y_train_dataset).reshape(300,1)
#print(y_train_dataset.shape)
#print(y_train_dataset)

#将训练集x进行转置
xt_train_dataset=np.transpose(x_train_dataset)
print(xt_train_dataset.shape)
#print(xt_train_dataset)

(300, 2)
(2, 300)

# 利用公式求参数w
w_predict = np.matmul(xt_train_dataset,x_train_dataset)
#print(w_predict)
w_predict = np.linalg.inv(w_predict)
#print(w_predict)
w_predict = np.matmul(w_predict,xt_train_dataset )
w_predict = np.matmul(w_predict,y_train_dataset )
print(w_predict)

[[0.94732599]
 [0.62680273]]

我们已经计算出参数，现在利用matplotlib将其可视化(散点图和直线图)

#重新读取数据
x_train_dataset=train_dataset[:,0]
y_train_dataset=train_dataset[:,1]
#画散点图
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.scatter(x_train_dataset,y_train_dataset,color='red',marker='o',s=1)
ax.plot(x, 0.94732599*x+0.62680273,color='green',ls='-')
#画直线图

上数我们计算出所需要的参数，这里利用测试集去测试所得参数的效果。我们计算出偏差和方差来对此种算法进行评估：

x_test_dataset=test_dataset[:,0]
#print(y_test_dataset)
x_test_dataset=mat(x_test_dataset).reshape(200,1)
#print(y_test_dataset)
y_test_dataset=test_dataset[:,1]
y_test_dataset=mat(y_test_dataset).reshape(200,1)
#print(y_test_dataset)

import math
bias2 = 0
#测试
#bias2=bias2+math.pow((x_test_dataset[1]*0.94732599+0.62680273)-y_test_dataset[1],2)
#print(bias2)
for i in range (0,200):
    bias2=bias2+math.pow((x_test_dataset[i]*0.94732599+0.62680273)-y_test_dataset[i],2)
print(bias2)
#之前看吴恩达的视频其中讲到要尽量避免用到loop语句，
#下面我将对其进行修正

926.8174513840322