[机器学习-从入门到入土] 线性回归

个人导航

知乎:https://www.zhihu.com/people/byzh_rc

CSDN:https://blog.csdn.net/qq_54636039

linear-regression.ipynb

假设变量之间的关系是线性的
h θ ( x ) = θ 0 + θ 1 x h_{\theta}(x)= \theta_{0} + \theta_{1} x hθ(x)=θ0+θ1x

θ \pmb{\theta} θ 就是学习算法需要学习的参数,在线性回归的问题上,就是 θ 1 \theta_{1} θ1 θ 0 \theta_{0} θ0
x x x 是我们对于问题所选取的特征,也即输入
h h h表示算法得到的映射

根据“真实值与算法拟合值的误差”来表示算法的“合适程度”。
在线性回归中,我们经常使用最小二乘的思路构建代价函数 J ( θ ) J(\theta) J(θ)
J ( θ ) = 1 2 m ∑ i = 1 m ( h θ ( x ( i ) ) − y ( i ) ) 2 = 1 2 m ∑ i = 1 m ( ( θ 0 + θ 1 x ( i ) ) − y ( i ) ) 2 J(\pmb{\theta}) = \frac{1}{2m}\sum_{i=1}^{m} \Big( h_{\theta}(x^{(i)}) - y^{(i)} \Big)^2=\frac{1}{2m}\sum_{i=1}^{m} \Big( (\theta_{0} + \theta_{1} x^{(i)}) - y^{(i)} \Big)^2 J(θ)=2m1i=1m(hθ(x(i))y(i))2=2m1i=1m((θ0+θ1x(i))y(i))2

误差函数/代价函数(cost function) 的值越小,则代表算法拟合结果与真实结果越接近

J = np.sum(((X @ theta).flatten() - y) ** 2) / (2 * m)
# input_size=2 (bias和weight)
# X: (m, input_size) y: (m,)
# theta: (input_size, 1)
# -> J: 数

梯度下降算法沿着误差函数的反向来更新 θ \theta θ的值,直到代价函数收敛到最小值
θ j = θ j − α ∂ ∂ θ j J ( θ ) , j = 0 , 1 ∂ J ∂ θ 0 = 1 m ∑ i = 1 m ( h θ ( x ( i ) ) − y ( i ) ) , ∂ J ∂ θ 1 = 1 m ∑ i = 1 m ( h θ ( x ( i ) ) − y ( i ) ) x ( i ) . \theta_{j} = \theta_{j} - \alpha\frac{\partial }{\partial \theta_{j}}J(\pmb{\theta}), \qquad j=0, 1 \\ \begin{split}\frac{\partial J}{\partial \theta_{0} } &= \frac{1}{m}\sum_{i=1}^{m} \Big( h_{\theta}(x^{(i)}) - y^{(i)} \Big),\\ \frac{\partial J}{\partial \theta_{1} } &= \frac{1}{m}\sum_{i=1}^{m} \Big( h_{\theta}(x^{(i)}) - y^{(i)} \Big)x^{(i)}.\end{split} θj=θjαθjJ(θ),j=0,1θ0Jθ1J=m1i=1m(hθ(x(i))y(i)),=m1i=1m(hθ(x(i))y(i))x(i).

α \alpha α表示学习率

theta = theta - alpha / m * (((X @ theta).reshape(-1) - y) @ X).reshape(-1, 1)
# input_size=2 (bias和weight)
# X: (m, input_size) y: (m,)
# theta: (input_size, 1)

python代码

import numpy as np
import matplotlib.pyplot as plt

######################### cell 1 #########################
# 读取数据
def plot_data1(x, y):
    """绘制给定数据x与y的图像"""
    plt.figure()
    # ====================== 你的代码 ==========================
    plt.plot(x, y, 'rx', markersize=8, markeredgewidth=2, label='scatter plot')
    plt.xlabel('Population of City in 10,000s')
    plt.ylabel('Profit in $10,000s')
    plt.legend(loc='lower right')
    # =========================================================
    
