AI-机器学习-第一课

# 使用for循环实现
def my_dot(a, b):
    x=0
    for i in range(a.shape[0]):
        x = x + a[i] * b[i]
    return x

# 测试
a = np.array([1, 2, 3, 4])
b = np.array([-1, 4, 3, 2])
# 打印结果my_dot(a, b) = 24
print(f"my_dot(a, b) = {my_dot(a, b)}")

# 使用NumPy的dot函数
a = np.array([1, 2, 3, 4])
b = np.array([-1, 4, 3, 2])
c = np.dot(a, b)

# 简单的测试for循环和NumPy的点积运算速度
np.random.seed(1)
# 生成
a = np.random.rand(10000000)
b = np.random.rand(10000000)
# 使用dot函数
# 起始时间
tic = time.time()
c = np.dot(a, b)
# 结束时间
toc = time.time()
print(f"np.dot(a, b) =  {c:.4f}")
print(f"Vectorized version duration: {1000*(toc-tic):.4f} ms ")

# 使用for循环
tic = time.time()
c = my_dot(a,b)
toc = time.time()
print(f"my_dot(a, b) =  {c:.4f}")
print(f"loop version duration: {1000*(toc-tic):.4f} ms ")
# 释放大数组的内存
del(a);del(b)

# 二维数组
X = np.array([[1],[2],[3],[4]])
w = np.array([2])
c = np.dot(X[1], w)
# 打印结果
# X[1] has shape (1,)
# w has shape (1,)
# c has shape ()
print(f"X[1] has shape {X[1].shape}")
print(f"w has shape {w.shape}")
print(f"c has shape {c.shape}")

# 表示创建包含5个元素的1维数组
a = np.zeros((1, 5))
# 打印结果为a shape = (1, 5), a = [[0. 0. 0. 0. 0.]]
print(f"a shape = {a.shape}, a = {a}")        
# 表示创建包含1个元素的2维数组
a = np.zeros((2, 1))
# 打印结果为a shape = (2, 1), a = [[0.]
# [0.]]
print(f"a shape = {a.shape}, a = {a}")
# 生成含有1个元素的1维随机数数组
a = np.random.random_sample((1, 1))
# 打印结果为a shape = (1, 1), a = [[0.44236513]]
print(f"a shape = {a.shape}, a = {a}")

# 也可使用如下方法创建 
a = np.array([[5], [4], [3]]);

# 使用reshape可方便将1维数组切分为多维数组
# reshape第1个参数表示分成几个数组,如果为-1则表示通过计算自动求得,如果指定了,若和计算求得的不一致,则报错
# reshape第2个参数表示分成的数组中元素的数量
# 下面的式子等价于a = np.arange(6).reshape(3, 2)
a = np.arange(6).reshape(-1, 2)
# 打印结果为a.shape: (3, 2), 
# a= [[0 1]
#  [2 3]
#  [4 5]]
print(f"a.shape: {a.shape}, \na= {a}")

# 访问多维数组的元素,如果是2维数组,则使用两个下标
# 打印结果为4
print(f"\na[2,0] = {a[2,0]}")

# 访问多维数组中的一行,使用一个下标
# 打印结果为a[2].shape:   (2,), a[2]   = [4 5]
print(f"a[2].shape:   {a[2].shape}, a[2]   = {a[2]}")

# 多维数组的切分操作
# 2维数组,数组行数2个,每个数组有10个元素
a = np.arange(20).reshape(-1, 10)
# 打印结果为
# a = 
# [[ 0  1  2  3  4  5  6  7  8  9]
#  [10 11 12 13 14 15 16 17 18 19]]
print(f"a = \n{a}")
# 访问第一行的数组(下标为0), 切分规则和一维数组相同
print("a[0, 2:7:1] = ", a[0, 2:7:1], ",  a[0, 2:7:1].shape =", a[0, 2:7:1].shape, "a 1-D array")
# 访问多维数组中所有数组, 并且进行切分
print("a[:, 2:7:1] = \n", a[:, 2:7:1], ",  a[:, 2:7:1].shape =", a[:, 2:7:1].shape, "a 2-D array")
# 访问多维数组中的所有元素
print("a[:,:] = \n", a[:,:], ",  a[:,:].shape =", a[:,:].shape)
# 访问下标为1的数组的所有元素
print("a[1,:] = ", a[1,:], ",  a[1,:].shape =", a[1,:].shape, "a 1-D array")
# a[1,:]可以简写为a[1]
print("a[1]   = ", a[1],   ",  a[1].shape   =", a[1].shape, "a 1-D array")

