用python不用框架来手写一个线性回归的模型｜Pyth

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

目标

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

一点小要求

熟悉简单高等数学，例如求导
会用 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] \rightarrow [1,x]

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')
复制代码

找到一个函数集

接下里来就是找到一个函数集，函数形式为

y = mx + c

h_{\theta}(x) = \Theta^Tx = \theta_1x + \theta_0

\theta_1

$表示函数线性方程的斜率，$

\theta_1

$，这里$

\Theta = \begin{bmatrix} \theta_0 & \theta_1 \end{bmatrix}

h_{\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(\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

是更新梯度的步幅的大小。也被称为学习率。

\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])
复制代码

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')
复制代码

预测

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

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>]
复制代码

用python不用框架来手写一个线性回归的模型｜Pyth

目标

一点小要求

获取数据集

找到一个函数集

代价函数(损失函数)

参数初始值

开始训练

超参数

预测

近期文章

近期评论

标签

热门

文章归档

分类目录

功能