ICode9

精准搜索请尝试: 精确搜索
首页 > 编程语言> 文章详细

一元线性回归的Python实现

2022-06-17 17:36:32  阅读:199  来源: 互联网

标签:一元 plt frac Python sum np 线性 partial wx


目录

1 问题的提出

对于给定的数据集 \(D = \{(x_1,y_1),(x_2,y_2),\cdots,(x_m,y_m)\}\),线性回归 (linear regression) 试图学得一个线性模型以尽可能准确地预测是指输出标记.

2 原理

设给定的数据集 \(D = \{(x_i,y_i)\}_{i=1}^m, \ x_i,y_i \in \mathcal{R}\). 对于离散属性,如果属性值间存在“序”(order)的关系,可通过连续化将其转化为连续值,例如二值属性“身高”的取值“高”“矮”可转化为 \(\{1.0,0.0\}\),三值属性“高度”的取值“高”“中”“低”可转化为 \(\{1.0,0.5,0.0\}\);若属性之间不存在有序关系,假定有 \(k\) 个属性值,则通常转化为 \(k\) 维向量,例如属性“瓜类”的取值“西瓜”“南瓜”“黄瓜”可转化为 \((0,0,1),(0,1,0),(1,0,0)\).

线性回归试图学得

\[f(x_i) = wx_i + b, \ s.t. f(x_i) \simeq y_i \tag{1} \]

2.1 代价函数

我们需要确定 \(w\) 和 \(b\) 的值使得 \(f(x)\) 和 \(y\) 之间的差别尽可能小. 因此我们引入均方误差,均方误差是回归问题中最常用的性能度量工具,我们可以试图让均方误差最小化,即

\[\begin{aligned} (w^*,b^*) &= \arg\min\limits_{(w,b)} \sum_{i=1}^m(f(x_i)-y_i)^2\\ &=\arg\min\limits_{(w,b)} \sum_{i=1}^m(y_i-wx_i-b)^2 \end{aligned} \tag{2} \]

\(\arg\min\limits_{\theta}f(x;\theta)\) 意思就是找出一个 \(\theta\) 使得 \(f(x;\theta)\) 的值最小,即他的返回值是 \(f(x;\theta)\) 的最小值所对应的 \(\theta\) 的值.

均方误差有很好的几何意义,它对应了常用的“欧式距离”(Euclidean distance),基于均方误差最小化的进行模型求解的方法称为“最小二乘法”(least square method). 在线性回归中,最小二乘法就是试图寻找一条直线,使得所有样本点到直线的欧氏距离之和最小.

在均方误差的基础上进一步构造代价函数

\[J(w,b) = \frac{1}{2m}\sum_{i=1}^m(f(x_i)-y_i)^2 = \frac{1}{2m}\sum_{i=1}^m(wx_i + b - y_i)^2 \]

这里分母的 \(2\) 是为了后续求导的方便

求解 \(w\) 和 \(b\) 使 \(J(w,b)\) 最小化的过程,称为线性回归模型的最小二乘“参数估计”(parameter estimation). 我们可将 \(J(w,b)\) 分别对 \(w\) 和 \(b\) 求导,得到

\[\begin{aligned} &\frac{\partial J(w,b)}{\partial w} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)x_i\\ &\frac{\partial J(w,b)}{\partial b} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i) \end{aligned} \]

令上两式等于0,解得

\[w = \frac{\sum_{i=1}^my_i(x_i-\bar{x})}{\sum_{i=1}^mx_i^2-\frac{1}{m}(\sum_{i=1}^mx_i)^2}\\ b = \frac{1}{m}\sum_{i=1}^m(y_i-wx_i) \]

\(J(w,b)\) 是关于 \(w,b\) 的凸函数,根据凸函数的性质,其偏导为 \(0\) 时就是 \(w\) 和 \(b\) 的最优解.

其中 \(\bar{x} = \frac{1}{m}\sum_\limits{i=1}^mx_i\) 为 \(x\) 的均值.

2.2 模型的评价

2.2.1 皮尔逊相关系数

使用相关系数衡量线性相关性的强弱,皮尔逊相关系数的公式如下:

\[r_{xy} = \frac{COV(X,Y)}{\sqrt{Var(X)Var(Y)}} \]

相关度越高,皮尔逊相关系数的值就趋于 1 或 -1 (趋于 1 表示它们呈正相关, 趋于 -1 表示它们呈负相关);如果相关系数等于0,表明它们之间不存在线性相关关系.

2.2.2 决定系数

决定系数 \(R^2\) 也称拟合优度,反应了 \(y\) 的波动有多少百分比能被 \(x\) 的波动所描述. 决定系数越接近 1 ,说明拟合程度越好.

总平方和

\[SST = \sum_{i=1}^n(y_i-\bar{y})^2 \]

回归平方和

\[SSR = \sum_{i=1}^n(\hat{y}_i-\bar{y})^2 \]

残差平方和

\[SSE = \sum_{i=1}^n(y_i-\hat{y}_i)^2 \]

其中

\[SST = SSR + SSE\\ R^2 = \frac{SSR}{SST} = 1- \frac{SSE}{SST} \]

3 Python 实现

3.1 不调sklearn库

Step1. 调库

import numpy as np
from numpy import genfromtxt
import matplotlib.pyplot as plt

Step2. 数据导入并绘制散点图