X_train = np.array([[2104, 5, 1, 45], [1416, 3, 2, 40], [852, 2, 1, 35]])
y_train = np.array([460, 232, 178])

b_init = 785.1811367994083
# 矢量w有四个特征
w_init = np.array([ 0.39133535, 18.75376741, -53.36032453, -26.42131618])

def compute_cost(X, y, w, b): 
    """
    成本函数计算
    参数:
      X: 特征数组
      y: 目标数据数组
      w: 模型参数w数组
      b: 模型参数b数字
    返回:
      cost: 成本
    """
    m = X.shape[0]
    cost = 0.0
    for i in range(m):
        # 使用点积运算快速计算函数f
        f_wb_i = np.dot(X[i], w) + b
        cost = cost + (f_wb_i - y[i])2
    cost = cost / (2 * m)
    return cost
# 进行调用
cost = compute_cost(X_train, y_train, w_init, b_init)
# 打印结果为Cost at optimal w : 1.5578904045996674e-12
print(f'Cost at optimal w : {cost}')

梯度下降

$$\begin{align*} \text{repeat}&\text{ until convergence:} ; \lbrace \newline; & w_j = w_j - \alpha \frac{\partial J(\mathbf{w},b)}{\partial w_j} \tag{5} ; & \text{for j = 0..n-1}\newline &b\ \ = b - \alpha \frac{\partial J(\mathbf{w},b)}{\partial b} \newline \rbrace \end{align*}$$

其中导数部分 $$ \begin{align} \frac{\partial J(\mathbf{w},b)}{\partial w_j} &= \frac{1}{m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})x_{j}^{(i)} \tag{6} \ \frac{\partial J(\mathbf{w},b)}{\partial b} &= \frac{1}{m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)}) \tag{7} \end{align} $$ 公式(6), (7)代码:

def compute_gradient(X, y, w, b):
    # m记录数据数量, n是特征数量
    m, n = X.shape
    dj_dw = np.zeros((n,))
    dj_db = 0.
    
    for i in range(m):
        # 计算都使用到的部分
        item = (np.dot(X[i], w) + b) - y[i]
        for j in range(n):
            dj_dw[j] = dj_dw[j] + item * X[i, j]
        dj_db = dj_db + item
    dj_dw = dj_dw / m
    dj_db = dj_db / m
    
    return dj_dw, dj_db

公式(5)代码

def gradient_descent(X, y, w_in, b_in, cost_fuction, gradient_function, alpha, num_iters):
    """
    参数alpha是学习率
    num_iters是迭代次数
    
    返回值有:
      w: 特征参数数组
      b: 参数b的值
      J_history: 计算函数J的历史值数组
    """
    J_history = []
    w = copy.deepcopy(w_in)
    b = b_in
    
    for i in range(num_iters):
        # 取出计算的线性导数值
        dj_dw, dj_db = gradient_function(X, y, w, b)
        w = w - alpha * dj_dw
        b = b - alpha * dj_db
        
        if i < 100000:
            J_history.append(cost_function(X, y, w, b))
        
        if i % math.ceil(num_iters/10) == 0:
            print(f"Iteration {i:4d}: Cost {J_history[-1]:8.2f}   ")
    return w, b, J_history

  
# 调用和测试
initial_w = np.zeros_like(w_init)
initial_b = 0.
# 循环次数
iterations = 1000
# 学习率
alpha = 5.0e-7
# 调用
w_final, b_final, J_hist = gradient_descent(X_train, y_train, initial_w, initial_b, compute_cost, compute_gradient, alpha, iterations)

