前言

上一篇我們對數據進行了讀取并進行了可視化，今天我們來繼續實現算法。

完整代碼會在最后給出，如果你直接復制下面零散的代碼可能會運行不了。

這篇的代碼已經默認import了pandas,numpy等模塊。

數據標準化

首先我們先取要用的三列數據作為訓練數據集：

trainData = cReader[['mpg','displacement','acceleration']]

trainData.insert(0,'ones',1)

cols = trainData.shape[1]

X = trainData.iloc[:,0:cols-1]

Y = trainData.iloc[:,cols-1:cols]

X = np.mat(X.values)

Y = np.mat(Y.values)

for i in range(1,3):

X[:,i] = (X[:,i] - min(X[:,i])) / (max(X[:,i]) - min(X[:,I]))

Y[:,0] = (Y[:,0] - min(Y[:,0])) / (max(Y[:,0]) - min(Y[:,0]))

打印一下trainData前五行數據:

trainData數據

代碼說明

在訓練數據第一列加一列全為1的列矩陣trainData.insert(0,'ones',1)，目的是為了進行線性回歸的時候簡化計算，因為我們的目標方程是h(x)= θ0+ θ1𝓧1+ θ2𝓧2的一個常數項可以看成是θ0和1的乘積，其他項為θi和𝓧i的乘積

下圖βi即為我們的θi

解釋

cols得到trainData的列數，shape[0],shape[1]分別是trainData的行數和列數，下面把它的前三行'ones','mpg','displacement'分配給自變量矩陣𝓧最后一列分配給因變量矩陣Y

得到數據之后將它們標準化，關于做線性回歸是否需要標準化數據，這里有一個比較好的解釋

一般來說，我們再做線性回歸時并不需要中心化和標準化數據。大多數情況下數據中的特征會以不同的測量單位展現，無論有沒有中心化或者標準化都不會影響線性回歸的結果。

因為估計出來的參數值β會恰當地將每個解釋變量的單位x轉化為響應變量的單位y.

但是標準化數據能將我們的結果更具有可解釋性，比如β1=0.6

和β2=0.3, 我們可以理解為第一個特征的重要性是第二個特征的兩倍。

這里采用以下公式進行標準化，對應上面代碼的最后三行：

離差標準化公式

開始回歸

首先用最小二乘法計算回歸系數，并給出梯度下降的代價函數的代碼表示

最小二乘法公式:

最小二乘法

代碼:

theta_n = (X.T*X).I*X.T*Y

print(theta_n)

打印結果: [[ 0.58007057] [-0.0378746 ] [-0.35473364]]，從左到右分別是θ0，θ1，θ2

定義代價函數:

公式：

代價函數公式

代碼:

我是直接自己定義了一個新模塊linearRegrassion.py，在里面寫了線性回歸相關的函數，這也是上篇文章中開頭的import linearRegrassion as lg

記得在新文件linearRegrassion.py里要

import pandas as pd

import numpy as np

代價函數costFunc：

def costFunc(X,Y,theta):

inner = np.power((X*theta.T)-Y,2)

return np.sum(inner)/(2*len(X))

觀察公式，h(x)是預測值，y是實際值，通過他們倆的差來體現不同的擬合曲線的可靠程度，這個值越小說明對應的系數θ擬合出的曲線越接近實際，我們下面要做的就是用梯度下降的算法，取一系列θ值來逼近這個最優解，最后得到回歸曲線。

梯度下降算法

公式:

梯度下降公式

解釋

我們舉一個代價函數的例子，看下面這個函數圖像:

cost

我們在上面取一點(-30, 2500), 假設它對應一組我們待求的系數θ，此時代價函數值為2500，遠沒有達到最小的情況。

由于這個函數在(-∞, 20)區間內遞減，所以我們肯定需要增大θ值來逼近最優解，于是我們對θ求偏導，然后再用θ減去偏導值，這樣就可以使θ增大

(因為此時斜率k為負，所以減去之后肯定是增大的，相應的，如果這個點跑到對稱軸右邊，斜率變成正數，那么θ則會減少，仍然向最優解方向逼近)

說了這么多，我們取兩個點畫個圖來看看:

k

