【机器学习笔记】线性回归模型（单特征值与多特征值）

1. 引言

在机器学习的监督学习中最常用的算法就是回归算法，而线性回归又是其中之一

2. 一元线性回归模型（单特征值）

$$f_{w,b}(x^{(i)})=w{x}^{(i)}+b$$

上式即使我们数学中最常见的一元一次方程了，这里的函数对应的就是我们的模型。
我们模型的目标即是通过对数据的训练得到最佳 (w,b) 组合，然后就可以通过该模型来预估任意给定的$x_i$

3. 代价函数

理想情况下，我们所有的数据集可以连接成一条线。但是由于实际数据集庞大且分布无规律，所以无法达成这一点！！
这时候我们就需要通过其他的方法来近似估计模型，使得所有预测数据离原真实数据的差值最小

如下图所示，假设有十个节点的数据集（随机分布)，显然左边的模型更好（均匀分布在两侧）！更适合用来预测

在不同的问题下我们选用代价函数是不同的，但目标不变 — 求出使得整体函数值最小的参数组合，平方损失函数作为最常见的一种，应用还是很广的，即是把预估值与真实值做差的平方再累加后除以样本总数，这里 /m 或者 /2m 都可以，我们选用后者的原因是考虑到后面梯度下降需要求导，而且求导多出来一个2正好可以消掉

3.3 只考虑其中一个变量 (简化)

固定一个变量（假设这里固定b=0），就只剩下一个变量w，其构成的函数图像就是二维图
$$J(w)=\frac{1}{2m}\sum _{i=1}^m(f_w(x^{(i)})-y^{(i)}{)}^2$$

通过不同的w构成的函数（f = w*x）可以得出对应的估计值，将这些结果带入上述方程可以得到不同的函数值，将其连成一条线可以看出是一个抛物线，如下图所示

我们都知道抛物线开口向上有最小值，所以图示中当w等于200时对应的损失函数值最小
画图代码如下

import numpy as np
import matplotlib.pyplot as plt
from lab_utils_uni import plt_intuition, plt_stationary, plt_update_onclick, soup_bowl
plt.style.use('./deeplearning.mplstyle') #老师提供的样式

x_train = np.array([1.0, 2.0])
y_train = np.array([300.0, 500.0])

def compute_cost(x, y, w, b): 
    """
    Computes the cost function for linear regression.
    
    Args:
      x (ndarray (m,)): Data, m examples 
      y (ndarray (m,)): target values
      w,b (scalar)    : model parameters  
    
    Returns
        total_cost (float): The cost of using w,b as the parameters for linear regression
               to fit the data points in x and y
    """
    # number of training examples
    m = x.shape[0] 
    
    cost_sum = 0 
    for i in range(m): 
        f_wb = w * x[i] + b   
        cost = (f_wb - y[i]) ** 2  
        cost_sum = cost_sum + cost  
    total_cost = (1 / (2 * m)) * cost_sum  

    return total_cost

plt_intuition(x_train,y_train)    

3.2 同时考虑变量（w,b）

两个变量构成的元组（w,b），其构成的函数图像就是三维图
$$J(w,b)=\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)}{)}^2$$

三维图往往很难像平面图一样轻松的找到使得损失函数值最小的变量元组，如下图所示，通过plt_stationary可以做出平面等高线图，越往越中心走，值应该越低。但是由于两个变量相互影响所以我们还是很难定位到最小值所对应的（w,b）… 这时候就引出了梯度下降！

画图代码如下

如果出现No module named ‘ipympl’的错误，则需要下载个包先
1.搜索栏输入powershell，以管理员身份打开命令行
2.输入命令 conda install -c conda-forge ipympl

import numpy as np
%matplotlib widget
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D #解决3d包识别不到
from lab_utils_uni import plt_intuition, plt_stationary, plt_update_onclick, soup_bowl
plt.style.use('./deeplearning.mplstyle')

x_train = np.array([1.0, 1.7, 2.0, 2.5, 3.0, 3.2])
y_train = np.array([250, 300, 480,  430,   630, 730,])

plt.close('all') 
fig, ax, dyn_items = plt_stationary(x_train, y_train) #平面等高线图
updater = plt_update_onclick(fig, ax, x_train, y_train, dyn_items)

