简单线性回归的梯度下降与正则方法实现

一、单变量情况
首先我们先调用需要的方法，我们希望数据在算法中以数组或矩阵进行数据处理，引入numpy。我们对数据进行预处理时，需要使用pandas。对数据进行可视化时，需要使用matplotlib。所以首先调用需要使用的模块，并导入数据，查看数据形式。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path='F:\MachineLearning\data\ex1data1.txt'
data=pd.read_csv(path,header=None,names=['population','profit'])    #不以任何行作为列索引，利用参数`names`，可以自定义列索引。
data.head()                                                         #展示数据形式

我们可以看到数据形式：

接下来将数据集X与标记向量分开y，并查看X的情况：

m,d=data.shape[0],data.shape[1]-1
X=data.iloc[:,0:d]
X.insert(0,'Ones',1)   #由于theata中常数项的存在，必须插入全是1的列
y=data.iloc[:,d:d+1]   #对X，y，即数据集和标记向量分开并进行处理

X.head()               #查看X的情况
y.head()               #查看y的情况

数据可视化，以我们命名的'population'列为横坐标，profit为纵坐标，画散点图：

data.plot(kind='scatter',x='population',y='profit',color='red') #数据可视化
plt.show()

得到散点图：

最后，我们把X，y转化为np.array形式，并定义初始theta值：

X=np.array(X.values)       #将数据处理为numpy数组数据，并查看其维度
y=np.array(y.values)
theata=np.zeros((1,2))

接下来进行线性回归最核心部分，代价函数与梯度下降法：
首先，编写代价函数

def computeCost(X,y,theata):    #编写代价函数
    sub_Vector=np.dot(X,theata.T)-y
    sqr_sub=sub_Vector**2
    tot_error=sqr_sub.sum(axis=0)
    avg_error=tot_error/(2*m)
    return avg_error[0]

这里的代价函数形式为：

如果以矩阵形式则有：

接着时梯度下降算法：

def gradientDescent(X,y,alpha,theata,iters):    #梯度下降算法
    cost=[]
    for i in range(iters):
        sub_Vector=np.dot(X,theata.T)-y
        temp=np.dot(sub_Vector.T,X)/m
        theata=theata-alpha*temp
        error=computeCost(X,y,theata)
        cost.append(error)
    return theata,cost

这里我们对theata用向量形式一次性表示，不必一次计算一个theata值，减少一次循环。就是利用以下矩阵等式：

利用以下的迭代式对theta进行更新：

我们设置迭代步数为2000，alpha为0.01：

alpha=0.01
iters=2000
g,cost=gradientDescent(X,y,alpha,theata,iters)
print(g)

g即为theta向量，输出为：

[[-3.78806857  1.18221277]]

然后我们把散点图和我们计算出的theta值的线性回归函数画在一起，观察你和情况：

x=np.linspace(data.population.min(),data.population.max(),100)
f=g[0,0]+g[0,1]*x
fig,ax=plt.subplots()
ax.plot(x,f,color='red',label='predict')
ax.scatter(data.population,data.profit,label='real data')
ax.legend()
ax.set_xlabel('population')
ax.set_ylabel('profit')
ax.set_title('predit and real data')

所以有

可以看出拟合情况良好。当然，我们在编写梯度下降算法时，把每一步迭代所对应的损失函数取值都放在了cost列表里，借此我们可以通过画图大致了解其收敛情况：

x=np.arange(iters)     #代价函数衰减曲线（到2000步）
y=cost
fig,ax=plt.subplots(figsize=(12,8))
ax.plot(x,y,color='blue',label='cost function')
ax.set_xlabel('iters')
ax.set_ylabel('error')
ax.set_title('Iters VS Error')
plt.show()

可以看出2000步之后，损失函数的斜率趋近于0，即收敛情况良好。

二、多变量情况
如果我们采用的数据复杂一点，即标记由两个变量控制。
我们先观察数据情况：

path='F:\MachineLearning\data\ex1data2.txt'
data2=pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])
data2.head()

可以看出Size这一属性范围过大，采用梯度下降法时会需要很多迭代步数才能达到收敛，所以应采取归一化处理：

def autoNomal(data):
    data=(data-data.mean())/data.std()
    return data

接下来，对数据向上面一样进行预处理，把处理后的数据集，标记向量，初始化theta命名为X2，y2，theta2：

path='F:\MachineLearning\data\ex1data2.txt'
data2=pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])
data2=autoNomal(data2)
data2.head()

def text2data(data):                                 #数据处理函数
    m,d=data.shape[0],data.shape[1]-1
    data.insert(0,'ones',1)
    X=data.iloc[:,0:d+1]
    y=data.iloc[:,d+1:d+2]
    X=np.array(X.values)
    y=np.array(y.values)
    theata=np.zeros((1,d+1))
    return X,y,theata
X2,y2,theata2=text2data(data2)

将X2，y2和theta2代入梯度下降算法：

g2,cost2=gradientDescent(X2,y2,alpha,theata2,iters)
print(g)

得到最终的theta：
[[-1.03513990e-16 8.77575946e-01 -4.59890736e-02]]
这里需注意，我们如果要把新样品用于模型，根据参数得到初始标记后还需进行进一步处理，因为我开始对数据采用了归一化。
我们在看看损失函数情况：

x=np.arange(iters)     #代价函数衰减曲线（到2000步）
y=cost2
fig,ax=plt.subplots(figsize=(12,8))
ax.plot(x,y,color='blue',label='cost function')
ax.set_xlabel('iters')
ax.set_ylabel('error')
ax.set_title('Iters VS Error')
plt.show()

可以看出2000步时，损失函数趋近于收敛。

三、正规方法

正规方法从纯矩阵的角度给线性回归计算参数，其推导过程为：



代码实现为：

def normalEqn(X,y):
    A=np.dot(X.T,X)
    B=np.linalg.inv(A)
    C=np.dot(B,X.T)
    theata=np.dot(C,y)
    return theata

print(normalEqn(X,y))

得到theta为
[[-3.89578088]
[ 1.19303364]]
与我们用梯度下降法得到的答案相差不大。
那么什么时候采用梯度下降法，什么时候采用正规方法呢：只要特征变量的数目并不大，标准方程是一个很好的计算参数的替代方法。具体地说，只要特征变量数量小于一万，我通常使用标准方程法，而不使用梯度下降法。