# 从txt中加载数据
print('Plotting Data ...\n')
data = np.loadtxt('./PRML_LR_data.txt', delimiter=',')
x, y = data[:, 0], data[:, 1]
# print(x.shape) # (97,)
# print(y.shape) # (97,)

# 绘图
plot_data1(x, y)
plt.show()

######################### cell 2 #########################
# 数据预处理
m = len(y)
X = np.ones((m, 2))
X[:, 1] = data[:, 0] # X的第一列是偏置项的系数, 均为1

# 纯着theta0和theta1
theta = np.zeros((2, 1))

iterations = 1500 # 迭代次数
alpha = 0.01 # 学习率

######################### cell 3 #########################
# 代价函数内(cost fucntion)
def compute_cost(X, y, theta):
    """计算线性回归的代价。"""
    m = len(y)
    # ====================== 你的代码 ==========================
    # 计算给定 theta 参数下线性回归的代价
    J = np.sum(((X@theta).flatten() - y)**2) / (2*m)
    # =========================================================
    return J

J0 = compute_cost(X, y, theta)
print(J0) # 32.07

######################### cell 4 #########################
# 梯度下降(gradient descent)
def gradient_descent(X, y, theta, alpha, num_iters):
    """执行梯度下降算法来学习参数 theta。"""
    m = len(y)
    J_history = np.zeros((num_iters,))

    for iter in range(num_iters):
        # ====================== 你的代码 ==========================
        # 计算给定 theta 参数下线性回归的梯度,实现梯度下降算法
        theta = theta - alpha/m*(((X@theta).reshape(-1) - y)@X).reshape(-1, 1)
        # =========================================================
        # 将各次迭代后的代价进行记录
        J_history[iter] = compute_cost(X, y, theta)

    return theta, J_history

theta, J_history = gradient_descent(X, y, theta, alpha, iterations)
print(theta) # [-3.630291, 1.166362]

######################### cell 5 #########################
# 绘制代价函数(等高线图)
def plot_visualize_cost(X, y, theta_best):
    """可视化代价函数"""

    # 生成参数网格
    theta0_vals = np.linspace(-10, 10, 101)
    theta1_vals = np.linspace(-1, 4, 101)
    t = np.zeros((2, 1))
    J_vals = np.zeros((101, 101))
    for i in range(101):
        for j in range(101):
            # =============== 你的代码 ===================
            # 加入代码,计算 J_vals 的值
            t[0, 0] = theta0_vals[i]
            t[1, 0] = theta1_vals[j]
            J_vals[i, j] = compute_cost(X, y, t)
            # ===========================================

    plt.figure()
    plt.contour(theta0_vals, theta1_vals, J_vals,
                levels=np.logspace(-2, 3, 21)) # 等高线图
    plt.plot(theta_best[0], theta_best[1], 'rx',
             markersize=8, markeredgewidth=2)
    plt.xlabel(r'$\theta_0$')
    plt.ylabel(r'$\theta_1$')
    plt.title(r'$J(\theta)$')
  

plot_visualize_cost(X, y, theta)
plt.show()

######################### cell 6 #########################

def plot_visual_history(X, y):
    plt.figure()
    # =============== 你的代码 ===================
    # 生成参数网格
    theta0_vals = np.linspace(-10, 10, 101)
    theta1_vals = np.linspace(-1, 4, 101)
    t = np.zeros((2, 1))
    J_vals = np.zeros((101, 101))
    for i in range(101):
        for j in range(101):
            t[0, 0] = theta0_vals[i]
            t[1, 0] = theta1_vals[j]
            J_vals[i, j] = compute_cost(X, y, t)
    plt.contour(theta0_vals, theta1_vals, J_vals,
                levels=np.logspace(-2, 3, 21))
    
    num_iters, alpha = 1500, 0.01
    theta = np.zeros((2, 1))
    m = len(y)
    J_history = np.zeros((num_iters,))
    for iter in range(num_iters):
        if iter%100 == 0:
            plt.plot(theta[0], theta[1], 'rx', markersize=8, markeredgewidth=2)
        theta = theta - alpha/m*(((X@theta).reshape(-1) - y)@X).reshape(-1, 1)
        J_history[iter] = compute_cost(X, y, theta)
    # ===========================================
    plt.xlabel(r'$\theta_0$')
    plt.ylabel(r'$\theta_1$')
    plt.title(r'$J(\theta)$')
    