# 显示最终的参数b和w
print(f"b,w found by gradient descent: {b_final:0.2f},{w_final} ")
m,_ = X_train.shape
# 打印使用训练好的参数m和b和实际的值进行比较
for i in range(m):
    print(f"prediction: {np.dot(X_train[i], w_final) + b_final:0.2f}, target value: {y_train[i]}")

特征缩放

为何使用特征缩放

当有不同的特征, 且它们的取值范围非常不同时, 可能会导致梯度下降运行缓慢. 如果重新缩放不同的特征, 使它们都具有可比较的取值范围, 从而让运行速度, 升级, 特征显著.

如果预测房价的模型是如下图:

尺寸和房价数量的取值范围相差巨大, 直观感受图:

通过上图可看出, 特征不同会造成数据点不均匀, 也导致成本函数J的等高线图的同心曲线不怎么接近圆形, 一侧较短, 一侧较长.

对w1的非常小的变化会对估计价格产生非常大的影响, 也对成本函数J的影响大, 因为w1的往往会乘以一个非常大的数(特征的取值范围大).

相比之下, w2需要更大的变化才能大大改变预测的价格, 对成本函数的影响比较小.

这样会导致:

等高线又高又瘦, 梯度下降可能会在它最终找到全局最小值之前来回弹跳很长时间.

使用特征缩放后, 梯度下降可以找到一条更直接的通往全局最小值的路径.

特征缩放的几种方式

最大值法

使用特征的取值范围的最大值, 作为分布, 对特征进行缩放, 如下图

均值归一化

Mean normalization

找到特征的平均值$\mu$, 使用公式$x = \frac {x - \mu} {max - min}$进行缩放

Z-Score归一化

要实现z-score归一化, 需要计算每个特征的标准差$\sigma$和平均值$\mu$, 使用公式$x = \frac {x - \mu} {\sigma}$进行缩放

如何计算标准差?

使用图上所示的正态分布

代码

实现Z-Score归一化

$$x^{(i)}_j = \dfrac{x^{(i)}_j - \mu_j}{\sigma_j} \tag{4}$$

其中$\mu_j$和$\sigma_j$由以下公式获得 $$ \begin{align} \mu_j &= \frac{1}{m} \sum_{i=0}^{m-1} x^{(i)}j \ \sigma^2_j &= \frac{1}{m} \sum{i=0}^{m-1} (x^{(i)}_j - \mu_j)^2 \end{align} $$

def zscore_normalize_features(X):
    # 把训练集中的所有特征数据都进行计算
    # axis=0为纵向,相当于取二维数组中的一维数组中的下标相同的元素进行操作
    # axis=1为横向,相当于取二维数组中的一维数组中所有元素进行操作
    mu = np.mean(X, axis=0)
    sigma = np.std(X, axis=0)
    X_norm = (X - mu) / sigma
    
    return (X_norm, mu, sigma)

度量特征缩放

缩放后的特征数据在一个数量级时是可接受的, 如下图

梯度下降收敛

如何判断梯度下降工作正确

通过观察迭代次数和成本函数J的曲线图, 成本函数J应该在每次迭代后减少. 如果增加, 意味着学习率$\alpha$选择不当或者代码中可能存在错误.

如果如上图那样平滑的下降, 并且在迭代到某个次数时不再有大的变化, 则此时梯度下降工作正确.

如上图, 可判断出在迭代400次后, 成本函数J收敛(converged). 不同的应用程序的收敛的迭代次数千差万别.

还有一个使用自动收敛测试的方法(但不推荐):

给定一个阈值$\epsilon$, 如果每次迭代成本函数的减少小于此变量, 则可声明收敛. 但选择正确的阈值$\epsilon$非常困难.

选择正确的学习率

如果发现成本函数J不是随着迭代次数下降, 则有可能是学习率选择不当, 或是代码有错误, 如何区分呢?

正确实施梯度下降的一个调试技巧是, 如果学习率足够小, 成本函数应该在每次迭代中减少. 如果梯度下降不起作用, (调试经验)将学习率设置为一个非常小的数字, 看是否在每次迭代时导致成本下降.