梯度下降迭代到最后，斜率越來越接近0，代價函數也越來越接近最優解，此時我們得到的便是擬合度最高的系數θ

上兩圖的代碼

import numpy as np

import matplotlib.pyplot as plt

x = np.linspace(-100,+100,10)

y = pow(x-20,2)

k = 2*(x-20)

y1 = -20*(x-10)+100

y2 = -100*(x+30)+2500

plt.plot(x,y)

plt.plot(x,y1,label= 'k=-20')

plt.plot(x,y2,label= 'k=-100')

leg = plt.legend(loc='upper right',ncol=1)

plt.show()

梯度下降代碼編寫,還是在linearRegrassion.py模塊中:

def gradientDescent(X,Y,theta,alpha,iters):

temp = np.mat(np.zeros(theta.shape))

cost = np.zeros(iters)

thetaNums = int(theta.shape[1])

print(thetaNums)

for i in range(iters):

error = (X*theta.T-Y)

for j in range(thetaNums):

derivativeInner = np.multiply(error,X[:,j])

temp[0,j] = theta[0,j] - (alpha*np.sum(derivativeInner)/len(X))

theta = temp

cost[i] = costFunc(X,Y,theta)

return theta,cost

解釋一下幾個參數:

alpha是公式中的𝒂，我們通常稱之為步長，它決定了θ逼近最優解的速度, 過大會導致函數無法收斂得不到最優解，過小則會使逼近速度過慢。

iters是迭代次數，足夠大的情況下得到的θ才有說服力

theta則是初始化迭代時給的一組θ值，最后得到擬合度最高的θ并在函數中返回

開始回歸

開始回歸之前，我們先設置一下學習的參數:

theta = np.mat([0,0,0])

iters = 100000

alpha = 0.001

回歸結果及可視化

finalTheta,cost = lg.gradientDescent(X,Y,theta,alpha,iters)

print(finalTheta)

print(cost)

x1 = np.linspace(X[:,1].min(),X[:,1].max(),100)

x2 = np.linspace(X[:,2].min(),X[:,2].max(),100)

x1,x2 = np.meshgrid(x1,x2)

f = finalTheta[0,0] + finalTheta[0,1]*x1 + finalTheta[0,2]*x2

fig = plt.figure()

Ax = Axes3D(fig)

Ax.plot_surface(x1, x2, f, rstride=1, cstride=1, cmap=cm.viridis,label='prediction')

Ax.scatter(X[:100,1],X[:100,2],Y[:100,0],c='y')

Ax.scatter(X[100:250,1],X[100:250,2],Y[100:250,0],c='r')

Ax.scatter(X[250:,1],X[250:,2],Y[250:,0],c='b')

Ax.set_zlabel('acceleration') # 坐標軸

Ax.set_ylabel('displacement')

Ax.set_xlabel('mpg')

plt.show()

finalTheta，也就是最終的θ:

[[ 15.47017446 2.99065096 -3.31870705]]

最小二乘法估計的結果

[[ 17.74518558]

[ -0.63629322]

[ -5.95952521]]

代價函數的值:

[ 124.67279846 124.37076615 124.06949161 ..., 2.77449658 2.77449481

2.77449304]

可以看到迭代十萬次之后得到的θ是和之前估計的比較接近的，如果你把迭代次數改為100萬次，雖然計算時間長一點但可以看到結果是基本一樣的

可視化:

回歸結果

比較密集的部分還是基本分布在平面上下的

梯度下降的過程可視化:

代碼：

fig, bx = plt.subplots(figsize=(8,6))

bx.plot(np.arange(iters), cost, 'r')

bx.set_xlabel('Iterations')

bx.set_ylabel('Cost')

bx.set_title('Error vs. Training Epoch')

plt.show()

梯度下降

如圖，隨著迭代次數的增加，代價函數越來越小，趨勢越來越平穩，同時逼近最優解。

完整代碼:

linearRegression.py

import pandas as pd

import numpy as np

def costFunc(X,Y,theta):

inner = np.power((X*theta.T)-Y,2)

return np.sum(inner)/(2*len(X))

def gradientDescent(X,Y,theta,alpha,iters):

temp = np.mat(np.zeros(theta.shape))

cost = np.zeros(iters)

