今天与大家一步一步手动实现基础的线型回归模型，不调包实现。

线型回归模型框架包括的基本组件示意图如下所示：

定义模型基本组件

下面手动编码分别实现以上基本组件。

先导入本次实战使用的包：

import time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D 
from matplotlib import cm 

线性模型定义

# 定义线性模型
def model(theta,X):
    theta = np.array(theta)
    return X.dot(theta)

定义损失函数

# 定义损失函数
def cost(m,theta,X,y):
    'print(theta)'
    ele = y - model(theta,X)
    item = ele**2
    item_sum = np.sum(item)
    return item_sum/2/m

求梯度

### 求梯度
def gradient(m,theta,X,y,cols):
    grad_theta = []
    for j in range(cols):
        grad = (y-model(theta,X)).dot(X[:,j])
        grad_sum = np.sum(grad)    
        grad_theta.append(-grad_sum/m)

    return np.array(grad_theta)

参数更新

# 参数更新
def theta_update(grad_theta,theta,sigma):
    return theta - sigma * grad_theta

训练停止准则

# 训练停止准则
def stop_stratege(cost,cost_update,threshold):
    return cost-cost_update < threshold

搭建模型

OLS，全称 ordinary least square，搭建 OLS模型的过程示意图如下所示：

使用第 1 节介绍的基本组件和上面的示意图，搭建 OLS 模型也就不会太困难，完整代码如下所示：

# OLS模型
def OLS(X,y,threshold):
    start = time.time()

    m=100 # 样本个数
    theta = [0,0,0] # 设置权重参数的初始值
    iters = 0; # 迭代步数
    cost_record=[] # 记录迭代过程中损失函数值
    sigma = 0.0001 # 学习率

    # 模型初始损失函数值
    cost_val = cost(m,theta,X,y)
    cost_record.append(cost_val)

    while True:
        # 模型迭代
        grad = gradient(m,theta,X,y,3)
        theta = theta_update(grad,theta,sigma)
        cost_update = cost(m,theta,X,y)
        if stop_stratege(cost_val,cost_update,threshold):
            break

        if iters%2000==0:
            print('iter:{}, loss:{:.6f}'.format(iters,cost_update))

        iters=iters+1
        cost_val = cost_update
        cost_record.append(cost_val)

    end = time.time()
    print('迭代次数:%d'%(iters))
    print("OLS收敛历经时长: %f s" % (end - start))
    print("参数收敛值：")
    print(theta)

    return cost_record, iters,theta

应用模型

使用 OLS 模型实战的基本过程如下所示，首先生成训练使用的数据，然后调用模型，最后得到结果、分析模型和绘图分析等。

创建数据

# 生成实战所需要的数据
def build_data(num_points = 100):
    x1 = np.random.normal(0.,0.55,size=(num_points,))
    x2 = np.random.normal(0.,0.75,size=(num_points,))
    y12 = x1*2.0 + x2*1.2 + np.random.normal(0.,0.01)
    data = np.c_[x1,x2,y12]

    X = data[:,[0,1]] # 两列特征(theta_1, theta_2; 不包含theta_0)
    y = data[:,2] # 标签值
    return X, y

默认创建 100 个含有两个特征的数据样本。如下调用创建数据的函数。因为线性回归模型会有一个截距列，所以使用 np.c_ 连接上这一列。

X, y = build_data()
b = np.repeat(1,100) #这就是参数theta_0,初始值都设置为 1
X = np.c_[b,X] # 扩展连接上theta_0特征列

绘制生成的数据样本，代码如下所示：

# 绘制数据点
def show_points(X, y):
    plt.figure(figsize=(12,8))
    plt.scatter(X[:,1],y,color='blue')
    plt.xlabel('X1')
    plt.ylabel('y')
    plt.show()

    plt.figure(figsize=(12,8))
    plt.scatter(X[:,2],y,color='green')
    plt.xlabel('X2')
    plt.ylabel('y')
    plt.show()

    fig=plt.figure(figsize=(12,8))
    ax = Axes3D(fig)
    ax.scatter(X[:,1],X[:,2],y,c='red')
    ax.set_zlabel('y') #坐标轴
    ax.set_ylabel('X2')
    ax.set_xlabel('X1')
    plt.show() 

调用绘制样本点的函数：

show_points(X,y)

第一个特征与 y 关系散点图：

第二个特征与 y 关系散点图：

第一、二个特征与 y 散点图：

打印显示前 5 个样本点：

X[:5]

结果：

array([[ 1.        ,  1.16544904,  0.85327085],
       [ 1.        ,  0.19225326, -0.47360726],
       [ 1.        , -0.47874836,  1.2603231 ],
       [ 1.        , -0.72399195,  0.40869013],
       [ 1.        , -0.07579486,  0.18751382]])

调用 OLS 模型

模型的损失值设为 1e-10：

costs, iters,theta = OLS(X,y,1e-10)

模型收敛后得到损失函数值，迭代次数，模型参数值 theta。每 2000 时步打印一次损失值，如下打印前 10000 时步模型的损失值：

iter:0, loss:1.205221
iter:2000, loss:0.993409
iter:4000, loss:0.822327
iter:6000, loss:0.683437
iter:8000, loss:0.570121
iter:10000, loss:0.477226

模型在 178000 时步收敛，OLS 收敛历经时长 9.023159 s，截距和两个特征参数收敛值分别为 [-0.00388289 1.99734008 1.1997594]。

iter:178000, loss:0.00001
迭代次数:178282
OLS收敛历经时长: 9.023159 s
参数收敛值：
[-0.00388289  1.99734008  1.1997594 ]

结果分析

训练得到的参数值，与我们理论模型设定的 x1*2.0 + x2*1.2 + np.random.normal(0.,0.01) 已经非常接近，验证确实已经收敛。

绘制时步与损失函数值的折线图：

plt.figure(figsize=(12,8))
plt.plot(costs,color='blue')
plt.xlabel('iter')
plt.ylabel('loss value')
plt.show()

绘制收敛时得到的线性回归曲线：

fig=plt.figure(figsize=(12,8))
ax = Axes3D(fig)
ax.scatter(X[:,1],X[:,2],y,c='red')


#绘制两个特征和y的三维图
theta2=theta[1:].reshape(2,1)
x1=np.linspace(-2,2,100)
x2=x1
xtmp=np.c_[x1,x2]
yn =xtmp.dot(theta2) #拟合值
ax.plot(xtmp[:,0],xtmp[:,1],yn.flatten(),c='blue')

ax.set_zlabel('yn') #坐标轴
ax.set_ylabel('X2')
ax.set_xlabel('X1')
plt.show()

plt.figure(figsize=(12,8))
plt.scatter(X[:,1],y,color='blue')
plt.plot(x1,x1*theta[1])
plt.xlabel('X1')
plt.ylabel('y')
plt.show()

plt.figure(figsize=(12,8))
plt.scatter(X[:,2],y,color='green')
plt.plot(x2,x2*theta[2])
plt.xlabel('X2')
plt.ylabel('y')
plt.show()

小结

今天与大家完整手动不调包实现了 OLS 模型，并且使用 OLS 模型创建数据，并且训练得到参数。整个过程都配有示意图。相信大家通过今天的学习已经掌握基本的机器学习模型创建和实践的过程。其实，对于后面我们会学到的更加复杂的模型，整个过程也大同小异。

恭喜大家，已经入门机器学习。下面我们继续探索机器学习其他更多的实用模型。

附今天课程 Notebook 的下载链接：

https://pan.baidu.com/s/1DoLTfnvOR63rCcRUpPrxMQ

提取码：ke89