一元線性回歸的Python實現

来源:https://www.cnblogs.com/hznudmh/archive/2022/06/17/16382786.html
-Advertisement-
Play Games

1 問題的提出 對於給定的數據集 \(D = \{(x_1,y_1),(x_2,y_2),\cdots,(x_m,y_m)\}\),線性回歸 (linear regression) 試圖學得一個線性模型以儘可能準確地預測是指輸出標記. 2 原理 設給定的數據集 \(D = \{(x_i,y_i)\} ...


目錄

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


您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • 這裡給大家分享我在OpenLayers 地圖開發工作中總結出的一下代碼和註意點,希望對大家有所幫助 效果如下: 核心代碼展示:附帶講解註釋 var map = new ol.Map({ // 初始化地圖 target: 'map',// 選擇地圖對象 layers: [ new ol.layer.T ...
  • springboot整合ueditor實現圖片上傳和文件上傳功能 寫在前面: 在閱讀本篇之前,請先按照我的這篇隨筆完成對ueditor的前期配置工作: springboot+layui 整合百度富文本編輯器ueditor入門使用教程(踩過的坑)https://www.cnblogs.com/rain ...
  • 1、前言 作為嵌入式軟體開發,可能經常會使用命令行或者顯示屏等設備實現人機交互的功能,功能中通常情況都包含 UI 菜單設計;很多開發人員都會有自己的菜單框架模塊,防止重覆造輪子,網上有很多這種菜單框架的代碼,但是大多耦合性太強,無法獨立出來適配不同的菜單設計。 本文介紹一個降低了耦合性,完全獨立的菜 ...
  • 目錄 一.簡介 二.效果演示 三.源碼下載 四.猜你喜歡 零基礎 OpenGL (ES) 學習路線推薦 : OpenGL (ES) 學習目錄 >> OpenGL ES 基礎 零基礎 OpenGL (ES) 學習路線推薦 : OpenGL (ES) 學習目錄 >> OpenGL ES 轉場 零基礎 O ...
  • 分類方式 按參數分: 有參構造(預設構造) & 無參構造 按類型分: 普通構造 & 拷貝構造 調用方式 括弧法 顯示法 隱式轉換法 PS:下方所有文本均以此代碼為基礎 1 class Person { 2 public: 3 //無參構造函數 4 Person() { 5 std::cout << ...
  • 引言:今天工作遇到了一個需要按行讀取txt文件數據的需求,查詢了一下getline()函數,發現這竟然是一個C++的標準庫函數,而且設計的很好,特地做一下記錄。getline本質是一個定界流輸入截取函數,預設是換行符‘/n’ 個人技術博客(文章整理+源碼): https://zobolblog.gi ...
  • 大佬的理解->《IO流和File》 1、File類 File類是IO包中唯一代表磁碟文件本身的對象,File類定義了一些與平臺無關的方法來操作文件。通過調用File類提供的各種方法,能夠完成創建、刪除文件、重命名文件、判斷文件的讀寫許可權許可權是否存在、設置和查詢文件的最近修改時間等操作。 ​ File ...
  • 前言 大部分情況下,地理繪圖可使用 Arcgis 等工具實現。但正版的 Arcgis 並非所有人可以承受。本文基於 Python 的 cartopy 和 matplotlib 等庫,為地理空間繪圖的代碼實現提供參考。 所有所需庫如下: gma、cartopy、matplotlib、numpy 更多內 ...
一周排行
    -Advertisement-
    Play Games
  • 前言 當別人做大數據用Java、Python的時候,我使用.NET做大數據、數據挖掘,這確實是值得一說的事。 寫的並不全面,但都是實際工作中的內容。 .NET在大數據項目中,可以做什麼? 寫腳本(使用控制台程式+頂級語句) 寫工具(使用Winform) 寫介面、寫服務 使用C#寫代碼的優點是什麼? ...
  • 前言 本文寫給想學C#的朋友,目的是以儘快的速度入門 C#好學嗎? 對於這個問題,我以前的回答是:好學!但仔細想想,不是這麼回事,對於新手來說,C#沒有那麼好學。 反而學Java還要容易一些,學Java Web就行了,就是SpringBoot那一套。 但是C#方向比較多,你是學控制台程式、WebAP ...
  • 某一日晚上上線,測試同學在回歸項目黃金流程時,有一個工單項目介面報JSF序列化錯誤,馬上升級對應的client包版本,編譯部署後錯誤消失。 線上問題是解決了,但是作為程式員要瞭解問題發生的原因和本質。但這都是為什麼呢? ...
  • 本文介紹基於Python語言中TensorFlow的Keras介面,實現深度神經網路回歸的方法。 1 寫在前面 前期一篇文章Python TensorFlow深度學習回歸代碼:DNNRegressor詳細介紹了基於TensorFlow tf.estimator介面的深度學習網路;而在TensorFl ...
  • 前段時間因業務需要完成了一個工作流組件的編碼工作。藉著這個機會跟大家分享一下整個創作過程,希望大家喜歡,組件暫且命名為"easyFlowable"。 接下來的文章我將從什麼是工作流、為什麼要自研這個工作流組件、架構設計三個維度跟大家來做個整體介紹。 ...
  • 1 簡介 我們之前使用了dapr的本地托管模式,但在生產中我們一般使用Kubernetes托管,本文介紹如何在GKE(GCP Kubernetes)安裝dapr。 相關文章: dapr本地托管的服務調用體驗與Java SDK的Spring Boot整合 dapr入門與本地托管模式嘗試 2 安裝GKE ...
  • 摘要:在jvm中有很多的參數可以進行設置,這樣可以讓jvm在各種環境中都能夠高效的運行。絕大部分的參數保持預設即可。 本文分享自華為雲社區《為什麼需要對jvm進行優化,jvm運行參數之標準參數》,作者:共飲一杯無。 我們為什麼要對jvm做優化? 在本地開發環境中我們很少會遇到需要對jvm進行優化的需 ...
  • 背景 我們的業務共使用11台(阿裡雲)伺服器,使用SpringcloudAlibaba構建微服務集群,共計60個微服務,全部註冊在同一個Nacos集群 流量轉發路徑: nginx->spring-gateway->業務微服務 使用的版本如下: spring-boot.version:2.2.5.RE ...
  • 基於php+webuploader的大文件分片上傳,帶進度條,支持斷點續傳(刷新、關閉頁面、重新上傳、網路中斷等情況)。文件上傳前先檢測該文件是否已上傳,如果已上傳提示“文件已存在”,如果未上傳則直接上傳。視頻上傳時會根據設定的參數(分片大小、分片數量)進行上傳,上傳過程中會在目標文件夾中生成一個臨 ...
  • 基於php大文件分片上傳至七牛雲,使用的是七牛雲js-sdk V2版本,引入js文件,配置簡單,可以暫停,暫停後支持斷點續傳(刷新、關閉頁面、重新上傳、網路中斷等情況),可以配置分片大小和分片數量,官方文檔https://developer.qiniu.com/kodo/6889/javascrip ...