用 python 不用框架来手写一个线性回归的模型|Python 主题月

1,877 阅读4分钟

本文正在参加「Python主题月」,详情查看活动链接

目标

仅用 python 以及其常用操作数组的库 numpy 和可视化库 matplotlib 来完成一个线性模型,从读取数据集、处理数据、找到函数集、设定目标函数到通过梯度下降算法来更新参数优化损失函数一步一步实现,让你了解实现一个模型的全部过程。在整个过程中,对于一些关键点或者难点给予必要说明,如果文章中出现偏差或者不够详细,欢迎留言来进行探讨。

002.png

一点小要求

  • 熟悉简单高等数学,例如求导
  • 会用 python 写代码
  • 了解一点 numpy 和 matplotlib 而已
%matplotlib inline
# import libs
import numpy as np
import matplotlib.pyplot as plt

获取数据集

通常我们是先分析数据集,数据只有两个字段,数据集比较简单,仅有 97 条记录,每条记录有 2 两列,第 1 列作为样本属性,而第 2 列作为样本的标签,我们就是学一个函数表示从特征(x)到标签(y)的映射。需要数据集的人可以给我在评论区留言,我会将数据集发送给你。

datafile = 'data/dataset.txt'

读取数据用到 np.loadtxt 这个方法,参数可以指定字段间分格符号,每行记录不同字段用逗号分隔 , unpack=True 表示按列读取

# load dataset comma , first col is x second is y(label)
cols = np.loadtxt(datafile,delimiter=',',usecols=(0,1),unpack=True)

因为 unpack=True 表示按列读取,所以读取数据是每行表示一个属性,每列对应一个记录,所以这里用 cols 变量来接受数据,表示是按列读取数据

cols.shape()#(2, 97)
(2, 97)

接下来就是将数据集第 1 行(也就是每个样本第一列属性值)做样本 X,然后将数据集第 2 行提取出做样本的标签 y 来使用

X = np.transpose(np.array(cols[:-1]))
X.shape # (97, 1)
# 形状为 [[x],[],[]]
(97, 1)
y = np.transpose(np.array(cols[-1:]))
y.shape #(97, 1)
(97, 1)
# 样本数量
m = y.size
m
97

这样做的原因是为向每一个记录添加一个维度,让每一个记录变为增广向量[x][1,x][x] \rightarrow [1,x]X = np.insert(X,0,1,axis=1) 解释一下,现在 X 形式是[[x1],[x2],..., [xn]] 用 insert 方法想轴 1 也就是每条记录第一个 0 位置添加一个 1 然后就编程[[1,x1],[1,x2],..., [1,xn]]

X = np.insert(X,0,1,axis=1)
X.shape#(97, 2)
(97, 2)
X[:2]
array([[1.    , 6.1101],
       [1.    , 5.5277]])

我们查看一下 X 和 y 的形状,然后利用 matplotlib 显示一下样本点在坐标系分别情况

X.shape, y.shape
((97, 2), (97, 1))
plt.figure(figsize=(10,6))
plt.plot(X[:,1],y[:,0],'rx',markersize=10)
plt.grid(True)
plt.ylabel('y')
plt.xlabel('x')
Text(0.5, 0, 'x')

output_17_1.png

找到一个函数集

接下里来就是找到一个函数集,函数形式为 y=mx+cy = mx + c 形式,将的输出是由x、m 和 c 的线性组合,我们进一步表示为 hθ(x)=ΘTx=θ1x+θ0h_{\theta}(x) = \Theta^Tx = \theta_1x + \theta_0,其中 θ1\theta_1 表示函数线性方程的斜率,θ1\theta_1 ,这里 Θ=[θ0θ1]\Theta = \begin{bmatrix} \theta_0 & \theta_1 \end{bmatrix}

hθ=θTx=θ0+θ1xh_{\theta} = \theta^T x = \theta_0 + \theta_1 x

iterations = 1500
alpha = 0.01
# y = \theta_0 x + \theta_1
def hypothesis_fxn(theta,X):
#     (97,2),(2,1)
    return np.dot(X,theta)

代价函数(损失函数)