plot_visual_history(X, y)
plt.show()

######################### cell 7 #########################
# 拟合曲线
def plot_data2(x, y):
    """绘制给定数据x与y的图像"""
    plt.figure()
    # ====================== 你的代码 ==========================
    plt.plot(x, y, 'rx', markersize=8, markeredgewidth=2, label='scatter plot')
    plt.plot(
        [x.min(), x.max()], [theta[0] + x.min() * theta[1], theta[0] + x.max() * theta[1]],
        markersize=8, markeredgewidth=2, label='linear regression'
    )  # 取两个点来画图, [x1, x2], [y1, y2]
    plt.xlabel('Population of City in 10,000s')
    plt.ylabel('Profit in $10,000s')
    plt.legend(loc='lower right')
    # =========================================================

plot_data2(x, y)
plt.show()

PRML_LR_data.txt

shape=(97, 2)
X: (97,)
Y: (97,)

6.1101,17.592
5.5277,9.1302
8.5186,13.662
7.0032,11.854
5.8598,6.8233
8.3829,11.886
7.4764,4.3483
8.5781,12
6.4862,6.5987
5.0546,3.8166
5.7107,3.2522
14.164,15.505
5.734,3.1551
8.4084,7.2258
5.6407,0.71618
5.3794,3.5129
6.3654,5.3048
5.1301,0.56077
6.4296,3.6518
7.0708,5.3893
6.1891,3.1386
20.27,21.767
5.4901,4.263
6.3261,5.1875
5.5649,3.0825
18.945,22.638
12.828,13.501
10.957,7.0467
13.176,14.692
22.203,24.147
5.2524,-1.22
6.5894,5.9966
9.2482,12.134
5.8918,1.8495
8.2111,6.5426
7.9334,4.5623
8.0959,4.1164
5.6063,3.3928
12.836,10.117
6.3534,5.4974
5.4069,0.55657
6.8825,3.9115
11.708,5.3854
5.7737,2.4406
7.8247,6.7318
7.0931,1.0463
5.0702,5.1337
5.8014,1.844
11.7,8.0043
5.5416,1.0179
7.5402,6.7504
5.3077,1.8396
7.4239,4.2885
7.6031,4.9981
6.3328,1.4233
6.3589,-1.4211
6.2742,2.4756
5.6397,4.6042
9.3102,3.9624
9.4536,5.4141
8.8254,5.1694
5.1793,-0.74279
21.279,17.929
14.908,12.054
18.959,17.054
7.2182,4.8852
8.2951,5.7442
10.236,7.7754
5.4994,1.0173
20.341,20.992
10.136,6.6799
7.3345,4.0259
6.0062,1.2784
7.2259,3.3411
5.0269,-2.6807
6.5479,0.29678
7.5386,3.8845
5.0365,5.7014
10.274,6.7526
5.1077,2.0576
5.7292,0.47953
5.1884,0.20421
6.3557,0.67861
9.7687,7.5435
6.5159,5.3436
8.5172,4.2415
9.1802,6.7981
6.002,0.92695
5.5204,0.152
5.0594,2.8214
5.7077,1.8451
7.6366,4.2959
5.8707,7.2029
5.3054,1.9869
8.2934,0.14454
13.394,9.0551
5.4369,0.61705

Logo

码道开发者社区,聚焦华为云码道 CodeArts 代码智能体,沉淀 Agent、Skill、鸿蒙开发实战内容,供开发者查阅资料、交流技术、分享工程实践

更多推荐