最近一直在学习机器学习的几个入门模型,在看了吴恩达教授的机器学习部分视频后(在网易云课堂可以免费看所有视频,良心厂商),终于开始动手实现我的第一个机器学习算法,线性回归模型。
废话不多说 贴代码(完整代码在文末) (先贴代码 再讲解结构 可能理解有误 欢迎大家互相纠正学习)
下面代码是数据读取与数据预处理。
path = 'ex1data1.txt' #其实就是路径名称 我是放在同一个文件夹下的 所以直接写名字
data = pd.read_csv(path,names=['populations','profits']) #读取这个csv格式的数据文件
data.insert(0,'Ones',1) #在数据的第0列 插入一列 1
cols = len(data.columns) #将数据中的第0列到第n-1列记为X 最后一列记为 y
x_cols = data.columns[0:cols-1]
y_cols = data.columns[cols-1:cols]
X = data[x_cols] #获取数据集的数据
y = data[y_cols] #获取训练集的结果
X= np.matrix(X.values) #转换成numpy矩阵
y = np.matrix(y.values)
我使用的是吴恩达机器学习里面提到的数据,数据文件和代码文件会挂出来,大家可以直接理解为一个97行 2列的数据 列代表参数数量 populations代表参数1,profits代表这个房子的价格。
这里我们以profits为纵轴,populations为横轴绘制散点图,看看数据的分布。
可以看出,数据的分布是线性相关的,所以我们使用线性函数来进行拟合。
也就是说 共有2列, populations列为影响房子价格的因素,profits列是房屋价格,假设函数为y = k0 + k1x1 + k2x2 +...+knxn,当然我们这里是最简单的单变量,所以假设函数是 y = k0 + k1x
我们最终求的参数theta就是 k0 和 k1
我们这里先插入一列将矩阵构造成如下形式,后面我会讲解##注释1##这一步的意义,插入完后矩阵如下。
alpha = 0.01 #学习率
maxcyles = 1000 #迭代次数
theta = np.matrix(np.array([1,1])) #[[1 1]] 创建一个一行两列的元素均为1的矩阵 实际上代表 系数 k0 + k1
接下来就是重头戏,核心函数的部分了,我们一共用到了两个函数,分别是 COST(X,y,theta) ,gradienDescent(X,y,theta,maxcyles)
到了这里我们就需要把几个公式摆出来了
首先是我们的假设函数 也就是一个简单的一元一次方程
其次代价函数
更新公式
有了这三个公式 就可以贴代码了
COST(X,y,theta)
def COST(X,y,theta): #求代价函数
pow = np.power((X*theta.T - y),2)
return np.sum(pow)/(2*len(X))
简单解释下函数运算过程 我们分别看看X矩阵 与theta的转置矩阵即theta.T的形式
X矩阵(97行2列) theta.T矩阵(2行1列)
所以我们可以简单的理解 临时变量 pow = ( (97*2)矩阵* (2*1)矩阵)的平方
在线性代数中我们知道,一个m*n的矩阵乘以一个n*1的矩阵,结果是一个m行1列的矩阵
所以我们打印一下pow看下 果然是个 m * 1的矩阵
那么最后返回的就是简单求个和 再乘以1/2m就行了。
打印一下这个返回值 np.sum(pow)/(2*len(X))
确实是一个数值,COST函数实际上就是求出第i行的参数与theta中各个系数相乘以后的结果 减去 对应y的第i行的结果 代表的是当前系数值下,预测的第i行与训练集第i行结果的误差。
gradienDescent(X,y,theta,maxcyles)
def gradienDescent(X,y,theta,maxcyles):
temp = np.matrix(np.zeros(theta.shape)) #创建一个和theta同样维度的零矩阵
parameters = theta.shape[1] #参数的数目
cost = np.zeros(maxcyles) #创建一个1000长度的cost误差数组
for i in range(maxcyles):
error = X*theta.T - y
for j in range(parameters):
term = np.multiply(error,X[:,j])
temp[0,j] = theta[0,j] - (alpha/len(X))*np.sum(term)
theta = temp
cost[i] = COST(X,y,theta)
return theta,cost
简单讲解一下
我们首先创建一个临时变量temp 它和theta是一模一样的 等于复制出了一个新的theta 用来在更新时作为临时变量
这里很关键 但博主也了解的不是很深刻 只知道在吴恩达教授的视频中,他提到过,参数更新时,是要同时更新的,也就是说参数theta[0,0]和theta[0,1]是要同时更新的,我认为这就是为什么需要这个临时变量temp 而不是直接写成这个形式
theta[0,j] = theta[0,j] - (alpha/len(X))*np.sum(term) 之后如果我理解的更深刻我会修改此处的理解 如果有网友对这里有更好的理解 欢迎评论 我会及时修改
term = np.multiply(error,X[:,j])
temp[0,j] = theta[0,j] - (alpha/len(X))*np.sum(term) #结合更新公式理解
这两步是关键的更新步骤 我们刚刚已经贴出了我们的更新公式 现在对应公式来理解每一项
代码中的term变量即为
而 temp[0,j] = theta[0,j] - (alpha/len(X))*np.sum(term)对应的就是
即更新公式
最后完成更新 并记录每次迭代的损失或者说代价
迭代完成后,我们便获得了一个theta 这个theta就是一个1*n的矩阵 n为系数个数
我们最终是通过得到theta[0,0] 与theta[0,1]即 y = k0 + k1x 的 系数k0 , k1
y = theta[0,0] + theta[0,1]*x
这样我们就可以将给的测试集的数据通过假设函数的到一个结果矩阵即预测结果集
最后我们看下拟合的效果
代价函数的迭代图像
可以观察到 整个损失是呈递减的趋势,且已经收敛。
下面是整个线性回归的代码
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def COST(X,y,theta): #求代价函数
pow = np.power((X*theta.T - y),2)
return np.sum(pow)/(2*len(X))
def gradienDescent(X,y,theta,maxcyles):
temp = np.matrix(np.zeros(theta.shape)) #创建一个和theta同样维度的零矩阵
parameters = theta.shape[1] #参数的数目
cost = np.zeros(maxcyles) #创建一个1000长度的cost误差数组
for i in range(maxcyles):
error = X*theta.T - y
for j in range(parameters):
term = np.multiply(error,X[:,j])
temp[0,j] = theta[0,j] - (alpha/len(X))*np.sum(term)
theta = temp
cost[i] = COST(X,y,theta)
return theta,cost
##### 数据读取与数据预处理 #####
path = 'ex1data1.txt' #其实就是路径名称 我是放在同一个文件夹下的 所以直接写名字
data = pd.read_csv(path,names=['populations','profits']) #读取这个csv格式的数据文件
data.insert(0,'Ones',1) #在数据的第0列 插入一列 1
cols = len(data.columns) #将数据中的第0列到第n-1列记为X 最后一列记为 y
x_cols = data.columns[0:cols-1]
y_cols = data.columns[cols-1:cols]
X = data[x_cols]
y = data[y_cols] #获取训练集的数据和结果
X= np.matrix(X.values)
y = np.matrix(y.values)
##### 变量初始化 #####
alpha = 0.01 #学习率
maxcyles = 1000 #迭代次数
theta = np.matrix(np.array([1,1])) #[[1 1]] 创建一个一行两列的元素均为1的矩阵 实际上代表 系数 k0 + k1
##### 将theta丢进去迭代得到预测结果与y误差最小的系数矩阵theta #####
g,cost = gradienDescent(X,y,theta,maxcyles)
##### x的数据范围 #####
x = np.linspace(data.populations.min(),data.populations.max(),100)
##### 假设函数 #####
f = g[0,0] + g[0,1]*x
##### 先画出数据的散点图看看数据分布 #####
fig,ax = plt.subplots(figsize = (12,8))
ax.scatter(data.populations,data.profits)
ax.set_xlabel('populations')
ax.set_ylabel('profits')
plt.show()
##### 画出假设函数的拟合效果 #####
fig,ax = plt.subplots(figsize = (12,8))
ax.plot(x,f,'r')
ax.scatter(data.populations,data.profits)
ax.set_xlabel('populations')
ax.set_ylabel('profits')
plt.show()
##### 绘制代价函数的迭代图像 #####
fig,ax = plt.subplots(figsize = (12,8))
ax.plot(np.arange(maxcyles),cost,'b')
plt.show()
这里是数据文件(我看看怎么放进来)