代价函数是用来衡量找到参数所对应的函数的误差。在代价函数中,ℎ𝜃(𝑥(𝑖)) 是模型的输出值。而 y(i)是真实值。用于计算预测值和真实值之间差距。这是一个线性例子。我们知道,(ℎ𝜃(𝑥(𝑖))-𝑦(𝑖))的值可以是正的或者负的。通过对计算误差取平方然后求和可以抵消掉负数值。还可以求平均数。回归的代价函数有一些变量。如:平均误差(MSE),平均绝对误差(MAE)等。 J(θ)=12mi=1m(hθ(x(i))y(i))2J(\theta) = \frac{1}{2m} \sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)})^2

def cost_fxn(theta, X, y):
    h = hypothesis_fxn(theta, X)
#     print(h)
    d = h-y
    
    c = 1/(2 * m)
    loss = c * np.sum(d **2)
    return loss

参数初始值

将模型要学习的参数(权重、偏置)都初始化为 0 然后我们计算一下在没有进行任何学习情况下损失函数值

# \theta(2,1)
intial_theta = np.zeros((X.shape[1],1))
intial_theta.shape
(2, 1)
cost_fxn(intial_theta,X,y)
32.072733877455676

开始训练

我们从一个简单的参数值(通常是0)开始,计算梯度,使用当前参数的梯度值来更新参数。为了让代价函数最小化,所以计算参数相对代价函数相对于的导数。𝐽(𝜃)相对于𝜃的导数是调整参数方向,而大小又和 α\alpha 有关, α\alpha是更新梯度的步幅的大小。也被称为学习率。

θj=θjα1mi=1m(h(θ)(x(i))y(i))x(i)\theta_j = \theta_j - \alpha \frac{1}{m} \sum_{i=1}^m (h_{(\theta)}(x^{(i)}) - y^{(i)})x^{(i)}

超参数

  • 迭代次数 iterations
  • 学习率 alpha
def gradient_descent(X, theta=np.zeros(2)):
#     cost list plot cost 
    costs = []
    theta_history = []
    for i in range(iterations):
        temp_theta = theta
#         cal loss cost
        c = cost_fxn(theta, X,y)
    
        costs.append(c)
        
        theta_history.append(list(theta[:,0]))
        
        for j in range(len(temp_theta)):
#         
            temp_theta[j] = theta[j] - (alpha / m) * np.sum((hypothesis_fxn(theta,X) -y )*np.array(X[:,j]).reshape(m,1))
        theta = temp_theta
    return theta, theta_history,costs
        
costs 和 theta_history 用于将每次迭代中损失值和参数值保存用于绘制图像,通过观察图像来了解训练过程。
initial_theta = np.zeros((X.shape[1],1))
initial_theta
array([[0.],
       [0.]])
list(initial_theta[:,0])
[0.0, 0.0]
theta, thetahisotry,jvec = gradient_descent(X, initial_theta)
jvec = np.array(jvec).reshape(-1,1)
def plot_convergence(jvec):
    plt.figure(figsize=(10,6))
    plt.plot(range(len(jvec)),jvec,'co')
    
    plt.grid(True)
    plt.title('Convergence of Cost Function')
    plt.xlabel("Iteration number")
    plt.ylabel("cost function")
    dummy = plt.xlim([-0.05 * iterations, 1.05 * iterations])
plot_convergence(jvec)
# dummy = plt.ylim([4,7])

output_36_0.png

thetas = np.array(thetahisotry)
plt.plot(thetas)
plt.grid(True)
plt.legend([r"$\Theta0$",r"$\Theta1$"])
plt.ylabel(r"$\Theta")
Text(0, 0.5, '$\\Theta')

output_37_1.png

预测

训练好模型之后我们可以绘制一些模型对应直线,然后观察这条直线拟合样本情况

def prediction(theta,x):
    return theta[0] + theta[1] * x
pred = prediction(theta,X[:,1])
plt.plot(X[:,1],pred)
plt.plot(X[:,1],y[:,0],'rx',markersize=10,label='Training Data')
[<matplotlib.lines.Line2D at 0x124e36a90>]

output_41_1.png