如果将学习率设置为非常小的数字, 成本函数J也不会在每次迭代中减少, 而有时会增加, 那么通常意味着代码中某处存在错误(如图中的把减号写成加号).

依次尝试学习率为...0.001, 0.01, 0.1, 1... 如果找到了让成本函数J下降到上升的边界, 则尝试用最小学习率. 通过每次乘以3倍的重试, 找到比较符合成本函数下降最快的学习率.

特征工程

Feature engineering

对于许多实际应用而言, 选择或输入正确的特征是使算法良好运行的关键.

可以利用对问题的知识或直觉来设计新特征, 通常是转换或组合问题的原始特征来使学习算法更容易做出准确的预测.

英文定义:

Using intuition(直觉) to design new features, by transforming or combining original features.

多项式回归

Polynomial regression

为了更切合训练数据, 可能会使用多项式, 而二次幂会随着特征值的增加而降低, 这显然是不符合逻辑的, 因此加入了三次幂来预测.

如果是图上的多次项, 那么特征缩放会越来越重要, 因为特征数据的范围会随着多次幂急剧增加.

还有一种合理的替代方案, 使用特征的平方根, 平方根函数随着特征的增加, 曲线不会那么陡峭, 但永远不会变平, 而且肯定永远不会回落.

使用Scikit-Learn

文档地址: https://scikit-learn.org/stable/index.html

代码

# 导入库
import numpy as np
from sklearn.linear_model import LinearRegression, SGDRegressor
from sklearn.preprocessing import StandardScaler

# 加载训练集
X_train, y_train = load_house_data()

# 特征缩放
scaler = StandardScaler()
X_norm = scaler.fit_transform(X_train)

# 创建和调整回归模型
sgdr = SGDRegressor(max_iter=1000)
sgdr.fit(X_norm, y_train)

# 训练出来的模型参数
b_norm = sgdr.intercept_
w_norm = sgdr.coef_

# 做出预测,y_pred_sgd是预测值,y_pred是目标值
y_pred_sgd = sgdr.predict(X_norm)
y_pred = np.dot(X_norm, w_norm) + b_norm  

分类算法

只有两种可能输出的分类问题称为二元分类(binary classification).

通常将否定结果的称为negative class, 正确结果的称为positive class.

线性回归不能很好的运用于分类问题, 如下图:

如上图使用阈值0.5来预测是否是肿瘤. 阈值0.5称为决策边界(decision boundary).

如图所示, 如果右边添加了一个训练数据, 会导致线性模型函数f向右倾斜, 决策边界向右移动. 导致原来预测为恶性的, 现在变为良性.

使用逻辑回归(logistic regression)算法解决分类问题.

逻辑回归

Logistic regression

描述

使用sigmoid 函数(也叫逻辑函数-logistic function)解决分类问题.

$g(z) = \frac{1}{1+e^{-z}}$ 其中$e$是自然参数, 此函数的曲线如上图所示.

若有线性模型f, 使$z = \mathbf{w} \cdot \mathbf{x} + b$ 代入到sigmoid函数中, 可得: $$ f_{\mathbf{w},b}(\mathbf{x}) = g(\mathbf{w} \cdot \mathbf{x} + b ) = \frac{1}{1+e^{-(\mathbf{w} \cdot \mathbf{x} + b)}} $$ 逻辑回归的输出描述:

在给定输入特征数据的情况下输出类或标签等于1(正确结果)的概率.

概率公式: $P(y = 0) + P(y = 1) = 1$, 其中P是概率.

假如输出的结果是0.8, 则说明有80%的几率有恶性肿瘤的风险.

有时会表示为如下公式: $$ f_{\mathbf{w},b}(\mathbf{x}) = P(y=1|\mathbf{x};\mathbf{w},b) $$ 公式描述: f的值为给定输入的特征x和参数w, b, y=1的概率.

代码

sigmoid函数的实现: $g(z) = \frac{1}{1+e^{-z}}$

# 导入NumPy库
import numpy as np