soup_bowl() #三维立体图

4. 梯度下降算法

用于解决在两个或者两个以上的参数下求得损失函数最优解的方法，在某一时刻同时变换变量的值，使得损失函数更新到一个新的状态。
如下是算法的数学公式，注意该算法更新变量值是两个或者多个同时改变！！！所以需要使用零时变量暂存以防止出现如 用改变后的$w$值带入到第二个公式的情况
$$tmp\_w=w-\alpha \frac{\partial }{{\partial }_w}J(w,b)tmp\_b=b-\alpha \frac{\partial }{{\partial }_b}J(w,b)$$
式中的$\alpha$为学习率（$\alpha$越高，迭代越快，可能导致跳过局部最小值）

代码实现

导包并初始化数据

import math, copy
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('./deeplearning.mplstyle')
from lab_utils_uni import plt_house_x, plt_contour_wgrad, plt_divergence, plt_gradients

# Load our data set
x_train = np.array([1.0, 2.0])   #features
y_train = np.array([300.0, 500.0])   #target value

代价函数

#Function to calculate the cost
def compute_cost(x, y, w, b):
   
    m = x.shape[0] 
    cost = 0
    
    for i in range(m):
        f_wb = w * x[i] + b
        cost = cost + (f_wb - y[i])**2
    total_cost = 1 / (2 * m) * cost

    return total_cost

计算上面讲到的梯度下降数学公式

def compute_gradient(x, y, w, b): 
    """
    Computes the gradient for linear regression 
    Args:
      x (ndarray (m,)): Data, m examples 
      y (ndarray (m,)): target values
      w,b (scalar)    : model parameters  
    Returns
      dj_dw (scalar): The gradient of the cost w.r.t. the parameters w
      dj_db (scalar): The gradient of the cost w.r.t. the parameter b     
     """
    
    # Number of training examples
    m = x.shape[0]    
    dj_dw = 0
    dj_db = 0
    
    for i in range(m):  
        f_wb = w * x[i] + b 
        dj_dw_i = (f_wb - y[i]) * x[i] #（求导，x平方导数为2x）
        dj_db_i = f_wb - y[i] 
        dj_db += dj_db_i #累加（求和部分）
        dj_dw += dj_dw_i 
    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, cost_function, gradient_function): 
    """
    Performs gradient descent to fit w,b. Updates w,b by taking 
    num_iters gradient steps with learning rate alpha
    
    Args:
      x (ndarray (m,))  : Data, m examples 
      y (ndarray (m,))  : target values
      w_in,b_in (scalar): initial values of model parameters  
      alpha (float):     Learning rate
      num_iters (int):   number of iterations to run gradient descent
      cost_function:     function to call to produce cost
      gradient_function: function to call to produce gradient
      
    Returns:
      w (scalar): Updated value of parameter after running gradient descent
      b (scalar): Updated value of parameter after running gradient descent
      J_history (List): History of cost values
      p_history (list): History of parameters [w,b] 
      """
    
    w = copy.deepcopy(w_in) # avoid modifying global w_in
    # An array to store cost J and w's at each iteration primarily for graphing later
    J_history = []
    p_history = []
    b = b_in
    w = w_in
    
    for i in range(num_iters):
        # Calculate the gradient and update the parameters using gradient_function
        dj_dw, dj_db = gradient_function(x, y, w , b)     

        # Update Parameters using equation above
        b = b - alpha * dj_db                            
        w = w - alpha * dj_dw                            

        # Save cost J at each iteration
        if i<100000:      # prevent resource exhaustion 
            J_history.append( cost_function(x, y, w , b))
            p_history.append([w,b])
        # Print cost every at intervals 10 times or as many iterations if < 10
        if i% math.ceil(num_iters/10) == 0:
            print(f"Iteration {i:4}: Cost {J_history[-1]:0.2e} ",
                  f"dj_dw: {dj_dw: 0.3e}, dj_db: {dj_db: 0.3e}  ",
                  f"w: {w: 0.3e}, b:{b: 0.5e}")
 
    return w, b, J_history, p_history #return w and J,w history for graphing

