机器学习——单变量线性回归模型

本文主要探讨一元线性回归模型，假设函数，损失函数和梯度下降，先来看一下原数据长什么样？

• 原数据

ex1data1.txt有两列数据，第一列是自变量(x)，第二列为因变量(y)，都是连续性变量，维度是97*2，（如果你需要数据可以去三行科创微信公众号交流群要）首先，我们要对数据进行初探式的探索。

• 探索性分析

绘制出因变量和自变量的散点图，观察散点图呈现主要关系是符合什么模式？

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
data = pd.read_csv("CourseraML/ex1/data/ex1data1.txt", header = None, names = ['area', 'price']) #读取数据并加上表头
x = np.array(data['area']) #自变量
y = np.array(data['price'])#因变量
m = len(x) #样本量
fig = plt.figure(figsize = (6,4)) #建立画布
plt.scatter(x, y, marker = 'x') #散点图
plt.xlabel("Population of City in 10,000s") #横坐标标签
plt.ylabel("Profit in $10,000s") #纵坐标标签

• 定义假设函数和损失函数

从散点图，我们看到自变量和因变量有比较好的线性关系，故可以进行假设，假设函数为

$$h_{\theta }(x)={\theta }^{T}x$$

损失函数
损失函数主要是指预测值和真实值的差异，可以定义总体损失函数如下

$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^i)-y^i{)}^2$$

x = x.reshape(m,1) #重塑成97*1的矩阵形式
y = y.reshape(m,1)
X= np.insert(x,0,1, axis =1) #增加常数项参数，np.insert()将向量插入某一行或列,形成97*2维
initial_theta = np.zeros((X.shape[1],1)) #构造初始theta值, 参数theta个数等于特征的个数2*1维

def h(theta, X): #定义假设函数h =theta^TX
    return np.dot(X, theta)  #97*2 x 2*1 = 97*1

def costFunction(mytheta, X, y): #定义损失函数
    return float(1./(2*m)*np.dot((h(mytheta, X)-y).T, (h(mytheta, X)-y)))
print(costFunction(initial_theta, X,y)) #32.07

• 定义梯度下降函数，并可视化损失函数

梯度下降是一种迭代的参数求解方法，其梯度下降方程为

$${\theta }_j:={\theta }_j-\frac{\alpha }{m}(h_{\theta }(x^{(i)}-y^{(i)})x_j^{(i)}$$

alpha = 0.01
max_iteration = 2000

def gradientDescent(X, theta_start = np.zeros(2)): #定义梯度下降函数
    theta = theta_start
    jvec = [] #用来存放损失函数值
    thetahistory = [] #用来存放theta历史值
    for iter in range(max_iteration):
        tmptheta = theta
        jvec.append(costFunction(theta,X,y))
        thetahistory.append(list(theta[:,0]))
        for j in range(len(tmptheta)): #同步更新theta参数
            tmptheta[j] = theta[j]- (alpha/m)*np.sum((h(theta, X)-y)*np.array(X[:,j]).reshape(m, 1))
        theta = tmptheta
    return theta, thetahistory, jvec
initial_theta = np.zeros((X.shape[1],1)) #构造初始theta值, 参数theta个数等于特征的个数2*1维

theta, thetahistory, jvec = gradientDescent(X,initial_theta)

def plotConvengence(jvec):
    plt.figure(figsize = (6,4))
    plt.plot(range(len(jvec)), jvec)
    plt.title("Convengence of cost function")
    plt.xlabel("Iteration")
    plt.ylabel("Cost function")
    plt.ylim([4,6]) #纵坐标压缩到4到6范围

plotConvengence(jvec)

• 绘制最佳拟合曲线

通过梯度下降法在获得最佳参数后，可以得到具体的假设函数，从而可以绘制出拟合曲线

def myfit(xval): # 绘制最佳拟合曲线
    return theta[0]+theta[1]*xval

plt.figure(figsize=(6,4))
plt.scatter(x,y,marker='x', label ="Training Data")
plt.plot(x, myfit(x), color ='r', label = "Hypothesis:h(x)=%0.2f+%0.2fx"%(theta[0],theta[1]))
plt.xlabel("Profit in $10,000s")
plt.ylabel("Population of city in 10,000s")
plt.legend()

• 损失函数的可视化

损失函数是关于参数${\theta }_0,{\theta }_1$的函数，在3维的情况下可以绘制出梯度下降损失函数变化情况

from mpl_toolkits.mplot3d import axes3d, Axes3D
from matplotlib import cm
import itertools

fig = plt.figure(figsize =(16,16))
ax =fig.gca(projection = '3d')
xvals = np.arange(-10, 10, 0.5)
yvals = np.arange(-1, 4, 0.1)
myxs, myys,myzs = [], [], []
for david in xvals:
    for kaleko in yvals:
        myxs.append(david)
        myys.append(kaleko)
        myzs.append(costFunction(np.array([[david],[kaleko]]), X, y))
scat = ax.scatter(myxs, myys, myzs, c =np.abs(myzs), cmap = plt.get_cmap('YlOrRd'))
plt.xlabel(r"$\theta_0$")
plt.ylabel(r"$\theta_1$")
plt.title("GradientDescent")
plt.plot([x[0] for x in thetahistory], [x[1] for x in thetahistory], jvec, 'bo-', marker = "x")
plt.show()