# 使用np.exp()函数获取自然参数e的值
# 但只有一个输入值,且为1时
input_val = 1
# 相当于e的input_val幂
exp_val = np.exp(input_val)
# 输出结果为Output of exp: 2.718281828459045
print("Output of exp:", exp_val)
# 传入数组
input_array = np.array([1,2,3])
exp_array = np.exp(input_array)
# 输出结果为Output of exp: [ 2.72  7.39 20.09]
print("Output of exp:", exp_array)

# 实现sigmoid函数
def sigmoid(z):
    g = 1 / (1 + np.exp(-z))
    return g

决策边界

Decision boundary

描述

如上图, 如果选择阈值为0.5, 则当模型函数$f >= 0.5$ 时, 预测的y值为1(yes, 正确结果). 有以下公式: $$ z = \mathbf{w} \cdot \mathbf{x} + b \newline f_{\mathbf{w},b}(\mathbf{x}) = g(z) = \frac{1}{1+e^{-(\mathbf{w} \cdot \mathbf{x} + b)}} $$ 也就是$g(z) >= 0.5$ 于是得出 $z >= 0$ 因此有 $\mathbf{w} \cdot \mathbf{x} + b >= 0$

可视化

决策边界就是让 $z = \mathbf{w} \cdot \mathbf{x} + b = 0$ 的曲线

上图时两个特征的决策边界, 是一条直线.

上图时二次多项式下的两个特征的决策边界, 是一个圆.

上图时多次多项式下的多个特征的决策边界, 是复杂的曲线.

逻辑回归的成本函数

对于逻辑回归使用误差平方成本函数(squared error cost)不是很好的选择, 如下图:

如果使用误差平方成本函数, 则在逻辑回归中, 曲线如上图所示, 则成本函数J将会有很多的局部最小值. 不是一个凸函数曲线.

图中的 $\frac{1}{2}(f_{w,b}(x^{(i)}) - y^{(i)})^2$ 称作单个训练数据中的损失(loss), 记作:$L(f_{w,b}(x^{(i)}), y^{(i)})$

损失函数衡量的是在一个训练样例的表现如何, 成本函数衡量的是在整个训练集是那个的表现.

通过给损失函数选择不同的形式, 将能够保持总成本函数成为一个凸函数.

公式

上图展示的是损失函数当 $y^{(i)} = 1$ 的情况

函数$-log(f)$的曲线如图所示, 因模型函数f的取值范围是[0, 1], 因此只看此函数$-log(f)$在区间[0, 1]的曲线.

如果模型函数f的值趋近于1, 则损失函数趋近于0.

如果模型函数f的值取经于0, 则损失函数趋近于无穷.

因此当 $y^{(i)} = 1$ 时, 损失函数会帮助算法做出更准确的预测, 因为当预测值接近1时, 此时的损失最小.

上图展示的是损失函数当 $y^{(i)} = 0$ 的情况

函数$-log(1-f)$的曲线如图所示, 因模型函数f的取值范围是[0, 1], 因此只看此函数$-log(1-f)$在区间[0, 1]的曲线.

如果模型函数f的值趋近于1, 则损失函数趋近于无穷.

因此当 $y^{(i)} = 0$ 时, 损失函数会帮助算法做出更准确的预测, 因为当预测值接近1时, 此时损失非常高.

通过上面两种损失函数的选择, 整体成本函数将是凸函数, 可以可靠地使用梯度下降达到全局最小值.

损失函数

$loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), y^{(i)})$ $$ \begin{equation} loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), y^{(i)}) = \begin{cases} - \log\left(f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) & \text{if $y^{(i)}=1$}\ \log \left( 1 - f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) & \text{if $y^{(i)}=0$} \end{cases} \end{equation} $$ 通用公式 $$ loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), y^{(i)}) = (-y^{(i)} \log\left(f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) - \left( 1 - y^{(i)}\right) \log \left( 1 - f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) $$ 当 $y^{(i)} = 0$ 时 $$ \begin{align} loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), 0) &= (-(0) \log\left(f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) - \left( 1 - 0\right) \log \left( 1 - f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) \ &= -\log \left( 1 - f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) \end{align} $$ 当 $y^{(i)} = 1$ 时 $$ \begin{align} loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), 1) &= (-(1) \log\left(f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) - \left( 1 - 1\right) \log \left( 1 - f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right)\ &= -\log\left(f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) \end{align} $$

成本函数

成本函数

$$ J(\mathbf{w},b) = \frac{1}{m} \sum_{i=0}^{m-1} \left[ loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), y^{(i)}) \right] \tag{1}$$

