當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

最小二乘法拟合直线（问题i）

發布時間：2024/1/1 编程问答 25 豆豆

生活随笔收集整理的這篇文章主要介紹了最小二乘法拟合直线（问题i）小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

問題描述

有一個函數，現在不知道函數的的具體形式，給定滿足函數關系的一組訓練樣本，請使用線性回歸模型擬合出函數y=f(x)。
(可嘗試一種或幾種不同的基函數，如多項式、高斯或sigmoid基函數）

要求：

i. 先完成最小二乘法的優化 (參考書中第二章 2.3節中的公式)

ii. 附加題：實現“多項式基函數”以及“高斯基函數”（可參考PRML）

iii. 附加題：完成梯度下降的優化 (參考書中第二章 2.3節中的公式)
使用訓練集train.txt 進行訓練，使用測試集test.txt 進行評估（標準差），訓練模型時請不要使用測試集。

讀取數據

我們先將采集到的數據個讀出來：這里以txt的格式保存，利用numpy庫從文本中讀取數據（注意jupyter notebook中數據必須保存到所建立）：

import numpy as np #導入numpy庫 import matplotlib.pyplot as plt train_dataset = np.loadtxt('train.txt') test_dataset = np.loadtxt('test.txt') #print(len(train_dataset)) #print(len(test_dataset)) #print(train_dataset) #print(test_dataset)

在讀出所給的數據之后，我們按照題目所給的要求完成最小二乘法的優化（參照書上2.3節公式）這里單單我們看書是很難理解的，因為合理給出了優化參數的方法，但是我拿到這道題目的時候被難住的是最開始的基本線性回歸的知識。題目中提到了用到基函數（內心很是mengbi）下面將一些基礎知識整理一下，然后自己再試著編程看一看。

線性回歸（線性基函數）

最簡單的線性回歸模型的形式如下所示：
$y(x,w)=w0+w1x+???+wnxny\left( \mathbf{x,w} \right) =w_0+w_1x+···+w_nx_n$
其中 $x=(x1,x2,?,xn)T\mathbf{x}=\left( x_1,x_2,\cdots ,x_n \right) ^T$ , $y$ 是輸入變量 $x_1$ 的線性組合，這種模型具有很大的局限性。
一般用向量形式寫成： $f(x)=wTx+bf\left( x \right) =\mathbf{w}^Tx+b$
其中$\mathbf{w}=\left( w_1,w_2,\cdots ,w_n \right) $,$ w $和$ b$學得之后，模型就可以確定。

我們可以將其擴展為輸入變量通過有限個非線性函數轉換后的線性組合，如下式表示： $y(x,w)=w0+∑j=1n?1wj?j(x)y\left( \mathbf{x,w} \right) =w_0+\sum_{j=1}^{n-1}{w_j\phi _j\left( \mathbf{x} \right)}$
這里$\phi \left( \mathbf{x} \right) $就是基函數。我們可以添加一個基函數$ \phi 0=1 $，則上式就可以改為$ $y\left( \mathbf{x,w} \right) =\sum{j=0}^{M-1}{w_j\phi _j\left( \mathbf{x} \right)}=\mathbf{w}^T\phi \left( \mathbf{x} \right) $$
這里 $w=(w0,w1,?wn?1)T\mathbf{w}=\left( w_0,w_1,\cdots w_{n-1} \right) ^T$
, $?=(?0,?1,??n?1)T\phi =\left( \phi _0,\phi _1,\cdots \phi _{n-1} \right) ^T$ 基函數$\phi _j\left( x \right) $可以有多種不同類型，常見的基函數有：

·多項式基函數： $?j(x)=xj\phi _j\left( x \right) =x^j$

·高斯基函數： $?j(x)=e{?(x?uj)22s2}\phi _j\left( x \right) =e^{\left\{ -\frac{\left( x-u_j \right) ^2}{2s^2} \right\}}$
這里 $u$ 控制位置 $s$ 控制跨度

· $s i g m o i d$ 基函數 $?j(x)=σ(x?ujs)\phi _j\left( x \right) =\sigma \left( \frac{x-u_j}{s} \right)$
$σ(a)=11+e(?a)\sigma \left( a \right) =\frac{1}{1+e^{\left( -a \right)}}$

參數估計

一般的情況是數據集 $D_i$ 樣本是由 $d$ 個屬性描述，此時我們試圖學得 $f(xi)=wTxi+bf\left( x_i \right) =\mathbf{w}^Tx_i+b$ 使得 $f(xi)?yi\ f\left( x_i \right) \backsimeq y_i$ ,這稱為“多元線性回歸”(multivariate linear reegression)

類似的我們可以用最小二乘法對 $w$ 和 $b$ 進行估計，為了便于討論我們將 $w$ 和 $b$ 吸入向量形式$\mathbf{\hat{w}}=\left( \mathbf{w:b} \right) $，相應的，把數據集$ D $表示為一個$ m×(d+1) $大小的矩陣$ X$,其中每行對應于一個示例表示為下述矩陣形式：
$X=[lx11x12?x1d1x21x22?x2d1?????xm1xm2?xmd1]=[lx1T1x2T1??xmT1]\mathbf{X}=\left[ \begin{matrix}{l} x_{11}& x_{12}& \cdots& x_{1d}& 1\\ x_{21}& x_{22}& \cdots& x_{2d}& 1\\ \vdots& \vdots& \ddots& \vdots& \vdots\\ x_{m1}& x_{m2}& \cdots& x_{md}& 1\\ \end{matrix} \right] =\left[ \begin{matrix}{l} x_{1}^{T}& 1\\ x_{2}^{T}& 1\\ \vdots& \vdots\\ x_{m}^{T}& 1\\ \end{matrix} \right]$
當 $X^TX$ 為滿秩矩陣的時候或者為正定矩陣的時候令上式為零我們可以很快速的求得參數。 $w^?=(XTX)?1XTy\hat{w}^*=\left( X^TX \right) ^{-1}X^Ty$
這里我們可以很快的求得多元線性回歸模型為： $f(x^i)=x^iT(XTX)?1XTyf\left( \hat{x}_i \right) =\hat{x}_{i}^{T}\left( X^TX \right) ^{-1}X^Ty$ 然而現實中 $X^TX$ 往往不是滿秩矩陣，例如很多任務中我們會遇到大量的變量，其數目超過了樣例數，導致 $X$ 的列數多于行數， $X^TX$ 顯然不是滿秩矩陣此時有多個 $w^\hat{w}$ ，它們都能夠是均方差最小化，這里將引入正則化(regularization)項。這里將先不論述正則化，待以后再寫一篇正則化的博客。

下面我們根據上述的要求來進行編程實現。

按照最小二乘法進行處理

from numpy import mat #先將訓練集x數據導出來做成矩陣的形式 x_train_dataset=train_dataset[:,0] #將一維數組轉換為矩陣 x_train_dataset=mat(x_train_dataset).reshape(300,1) a = np.ones(300) x_train_dataset=np.column_stack((x_train_dataset,a)) #print(x_train_dataset) print(x_train_dataset.shape)#將訓練集y數據導出來 y_train_dataset=train_dataset[:,1] y_train_dataset=mat(y_train_dataset).reshape(300,1) #print(y_train_dataset.shape) #print(y_train_dataset)#將訓練集x進行轉置 xt_train_dataset=np.transpose(x_train_dataset) print(xt_train_dataset.shape) #print(xt_train_dataset) (300, 2) (2, 300) # 利用公式求參數w w_predict = np.matmul(xt_train_dataset,x_train_dataset) #print(w_predict) w_predict = np.linalg.inv(w_predict) #print(w_predict) w_predict = np.matmul(w_predict,xt_train_dataset ) w_predict = np.matmul(w_predict,y_train_dataset ) print(w_predict) [[0.94732599][0.62680273]]

我們已經計算出參數，現在利用matplotlib將其可視化(散點圖和直線圖)

#重新讀取數據 x_train_dataset=train_dataset[:,0] y_train_dataset=train_dataset[:,1] #畫散點圖 fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.scatter(x_train_dataset,y_train_dataset,color='red',marker='o',s=1) ax.plot(x, 0.94732599*x+0.62680273,color='green',ls='-') #畫直線圖

上數我們計算出所需要的參數，這里利用測試集去測試所得參數的效果。我們計算出偏差和方差來對此種算法進行評估：

x_test_dataset=test_dataset[:,0] #print(y_test_dataset) x_test_dataset=mat(x_test_dataset).reshape(200,1) #print(y_test_dataset) y_test_dataset=test_dataset[:,1] y_test_dataset=mat(y_test_dataset).reshape(200,1) #print(y_test_dataset) import math bias2 = 0 #測試 #bias2=bias2+math.pow((x_test_dataset[1]*0.94732599+0.62680273)-y_test_dataset[1],2) #print(bias2) for i in range (0,200):bias2=bias2+math.pow((x_test_dataset[i]*0.94732599+0.62680273)-y_test_dataset[i],2) print(bias2) #之前看吳恩達的視頻其中講到要盡量避免用到loop語句， #下面我將對其進行修正 926.8174513840322

總結

以上是生活随笔為你收集整理的最小二乘法拟合直线（问题i）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。