# initialize parameters
w_init = 0
b_init = 0
# some gradient descent settings
iterations = 10000
tmp_alpha = 1.0e-2
# run gradient descent
w_final, b_final, J_hist, p_hist = gradient_descent(x_train ,y_train, w_init, b_init, tmp_alpha, 
                                                    iterations, compute_cost, compute_gradient)
print(f"(w,b) found by gradient descent: ({w_final:8.4f},{b_final:8.4f})")            

得到最终拟合的final值（w,b）后就可以进行预测了！！

5. 多元线性回归模型（多特征值）

前面讲到的都是单一特征值，如房价等。但是实际上往往有多种多样的特征值（比如房贷，房大小，是否装修等等因素）这时候我们就需引入多元线性回归模型。
$$f_{\vec{w},b}(\vec{x})=\vec{w}\cdot \vec{x}+b=w_1{x}_1+w_2{x}_2+w_3{x}_3+...+w_n{x}_n+b$$
上式写成了向量点乘的方式（后面讲原因），这里每一个$x$都对应了一个$w$

5.1 向量化

为什么我们需要引入向量化而不是直接在代码中使用一个forloop遍历呢？
在Python中有NumPy包引入了像我们这里提到的向量等数据类型，方便我们进行运算操作

• 如果使用forloop的话效率大大降低🙃
• 使用包中提供的 np.dot() 方法，该方法可以调用硬件来并行计算（SIMD），然后通过并行加法器将他们都加起来，这在大特征值的条件下能显著提升效率！！

5.2 多元线性回归中的梯度下降

$$w_j=w_j-\alpha \frac{\partial J_{(\vec{w},b)}}{\partial w_j}b=b-\alpha \frac{\partial J_{(\vec{w},b)}}{\partial b}$$上式中j $\in$ (0,n-1)，注意跟之前的一样也需要同时更新！！此时损失函数J求偏导后结果如下所示
$$\frac{\partial J_{(\vec{w},b)}}{\partial w_j}=\frac{1}{m}\sum _{i=0}^{m-1}(f_{\vec{w},b}(x^{(i)})-y^{(i)})x_j^{(i)}\frac{\partial J_{(\vec{w},b)}}{\partial b}=\frac{1}{m}\sum _{i=0}^{m-1}(f_{\vec{w},b}(x^{(i)})-y^{(i)})$$ (根据上面提到的$f_{\vec{w},b}(\vec{x})$可以得出对$w_j{x}_j$求导就是其他看作常数，只剩下$x_j$)

6. 调参

上面公式这么多设定参数，我们应该怎么确定呢？

• 像是$\alpha$不能设置太高也不能设置太低，老师建议从0.001开始每次*3直到1慢慢根据图像确定$\alpha$选的是否合理，应该往大或者往小调整…
• 由于不同特征值的计量单位往往不同，在梯度下降的过程中会导致时间消耗增加，这时候我们就可以通过归一化来使所有特征值归转化为标准正态分布，提升数据的可比性…有多种方法可以实现这里介绍z-score标准化，即是每个特征值减去期望除以方差
  $$x_j^{(i)}=\frac{x_j^{(i)}-{\mu }_j}{{\sigma }_j}其中{\mu }_j=\frac{1}{m}\sum _{i=0}^{m-1}{x}_j^{(i)},{\sigma }_j^2=\frac{1}{m}\sum _{i=0}^{m-1}(x_j^{(i)}-{\mu }_j{)}^2$$ 下图例子是经过算法处理后的数据，可以清楚看出服从正态分布
  z-score代码如下

def zscore_normalize_features(X,rtn_ms=False):
  """
  returns z-score normalized X by column
  Args:
    X : (numpy array (m,n)) 
  Returns
    X_norm: (numpy array (m,n)) input normalized by column
  """
  mu     = np.mean(X,axis=0)  
  sigma  = np.std(X,axis=0)
  X_norm = (X - mu)/sigma      

  if rtn_ms:
      return(X_norm, mu, sigma)
  else:
      return(X_norm)

值得注意的是，当我们拟合出 $\vec{w}$ 和 $b$ 后就可以进行估计未知的数据了。但是对于位置的数据我们同样也需要对其经行归一化处理来保证一致性！！！