损失函数

$$loss(f_{\mathbf{w},b}(\mathbf{x}^{(i)}), y^{(i)}) = -y^{(i)} \log\left(f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) - \left( 1 - y^{(i)}\right) \log \left( 1 - f_{\mathbf{w},b}\left( \mathbf{x}^{(i)} \right) \right) \tag{2}$$

其中模型函数f和sigmoid函数 $$ \begin{align} f_{\mathbf{w},b}(\mathbf{x^{(i)}}) &= g(z^{(i)}) \ z^{(i)} &= \mathbf{w} \cdot \mathbf{x}^{(i)}+ b \ g(z^{(i)}) &= \frac{1}{1+e^{-z^{(i)}}} \end{align} $$ 这个特定的成本函数是使用称为最大似然估计(maximum likelihood)的统计原理从统计中推导出来的.

成本函数代码

def compute_cost_logistic(X, y, w, b):
    """
    参数:
    X: 训练集,包含多个特征
    y: 目标值
    w: 模型参数,有多个
    b: 模型参数
    """
    m = X.shape[0]
    cost = 0.0
    
    for i in range(m):
        # 计算函数g(z)的z值
        z_i = np.dot(X[i], w) + b
        # 代入函数g(z)中
        f_wb_i = 1 / (1 + np.exp(-z_i))
        # 计算损失函数
        loss_i = -y[i] * np.log(f_wb_i) - (1 - y[i]) * np(1 - f_wb_i)
        # 所有损失函数之和
        cost = cost + lost_i
    # 最后除以训练集数据数量m
    cost = cost / m
    return cost

实现梯度下降

如上图显示的公式, 是不是和线性回归的梯度下降表达式相似?

但逻辑回归的模型函数f是 $f_{\mathbf{w},b}(\mathbf{x}) = \frac{1}{1+e^{-(\mathbf{w} \cdot \mathbf{x} + b)}}$

而线性回归的模型函数f是 $$f_{\mathbf{w},b}(\mathbf{x}) = \mathbf{w} \cdot \mathbf{x} + b$$

如下图所示:

逻辑回归同样可以使用:

• 监控梯度下降(学习率曲线)
• 矢量化运算
• 特征缩放

梯度下降代码

实现梯度下降中的repeat代码块部分

def compute_gradient_logistic(X, y, w, b):
    m,n = X.shape
    # 注意参数w是特征数量n个,因此使用np的数组初始化
    dj_dw = np.zeros((n,))
    dj_db = 0.0
    
    for i in range(m):
        z_i = np.dot(X[i], w) + b
        # 计算函数f
        f_wb_i = 1 / (1 + np.exp(-z_i))
        # 计算函数f和目标值y的差值
        dj_db_i = f_wb_i - y[i]
        # 计算每个特征对应的参数w的值
        for j in range(n):
            dj_dw[j] = dj_dw[j] + dj_db_i * X[i, j]
        # 计算参数b的值
        dj_db = dj_db + dj_db_i
    # 最后不忘除以训练集的数量m
    dj_dw = dj_dw / m
    dj_db = dj_db / m
    
    return dj_dw, dj_db

实现逻辑回归的梯度下降

def gradient_descent(X, y, w_in, b_in, alpha, num_iters):
    # 存储成本函数J的历史值
    J_history = []
    w = copy.deepcopy(w_in)
    b = b_in
    
    for i in range(num_iters):
        # 通过调用函数获取dj_dw, dj_db
        dj_dw, dj_db = compute_gradient_logistic(X, y, w, b)
        # 梯度下降
        w = w - alpha*dj_dw
        b = b - alpha*dj_db
        
        # 记录成本函数
        if (i < 100000):
            J_history.append(compute_cost_logistic(X, y, w, b))
        # 打印一些统计信息
        if i% math.ceil(num_ites/10) == 0:
            print(f"Iteration {i:4d}: Cost {J_history[-1]}   ")
    return w, b, J_history