thetaNums = int(theta.shape[1])

print(thetaNums)

for i in range(iters):

error = (X*theta.T-Y)

for j in range(thetaNums):

derivativeInner = np.multiply(error,X[:,j])

temp[0,j] = theta[0,j] - (alpha*np.sum(derivativeInner)/len(X))

theta = temp

cost[i] = costFunc(X,Y,theta)

return theta,cost

lgScript.py

from io import StringIO

from urllib import request

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from matplotlib import cm

import ssl

import pandas as pd

import numpy as np

import linearRegrassion as lg

ssl._create_default_https_context = ssl._create_unverified_context

names =["mpg","cylinders","displacement","horsepower",

"weight","acceleration","model year","origin","car name"]

url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data'

s = request.urlopen(url).read().decode('utf8')

dataFile = StringIO(s)

cReader = pd.read_csv(dataFile,delim_whitespace=True,names=names)

ax = plt.subplot(111, projection='3d') # 創建一個三維的繪圖工程

ax.scatter(cReader["mpg"][:100],cReader["displacement"][:100],cReader["acceleration"][:100],c='y')

ax.scatter(cReader["mpg"][100:250],cReader["displacement"][100:250],cReader["acceleration"][100:250],c='r')

ax.scatter(cReader["mpg"][250:],cReader["displacement"][250:],cReader["acceleration"][250:],c='b')

ax.set_zlabel('acceleration') # 坐標軸

ax.set_ylabel('displacement')

ax.set_xlabel('mpg')

plt.show()

plt.scatter(cReader["mpg"],cReader["displacement"])

plt.xlabel('mpg')

plt.ylabel('displacement')

plt.show()

trainData = cReader[['mpg','displacement','acceleration']]

trainData.insert(0,'ones',1)

print(trainData.head(5))

cols = trainData.shape[1]

X = trainData.iloc[:,0:cols-1]

Y = trainData.iloc[:,cols-1:cols]

X = np.mat(X.values)

Y = np.mat(Y.values)

for i in range(1,3):

X[:,i] = (X[:,i] - min(X[:,i])) / (max(X[:,i]) - min(X[:,i]))

print(X[:5:,:3])

#Y[:,0] = (Y[:,0] - min(Y[:,0])) / (max(Y[:,0]) - min(Y[:,0]))

print(Y[:5,0])

theta_n = (X.T*X).I*X.T*Y

print(theta_n)

theta = np.mat([0,0,0])

iters = 100000

alpha = 0.001

finalTheta,cost = lg.gradientDescent(X,Y,theta,alpha,iters)

print(finalTheta)

print(cost)

x1 = np.linspace(X[:,1].min(),X[:,1].max(),100)

x2 = np.linspace(X[:,2].min(),X[:,2].max(),100)

x1,x2 = np.meshgrid(x1,x2)

f = finalTheta[0,0] + finalTheta[0,1]*x1 + finalTheta[0,2]*x2

fig = plt.figure()

Ax = Axes3D(fig)

Ax.plot_surface(x1, x2, f, rstride=1, cstride=1, cmap=cm.viridis,label='prediction')

Ax.scatter(X[:100,1],X[:100,2],Y[:100,0],c='y')

Ax.scatter(X[100:250,1],X[100:250,2],Y[100:250,0],c='r')

Ax.scatter(X[250:,1],X[250:,2],Y[250:,0],c='b')

Ax.set_zlabel('acceleration') # 坐標軸

Ax.set_ylabel('displacement')

Ax.set_xlabel('mpg')

plt.show()

fig, bx = plt.subplots(figsize=(8,6))

bx.plot(np.arange(iters), cost, 'r')

bx.set_xlabel('Iterations')

bx.set_ylabel('Cost')

bx.set_title('Error vs. Training Epoch')

plt.show()

總結

這個系列算是學習吳恩達機器學習的筆記與作業，由于平時課程較多，所以不定期更新，下一篇大概是用Python實現邏輯回歸，希望不足之處大家可以討論一下。

總結

以上是生活随笔為你收集整理的python求解多元方程最优解_Python实现梯度下降算法求多元线性回归(二)的全部內容，希望文章能夠幫你解決所遇到的問題。

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