单变量线性回归
个人blog https://verbf.github.io
在这部分的练习中,你将使用单变量的线性回归来预测食品卡车的利润。假设你是a公司的CEO并正在考虑在不同的城市开设一家连锁餐厅。该连锁店已经在多个城市拥有卡车,并且你拥有有关城市的人口和利润的数据。
你可以使用这些数据来帮助你选择接下来在那个城市发展。
文件ex1data1.txt包含了线性回归问题的数据集,其中第一列是城市人口,第二列是食品卡车在该城市的利润。
绘制数据图像
再开始任务之前,我们通过数据可视化来更好的理解数据。因为这个数据集只有两个属性,所以我们可以使用散点图来进行可视化。(我们在现实生活中遇到的很多问题往往是多维的,不能进行二维的绘制)
首先要载入数据,载入数据后我们将其打印在屏幕上
print("loading data ex1data1.txt...")
ex1data = loadtData('ex1data1.txt')
X = ex1data[0]
y = ex1data[1]
print(X)
print(y)
print()
loadData 函数实现
# 载入数据,返回一个二维的numpy数组,
# 第一维是 x轴数据,第二维是 y轴数据
def loadtData(file_path):
X = np.array([])
Y = np.array([])
for i in open(file_path):
# 根据逗号的位置取出数据
x = i[0:i.index(',')-1]
y = i[i.index(',')+1:len(i)-1]
# 读出的数据是字符串类型,需要转换为浮点类型
X = np.append(X,float(x))
Y = np.append(Y,float(y))
return np.array([X,Y])
绘制图像
接下来我们调用plotData 函数来绘制散点图,顺便设置一下横纵坐标的标题
x_label = "Population of City in 10,000s"#设置横纵坐标的标题
y_label = "Profit in $10,000s"
plotData(X,y,x_label,y_label)
plotData函数实现
# 将数据可视化,使用散点图
def plotData(X,y,x_label,y_label):
plt.scatter(X,y)
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.show()
绘制出来的结果应该与下图类似
梯度下降
在这一部分,你将使用梯度下降法来拟合单变量线性回归中的参数 θ \theta θ,使其与我们的数据集相符
更新方程
线性回归的目标是使代价函数达到最小值
J
(
Θ
)
=
1
2
m
∑
i
=
1
m
(
h
Θ
(
x
(
i
)
)
−
y
(
i
)
)
2
J(\Theta)=\frac{1}{2m}\sum_{i=1}^{m}{({h_\Theta}({x}^{(i)})-{y}^{(i)})}^{2}
J(Θ)=2m1i=1∑m(hΘ(x(i))−y(i))2
假设函数使用下面的线性模型
h
θ
(
x
)
=
θ
T
x
=
θ
0
+
θ
1
x
1
h_\theta(x) = \theta^Tx = \theta_0 +\theta_1x_1
hθ(x)=θTx=θ0+θ1x1
重新计算模型中参数
θ
\theta
θ的值,通过调整参数的值使代价函数的值最小化。我们使用批次梯度下降法来达到目的。在梯度下降法中,每一次迭代都会完成一次更新
θ
j
=
θ
j
−
α
1
m
∑
i
=
1
m
(
h
θ
(
x
(
i
)
)
−
y
(
i
)
)
x
j
(
i
)
\theta_j = \theta_j - \alpha\frac{1}{m}\sum_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}
θj=θj−αm1i=1∑m(hθ(x(i))−y(i))xj(i)
需要注意的是,所有的
θ
j
\theta_j
θj需要同时更新。
梯度下降的每一次更新都会使参数
θ
j
\theta_j
θj更接近最优的值,同时这也会使代价函数的值达到最小。
实现
在上面的步骤中,我们已经准备好了线性回归所需的数据。
接着,我们给数据增加一个维度,方便我们对参数
θ
0
\theta_0
θ0进行优化。
我们将参数
θ
\theta
θ都初始化为0,并设置学习率为0.01。
m = len(y)#样本数量
X = np.c_[np.ones(m),X]#为X增加一行 值为1
theta = np.zeros((2,1),float)#初始化参数theta
#一些梯度下降的设置
iterations = 1500 #迭代次数
alpha = 0.01 #学习率
计算代价函数
当我们使用梯度下降法来最小化代价函数的值时,我们可以通过计算代价函数的值来判断是否收敛。
我们接下来的任务就是实现computeCost()函数,该函数的功能就是计算代价函数的值。
当我们在实现该函数的时候,需要注意变量X,y是矩阵类型的变量,而不是标量。矩阵的行代表了训练集中的样本。
一旦你实现了这个函数,我们就使用初始化为0的参数
θ
\theta
θ来运行一次该函数。
该函数的运行结果应该是32.07
# 计算并显示初始的代价值
J = computeCost(X,y,theta)
print('With theta = [0 ; 0] ');
print("Cost computed = %f \n" % J)
print('Expected cost value (approx) 32.07\n');
# 继续测试代价函数
theta[0] = -1
theta[1] = 2
J = computeCost(X, y, theta);
print('\nWith theta = [-1 ; 2]\nCost computed = %f\n'% J);
print('Expected cost value (approx) 54.24\n');
computeCost函数实现
# 计算代价函数
def computeCost(X,y,theta):
m = len(y)
result = np.dot(X , theta)
result = result - y.reshape(97,1)
result = np.square(result)
result = np.sum( result,axis=0)
result = result/(2.0*float(m))
return result
梯度下降
接下来你要完成梯度下降的编码。
在你的程序中,你要清楚的理解优化目标是什么,什么是需要更新的。
你要记住,代价函数的参数是向量
θ
\theta
θ,而不是X,y。也就是说,我们需要更新向量
θ
\theta
θ的值来是代价函数最小化,而不是改变 X 或 y。如果你不确定的话,参考上面给出的方程,或是视频的课程。
我们可以通过观察代价函数的值在每一次更新中是否持续下降,以此来判断梯度下降法是否正常工作。
梯度下降在每一次迭代中都会调用computeCost函数,如果你正确实现了computeCost函数和梯度下降,那么你的代价函数的值绝不会增加,并且将会在算法的最后达到一个稳定的值。
print('\nRunning Gradient Descent ...\n')
# 运行梯度下降
theta = gradientDescent(X, y, theta, alpha, iterations);
# 将theta的值打印到屏幕上
print('Theta found by gradient descent:\n');
print('theta_0 : %f \ntheta_1 : %f'%(theta[0],theta[1])) ;
print('Expected theta values (approx)\n');
print(' -3.6303\n 1.1664\n\n');
gradientDescent 函数实现
# 梯度下降
def gradientDescent(X, y, theta, alpha, iterations):
m = len(y)
for i in range(0,iterations):
theta_0 = theta[0] - alpha * computePartialDerivative(X,y,theta,0)
theta_1 = theta[1] - alpha * computePartialDerivative(X,y,theta,1)
theta[0] = theta_0
theta[1] = theta_1
print( "iterations : ",i, " cost : ",computeCost(X,y,theta))
return theta
# 计算偏导数
def computePartialDerivative(X , y, theta , num):
m = len(y)
result = 0
for i in range(0,m):
result += (theta[0]*X[i][0] + theta[1]*X[i][1] - y[i])*X[i][num]
result /= m
return result
当你完成了以上任务时,使用最后得到的参数 θ \theta θ的值来绘制假设函数的图像,你会看到与下面类似的图
# 绘制线性拟合的图
plt.plot(X[:,1],np.dot(X,theta))
plt.scatter(X[:,1],y,c='r')
plt.show()
可视化 J( θ \theta θ)
为了更好的理解代价函数 J(
θ
\theta
θ) , 我们可以将
θ
0
\theta_0
θ0和
θ
1
\theta_1
θ1的值绘制在二维的网格上。
绘图代码
# 绘制三维的图像
fig = plt.figure()
axes3d = Axes3D(fig)
# 指定参数的区间
theta0_vals = np.linspace(-10, 10, 100)
theta1_vals = np.linspace(-1, 4, 100)
# 存储代价函数值的变量初始化
J_vals = np.zeros((len(theta0_vals), len(theta1_vals)))
# 为代价函数的变量赋值
for i in range(0,len(theta0_vals)):
for j in range(0,len(theta1_vals)):
t = np.zeros((2,1))
t[0] = theta0_vals[i]
t[1] = theta1_vals[j]
J_vals[i,j] = computeCost(X, y, t)
# 下面这句代码不可少,matplotlib还不熟悉,后面填坑
theta0_vals, theta1_vals = np.meshgrid(theta0_vals, theta1_vals) #必须加上这段代码
axes3d.plot_surface(theta0_vals,theta1_vals,J_vals, rstride=1, cstride=1, cmap='rainbow')
plt.show()
完整代码
最后附上完整代码
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 载入数据,返回一个二维的numpy数组,
# 第一维是 x轴数据,第二维是 y轴数据
def loadtData(file_path):
X = np.array([])
Y = np.array([])
for i in open(file_path):
# 根据逗号的位置取出数据
x = i[0:i.index(',')-1]
y = i[i.index(',')+1:len(i)-1]
# 读出的数据是字符串类型,需要转换为浮点类型
X = np.append(X,float(x))
Y = np.append(Y,float(y))
return np.array([X,Y])
# 将数据可视化,使用散点图
def plotData(X,y,x_label,y_label):
plt.scatter(X,y)
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.show()
# 计算代价函数
def computeCost(X,y,theta):
m = len(y)
result = np.dot(X , theta)
result = result - y.reshape(97,1)
result = np.square(result)
result = np.sum( result,axis=0)
result = result/(2.0*float(m))
return result
# 梯度下降
def gradientDescent(X, y, theta, alpha, iterations):
m = len(y)
for i in range(0,iterations):
theta_0 = theta[0] - alpha * computePartialDerivative(X,y,theta,0)
theta_1 = theta[1] - alpha * computePartialDerivative(X,y,theta,1)
theta[0] = theta_0
theta[1] = theta_1
print( "iterations : ",i, " cost : ",computeCost(X,y,theta))
return theta
# 计算偏导数
def computePartialDerivative(X , y, theta , num):
m = len(y)
result = 0
for i in range(0,m):
result += (theta[0]*X[i][0] + theta[1]*X[i][1] - y[i])*X[i][num]
result /= m
return result
#======================== 绘图 ====================================
print("loading data ex1data1.txt...")
ex1data = loadtData('ex1data1.txt')
X = ex1data[0]
y = ex1data[1]
print(X)
print(y)
print()
x_label = "Population of City in 10,000s"#设置横纵坐标的标题
y_label = "Profit in $10,000s"
#绘图
plotData(X,y,x_label,y_label)
#======================= 代价函数 和 梯度下降 ======================
m = len(y)#样本数量
X = np.c_[np.ones(m),X]#为X增加一行 值为1
theta = np.zeros((2,1),float)#初始化参数theta
#一些梯度下降的设置
iterations = 1500 #迭代次数
alpha = 0.01 #学习率
print("Testing the coust function ...\n")
# 计算并显示初始的代价值
J = computeCost(X,y,theta)
print('With theta = [0 ; 0] ');
print("Cost computed = %f \n" % J)
print('Expected cost value (approx) 32.07\n');
# 继续测试代价函数
theta[0] = -1
theta[1] = 2
J = computeCost(X, y, theta);
print('\nWith theta = [-1 ; 2]\nCost computed = %f\n'% J);
print('Expected cost value (approx) 54.24\n');
print('\nRunning Gradient Descent ...\n')
# 运行梯度下降
theta = gradientDescent(X, y, theta, alpha, iterations);
# 将theta的值打印到屏幕上
print('Theta found by gradient descent:\n');
print('theta_0 : %f \ntheta_1 : %f'%(theta[0],theta[1])) ;
print('Expected theta values (approx)\n');
print(' -3.6303\n 1.1664\n\n');
# 绘制线性拟合的图
plt.plot(X[:,1],np.dot(X,theta))
plt.scatter(X[:,1],y,c='r')
plt.show()
# 绘制三维的图像
fig = plt.figure()
axes3d = Axes3D(fig)
# 指定参数的区间
theta0_vals = np.linspace(-10, 10, 100)
theta1_vals = np.linspace(-1, 4, 100)
# 存储代价函数值的变量初始化
J_vals = np.zeros((len(theta0_vals), len(theta1_vals)))
# 为代价函数的变量赋值
for i in range(0,len(theta0_vals)):
for j in range(0,len(theta1_vals)):
t = np.zeros((2,1))
t[0] = theta0_vals[i]
t[1] = theta1_vals[j]
J_vals[i,j] = computeCost(X, y, t)
# 下面这句代码不可少,matplotlib还不熟悉,后面填坑
theta0_vals, theta1_vals = np.meshgrid(theta0_vals, theta1_vals) #必须加上这段代码
axes3d.plot_surface(theta0_vals,theta1_vals,J_vals, rstride=1, cstride=1, cmap='rainbow')
plt.show()