data = genfromtxt("data.csv", delimiter = ",")
x = data[:, 0, np.newaxis]
y = data[:, 1, np.newaxis]
plt.scatter(x, y)
plt.show()

Step3. 求回归系数
根据先前的推导,已经知道

\[w = \frac{\sum_{i=1}^my_i(x_i-\bar{x})}{\sum_{i=1}^mx_i^2-\frac{1}{m}(\sum_{i=1}^mx_i)^2}\\ b = \frac{1}{m}\sum_{i=1}^m(y_i-wx_i) \]

m = len(x)
x_bar = np.mean(x)
w = (np.sum((x - x_bar)*y))/(np.sum(x**2)-(1/m)*(np.sum(x))**2)
b = np.mean(y-w*x)
print(w,b)

1.3224310227553517 7.991020982270779

Step4. 拟合图像

plt.plot(x, y, 'b.')
plt.plot(x, w*x+b, 'r')
plt.show()

Step5. 计算相关系数和决定系数

COVxy = np.cov(x.T,y.T)
r = COVxy[0,1]/(x.std()*y.std())
print(r)

0.78154393928063

y_hat = w*x+b
SSR = np.sum((y_hat-np.mean(y))**2)
SST = np.sum((y-np.mean(y))**2)
R2 = SSR/SST
print(R2)

0.5986557915386548

3.2 调 sklearn 库

建模:

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x, y)

LinearRegression()

拟合图像的得出

plt.plot(x, y, 'b.')
plt.plot(x, model.predict(x), 'r')
plt.show()

回归系数

w = model.intercept_
b = model.coef_
print("截距为 {0}, 回归系数为 {1}".format(w[0], b[0][0]))

截距为 7.991020982270385, 回归系数为 1.32243102275536

决定系数

model.score(x, y)

0.598655791538662

4 梯度下降法

4.1 原理

由于代价函数是凸函数,因此只有全局最小值,梯度下降法的原理是先定一个初始值,然后利用导数就像阶梯一样慢慢逼近全局最小值

图片出处:https://zhuanlan.zhihu.com/p/36564434

已知代价函数 \(J(w,b)\),我们需要找一组 \(w,b\) 使得 \(J(w,b)\) 最小,给定一个算法:

\[\begin{aligned} &给定 \ w,b \ 的初始值 \ w_0,b_0\\ \\ &repeat until convergence \ \{\\ &\quad\quad temp0 = w - \alpha \frac{\partial}{\partial w}J(w,b)\\ &\quad\quad temp1 = b - \alpha \frac{\partial}{\partial b}J(w,b)\\ &\quad\quad w = temp0\\ &\} \end{aligned} \]

其中 \(\alpha\) 为学习率,学习率不能太大也不能太小,可以多次尝试 \(0.1,0.03,0.01,0.003,0.001,0.0003,\cdots\).

根据已知条件,有

\[J(w,b) = \frac{1}{2m}\sum_{i=1}^m(wx_i + b - y_i)^2\\ \frac{\partial J(w,b)}{\partial w} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)x_i\\ \frac{\partial J(w,b)}{\partial b} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i) \]

于是

\[w = w - \alpha \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)\\ b = b - \alpha \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)x_i \]

4.2 Python实现

# 学习率learning rate
lr = 0.0001
# 截距初值
b = 0
# 斜率初值
w = 0
# 最大迭代次数
epochs = 50
# 最小二乘法
def compute_error(b, w, x, y):
  totalError = 0
  for i in range(0, len(x)):
    totalError += (y[i] - (w * x[i] + b)) ** 2
  return totalError / float(len(x)) / 2.0

def gradient_descent_runner(x, y, b, w, lr, epochs):
  # 计算总数据量
  m = float(len(x))
  # 循环epochs次
  for i in range(epochs):
    b_grad = 0
    w_grad = 0
    # 计算梯度的总和再求平均
    for j in range(0, len(x)):
      b_grad += -(1/m) * (y[j] - ((w * x[j]) + b))
      w_grad += -(1/m) * x[j] * (y[j] - ((w * x[j]) + b))
     
    # 更新 b 和 w
    b = b - (lr * b_grad)
    w = w - (lr * w_grad)
    
    # 每迭代5次,输出一次图像
    if i % 5 == 0:
      print('epochs:', i)
      plt.plot(x, y, 'b.')
      plt.plot(x, w*x + b, 'r')
      plt.show()
  return b,w
print('Starting b = {0}, w = {1}, error = {2}'.format(b, w, compute_error(b, w, x, y)))
print('Running')
b, w = gradient_descent_runner(x, y, b, w, lr, epochs)
print('After {0} iterations b = {1}, w = {2}, error = {3}'.format(epochs, b, w, compute_error(b,w,x,y)))
# 画图
plt.plot(x, y, 'b.')
plt.plot(x, w * x + b, 'r')
plt.show()

迭代50次后得到 \(b = 0.03207192, w = 1.47886174\),最小二乘误差为 56.3244305.

参考

[1]. 周志华.《机器学习》.清华大学出版社
[2]. https://www.bilibili.com/video/BV1Rt411q7WJ?p=4&vd_source=08d2535b05740d396ec0fc720c52f36a

标签:一元,plt,frac,Python,sum,np,线性,partial,wx
来源: https://www.cnblogs.com/hznudmh/p/16382786.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有