使用Scikit-Learn

可使用Scikit-Learn直接封装好的逻辑回归

# 导入库
import numpy as np
from sklearn.linear_model import LogisticRegression

# 训练集
X = np.array([[0.5, 1.5], [1,1], [1.5, 0.5], [3, 0.5], [2, 2], [1, 2.5]])
# 目标值
y = np.array([0, 0, 0, 1, 1, 1])

# 创建逻辑回归模型
lr_model = LogisticRegression()
# 适配模型数据
lr_model.fit(X, y)

# 根据训练集进行预测
y_pred = lr_model.predict(X)
# 打印预测结果
print("Prediction on training set:", y_pred)
# 打印预测的分数-准确率
print("Accuracy on training set:", lr_model.score(X, y))

过拟合问题

描述

对于如上图的模型, 过度的切合训练数据会造成过拟合(overfit)问题.

上图左边的模型对训练数据是欠拟合(underfit)的, 或算法具有高偏差(high bias).

中间的模型对训练数据切合的比较不错, 并且此算法能够很好地进行泛化(generalization), 这意味着即使模型在全新从未见过的数据上也能够做出良好的预测.

右边的模型对训练数据极度的适合, 并且完美通过了所有的训练数据.

• 但当在全新的数据上做出的预测是不正确的(比如图中粉色部分, 房子尺寸大了, 反而价格降低了), 这个模型对数据过度拟合(overfit), 或说这个算法具有高方差(high variance).
• 这个模型不会推广到从未见过的新数据上. 如果训练集有一点点不同, 那么该算法拟合的函数最终可能会完全不同.
• 如果讲模型拟合到略有不同的数据集, 就会得到完全不同的预测或高度可变的预测, 即该算法具有高方差.

对于分类模型, 也有如下的问题:

解决方法

收集更多的训练数据

通过收集更多的训练数据, 学习算法将学会适应一个波动较小的函数.

但在实际中通常是比较困难的.

选择特征

使用太多多项式特征, 也可能会造成过拟合问题.

选择最合适的一组特征来使用, 称为特征选择(feature selection). 关于特征选择, 将在第二课中讲解.

特征选择的一个缺点: 假如所有的特征实际上都可用于预测, 仅使用特征的一个子集, 该算法会丢弃有关的一些信息.

正则化

Regularization

如果要消除某一个特征, 有可能把这个特征的参数设置为0, 正则化是一种更温和地减少某些特征影响的方法, 而不用像彻底消除那样严厉.

正则化的作用是鼓励学习算法缩小参数值, 而不必要求参数正好设置为0.

正则化的作用是可以让你保留所有特征, 但防止特征产生过大的影响, 因为有时会导致过度拟合.

是否对参数b进行正则化并没有太大区别, 推荐不正则化参数b.

应用正则化

原理

一般地说, 正则化的实现方式是, 如果有很多特征, 不知道哪些是最重要的特征, 哪些是要惩罚的特征, 正则化的典型实现方式是惩罚所有的特征, 即惩罚所有的$w_j$参数, 并且有可能表明这通常会导致拟合更平滑, 更简单, 周期更少的函数, 也不太容易过度拟合.

正则化公式

$$J(\mathbf{w},b) = \frac{1}{2m} \sum\limits_{i = 1}^{m} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})^2 + \frac{\lambda}{2m} \sum\limits_{j=1}^{n} w_j^2$$

其中

$$ f_{\mathbf{w},b}(\mathbf{x}^{(i)}) = \mathbf{w} \cdot \mathbf{x}^{(i)} + b \tag{2} $$

参数 $\lambda$ 称为正则化参数

$\frac{1}{2m} \sum\limits_{i = 1}^{m} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})^2$ 部分称为均方误差成本(mean squared error)

$\frac{\lambda}{2m} \sum\limits_{j=1}^{n} w_j^2$ 部分称为正则化项(regularization term)

