吴恩达机器学习课后习题——线性回归

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path = "ex1data1.txt"
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
# pandas.read_csv默认情况下会把数据内容的第一行默认为字段名标题，header=None指明读取的原始文件数据没有列索引
data.head() # 方法用于返回数据帧或序列的前n行(默认值为5)。返回类型：DataFrame，可理解为一个表格，或者说矩阵
# DataFrame本质上是二维矩阵，跟常规二维矩阵的差别在于前者额外指定了每一行和每一列的名称。

data.describe()
# describe() 函数可以查看数据的基本情况，
# 包括：count 非空值数、mean 平均值、std 标准差、max 最大值、min 最小值、（25%、50%、75%）分位数等。

# 画出图形
data.plot(kind='scatter',x='Population',y='Profit',figsize=(12,8))
# scatter表示散点图，figsize表示图片大小，单位：英寸
plt.show()

​

​

使用梯度下降来最小化代价函数,
创建一个以参数θ为特征函数的代价函数
J ( θ ) = 1 2 m ∑ i = 1 m ( h θ ( x ( i ) ) − y ( i ) ) 2 Jleft( theta right)=frac{1}{2m}sumlimits_{i=1}^{m}{{{left( {{h}_{theta }}left( {{x}^{(i)}} right)-{{y}^{(i)}} right)}^{2}}} J(θ)=2m1​i=1∑m​(hθ​(x(i))−y(i))2
其中：
h θ ( x ) = θ T X = θ 0 x 0 + θ 1 x 1 + θ 2 x 2 + . . . + θ n x n {{h}_{theta }}left( x right)={{theta }^{T}}X={{theta }_{0}}{{x}_{0}}+{{theta }_{1}}{{x}_{1}}+{{theta }_{2}}{{x}_{2}}+...+{{theta }_{n}}{{x}_{n}}\ hθ​(x)=θTX=θ0​x0​+θ1​x1​+θ2​x2​+...+θn​xn​

def computeCost(X,y,theta): # 输入X是列向量，y也是列向量，theta是行向量
    inner = np.power(((X*theta.T)-y),2) # x*θ的转置就是假设函数 power函数计算numpy数组中每个数值的指定幂次
    return np.sum(inner/(2*len(X))) # 将数组/矩阵中的元素全部加起来 len(X)矩阵行数

data.insert(0,'Ones',1)  
# 插入一列为1是为了更好的向量化，需要在训练集中添加一列x_0，以便我们可以使用向量化的解决方案来计算代价和梯度，见视频p18

# 变量的初始化
cols = data.shape[1] # shape[0]是行数  shape[1]是列数
X = data.iloc[:,0:cols-1] # 数据集 行全要  列从0到cols-1(不包括) 左闭右开
y = data.iloc[:,cols-1:cols] # 目标值
X.head()

获取的数据类型是DataFrame类型，所以需要进行类型转换。同时还需要初始化theta,即theta所有元素都设置为0

X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix(np.array([0,0]))

维度

X.shape, theta.shape, y.shape

((97, 2), (1, 2), (97, 1))

computeCost(X,y,theta) #计算代价函数

32.072733877455676

θ j : = θ j − α ∂ ∂ θ j J ( θ ) {{theta }_{j}}:={{theta }_{j}}-alpha frac{partial }{partial {{theta }_{j}}}Jleft( theta right) θj​:=θj​−α∂θj​∂​J(θ)
梯度下降

def gradientDescent(X,y,theta,alpha,iters): # iters是迭代次数 alpha是学习率 即步长
    temp_theta = np.matrix(np.zeros(theta.shape)) # 零值矩阵 暂存theta
    parameters = int(theta.ravel().shape[1]) # ravel计算需要求解的参数个数   功能将多维数组降至一维 这里值为2
    cost = np.zeros(iters) #构建iters个0的数组
    #迭代：
    for i in range(iters):
        difference = (X*theta.T) - y # 差值
        for j in range(parameters):
            term = np.multiply(difference,X[:,j])    #乘以x_i  X[:,j]表示所有行 j列，简言之就是抽出第j列  
             # 注意这里就不是矩阵乘法了 这里是求导后要求乘的 multiply是两个相同大小矩阵对应位置的数直接相乘
            temp_theta[0,j] = theta[0,j] - (alpha/len(X))*np.sum(term) # 更新theta_j
        
        theta = temp_theta # 更新全部theta值
        cost[i] = computeCost(X,y,theta) # 更新代价值
    
    return theta,cost

初始化一些附加变量 - 学习速率α和要执行的迭代次数。

alpha = 0.01
iters = 1000

获取theta值，代价最小值

g,cost = gradientDescent(X,y,theta,alpha,iters)
g

matrix([[-3.24140214,  1.1272942 ]])

cost[-1] # 获取数组的最后一个值

4.515955503078913

最后，用我们拟合过的theta值计算训练模型的代价函数

computeCost(X, y, g)

4.515955503078913

绘制线性模型以及数据，直观地看出它的拟合

x = np.linspace(data.Population.min(), data.Population.max(), 100) # 生成指定范围内指定个数的一维数组
f = g[0, 0] + (g[0, 1] * x) # x1=1  g[0, 0]表示theta_0 g[0, 1]表示theta_1

fig, ax = plt.subplots(figsize=(12,8)) # plt.subplot()函数用于直接指定划分方式和位置进行绘图。
# fig代表绘图窗口(Figure)；ax代表这个绘图窗口上的坐标系(axis)，一般会继续对ax进行操作。

# 设置横纵坐标的函数，并设置颜色为红色，在图标上标上'Prediction'标签  横坐标x 纵坐标f
ax.plot(x, f, 'r', label='Prediction') 
# 定义散点图的横纵坐标，并给出散点图点的'Traning Data'标签名
ax.scatter(data.Population, data.Profit, label='Traning Data') 
ax.legend(loc=2)  # 2表示在左上角


ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

​

​

由于梯度方程式函数也在每个训练迭代中输出一个代价的向量，所以我们也可以绘制。 注意，代价总是降低

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

​

​

该练习包括一个房屋价格数据集，其中包含两个变量（房屋大小、卧室数量）和目标（房子价格）。

path =  'ex1data2.txt'
data2 = pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])
data2.head()

data2.describe()

房屋大小约为卧室数量的1000倍。当特征相差几个数量级时，首先执行特征缩放（均值归一化）可以使梯度下降更快地收敛。

data2 = (data2 - data2.mean()) / data2.std() # 解决的方法是尝试将所有特征尺度都尽量缩放到 -1 ~ 1
data2.head()

data2.insert(0, 'Ones', 1) 
cols = data2.shape[1]
X2 = data2.iloc[:,0:cols-1]
y2 = data2.iloc[:,cols-1:cols]

X2 = np.matrix(X2.values)
y2 = np.matrix(y2.values)
theta2 = np.matrix(np.array([0,0,0]))

g2, cost2 = gradientDescent(X2, y2, theta2, alpha, iters)
computeCost(X2, y2, g2) # 计算代价

0.1307033696077189

代价函数收敛图

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost2, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

​

​

使用sklearn线性回归函数简化过程

from sklearn import linear_model
model = linear_model.LinearRegression()  #基于最小二乘法的线性回归。
model.fit(X, y)

x = np.array(X[:, 1].A1)
f = model.predict(X).flatten()

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

正规方程是通过求解下面的方程来找出使得代价函数最小的参数的： ∂ ∂ θ j J ( θ j ) = 0 frac{partial }{partial {{theta }_{j}}}Jleft( {{theta }_{j}} right)=0 ∂θj​∂​J(θj​)=0 。
假设我们的训练集特征矩阵为 X（包含了 x 0 = 1 {{x}_{0}}=1 x0​=1）并且我们的训练集结果为向量 y，则利用正规方程解出向量 θ = ( X T X ) − 1 X T y theta ={{left( {{X}^{T}}X right)}^{-1}}{{X}^{T}}y θ=(XTX)−1XTy 。
上标T代表矩阵转置，上标-1 代表矩阵的逆。设矩阵 A = X T X A={{X}^{T}}X A=XTX，则： ( X T X ) − 1 = A − 1 {{left( {{X}^{T}}X right)}^{-1}}={{A}^{-1}} (XTX)−1=A−1

梯度下降：需要选择学习率α，需要多次迭代，当特征数量n大时也能较好适用，适用于各种类型的模型
正规方程：不需要选择学习率α，一次计算得出，需要计算 ( X T X ) − 1 {{left( {{X}^{T}}X right)}^{-1}} (XTX)−1，如果特征数量n较大则运算代价大，因为矩阵逆的计算时间复杂度为 O ( n 3 ) O(n3) O(n3)，通常来说当 n n n小于10000 时还是可以接受的，只适用于线性模型，不适合逻辑回归模型等其他模型

# 正规方程函数
def normalEqn(X, y):
    theta = np.linalg.inv(X.T@X)@X.T@y#X.T@X等价于X.T.dot(X)  np.linalg.inv()：矩阵求逆
    return theta

final_theta2=normalEqn(X, y)
final_theta2  # 梯度matrix([[-3.24140214,  1.1272942 ]])

matrix([[-3.89578088],
        [ 1.19303364]])