正则化项和前面的部分都是按2m的比例缩放, 事实证明, 通过相同的方式缩放两项, 为参数 $\lambda$ 选择一个好的值更容易一些.

特别是, 即使你的训练集规模增加, 使用之前相同的参数 $\lambda$ 现在也更有可能继续工作.

对参数b是否进行正则化, 在实践中几乎没有什么区别.

新的成本函数会权衡下面两个目标:

1. 尝试最小化第一项, 会鼓励算法通过最小化预测值和实际值的平方差来很好地拟合训练数据.
2. 尽量减少第二项, 算法试图使参数 $w_j$ 保持较小, 这将倾向于减少过度拟合.

如果参数 $\lambda$ 是个非常大的数, 在正则化项上赋予了非常大的权重. 最小化这种情况的唯一方法是确保参数 $w_j$ 的所有值都非常接近0. 此时模型函数f等于参数b, 成为了一条水平直线, 并且欠拟合.

如果参数 $\lambda$ 为0, 则没有正则化项, 则会拟合为过于摇摆, 过于复杂的曲线, 并且会过度拟合.

因此需要一个合适的参数 $\lambda$, 适当地平衡第一项和第二项, 最小化均方误差并保持参数较小.

正则化线性回归

如上图所示公式, 和前面的基本一致, 只是梯度下降时, 参数 $w_j$ 多了 $\frac{\lambda}{m}w_j$ 项.

通过上图的转换, 可得知正则化在每次迭代所做的是将参数 $w_j$ 乘以一个略小于1的数组, 这会缩小参数 $w_j$ 的值.

公式的推导过程, 如下图

正则化逻辑回归

如上图所示公式, 和前面的基本一致, 只是梯度下降时, 参数 $w_j$ 多了 $\frac{\lambda}{m}w_j$ 项. 和线性回归不同的是模型函数 $f_{\mathbf{w},b}$

代码

正规化线性回归

成本函数

梯度下降

def compute_gradient_linear_reg(X, y, w, b, lambda_):
    m, n = X.shape
    dj_dw = np.zeros((n,))
    dj_db = 0.0
    
    for i in range(m):
        dj_db_i = (np.dot(X[i], w) + b) - y[i]
        for j in range(n):
            dj_dw[j] = dj_dw[j] + X[i, j] * dj_db_i
        dj_db = dj_db + dj_db_i
    dj_dw = dj_dw / m
    dj_db = dj_db / m
    
    for j in range(n):
        dj_dw[j] = dj_dw[j] + (lambda_ / m) * w[j]
    
    return dj_dw, dj_db

正规化逻辑回归

成本函数

def compute_cost_logistic_reg(X, y, w, b, lambda_ = 1):
    m, n = X.shape
    cost = 0.0
    for i in range(m):
        z_i = np.dot(X[i], w) + b
        f_wb_i = 1 / (1 + np.exp(-z_i))
        # 注意损失函数的前面的负号
        cost = cost + -y[i] * np.log(f_wb_i) - (1 - y[i]) * np.log(1 - f_wb_i)
    # 最后除以训练集数据数量m
    cost = cost / m
    
    reg_cost = 0.0
    for j in range(n):
        reg_cost = reg_cost + w[j]2
    reg_cost = (lambda_ / (2 * m)) * reg_cost
    
    total_cost = cost + reg_cost
    return total_cost

def compute_gradient_logistic_reg(X, y, w, b, lambda_):
    m, n = X.shape
    dj_dw = np.zeros((n,))
    dj_db = 0.0
    
    for i in range(m):
        z_i = np.dot(X[i], w) + b
        f_wb_i = 1 / (1 + np.exp(-z_i))
        dj_db_i = f_wb_i - y[i]
        for j in range(n):
            dj_dw[j] = dj_dw[j] + X[i, j] * dj_db_i
        dj_db = dj_db + dj_db_i
    dj_dw = dj_dw / m
    dj_db = dj_db / m
    
    for j in range(n):
        dj_dw[j] = dj_dw[j] + (lambda_ / m) * w[j]
    
    return dj_dw, dj_db