前言
- 上一篇我們對(duì)數(shù)據(jù)進(jìn)行了讀取并進(jìn)行了可視化狸相，今天我們來(lái)繼續(xù)實(shí)現(xiàn)算法。
- 完整代碼會(huì)在最后給出，如果你直接復(fù)制下面零散的代碼可能會(huì)運(yùn)行不了。
- 這篇的代碼已經(jīng)默認(rèn)import了pandas,numpy等模塊苔货。
數(shù)據(jù)標(biāo)準(zhǔn)化
- 首先我們先取要用的三列數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)集：
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前五行數(shù)據(jù):
trainData數(shù)據(jù)
代碼說(shuō)明
-
在訓(xùn)練數(shù)據(jù)第一列加一列全為1的列矩陣
trainData.insert(0,'ones',1)犀概，目的是為了進(jìn)行線性回歸的時(shí)候簡(jiǎn)化計(jì)算，因?yàn)槲覀兊哪繕?biāo)方程是h(x)= θ0+ θ1??1+ θ2??2的一個(gè)常數(shù)項(xiàng)可以看成是θ0和1的乘積夜惭，其他項(xiàng)為θi和??i的乘積
下圖βi即為我們的θi
解釋 cols得到trainData的列數(shù)姻灶，shape[0],shape[1]分別是trainData的行數(shù)和列數(shù)，下面把它的前三行
'ones','mpg','displacement'分配給自變量矩陣??最后一列分配給因變量矩陣Y得到數(shù)據(jù)之后將它們標(biāo)準(zhǔn)化诈茧，關(guān)于做線性回歸是否需要標(biāo)準(zhǔn)化數(shù)據(jù)产喉，這里有一個(gè)比較好的解釋
一般來(lái)說(shuō)，我們?cè)僮鼍€性回歸時(shí)并不需要中心化和標(biāo)準(zhǔn)化數(shù)據(jù)敢会。大多數(shù)情況下數(shù)據(jù)中的特征會(huì)以不同的測(cè)量單位展現(xiàn)曾沈，無(wú)論有沒(méi)有中心化或者標(biāo)準(zhǔn)化都不會(huì)影響線性回歸的結(jié)果。
因?yàn)楣烙?jì)出來(lái)的參數(shù)值β會(huì)恰當(dāng)?shù)貙⒚總€(gè)解釋變量的單位x轉(zhuǎn)化為響應(yīng)變量的單位y.
但是標(biāo)準(zhǔn)化數(shù)據(jù)能將我們的結(jié)果更具有可解釋性鸥昏，比如β1=0.6
和β2=0.3, 我們可以理解為第一個(gè)特征的重要性是第二個(gè)特征的兩倍塞俱。
這里采用以下公式進(jìn)行標(biāo)準(zhǔn)化，對(duì)應(yīng)上面代碼的最后三行：
離差標(biāo)準(zhǔn)化公式
開(kāi)始回歸
- 首先用最小二乘法計(jì)算回歸系數(shù)吏垮，并給出梯度下降的代價(jià)函數(shù)的代碼表示
-
最小二乘法公式:
最小二乘法 - 代碼:
theta_n = (X.T*X).I*X.T*Y
print(theta_n)
打印結(jié)果: [[ 0.58007057] [-0.0378746 ] [-0.35473364]]障涯，從左到右分別是θ0，θ1膳汪，θ2
- 定義代價(jià)函數(shù):
-
公式：
代價(jià)函數(shù)公式 - 代碼:
我是直接自己定義了一個(gè)新模塊linearRegrassion.py唯蝶，在里面寫(xiě)了線性回歸相關(guān)的函數(shù)，這也是上篇文章中開(kāi)頭的import linearRegrassion as lg
記得在新文件linearRegrassion.py里要
import pandas as pd
import numpy as np
代價(jià)函數(shù)costFunc：
def costFunc(X,Y,theta):
inner = np.power((X*theta.T)-Y,2)
return np.sum(inner)/(2*len(X))
- 觀察公式旅敷，h(x)是預(yù)測(cè)值生棍，y是實(shí)際值，通過(guò)他們倆的差來(lái)體現(xiàn)不同的擬合曲線的可靠程度媳谁，這個(gè)值越小說(shuō)明對(duì)應(yīng)的系數(shù)θ擬合出的曲線越接近實(shí)際涂滴，我們下面要做的就是用梯度下降的算法，取一系列θ值來(lái)逼近這個(gè)最優(yōu)解晴音，最后得到回歸曲線柔纵。
- 梯度下降算法
-
公式:
梯度下降公式 - 解釋
我們舉一個(gè)代價(jià)函數(shù)的例子，看下面這個(gè)函數(shù)圖像:
cost
我們?cè)谏厦嫒∫稽c(diǎn)(-30, 2500), 假設(shè)它對(duì)應(yīng)一組我們待求的系數(shù)θ锤躁，此時(shí)代價(jià)函數(shù)值為2500搁料，遠(yuǎn)沒(méi)有達(dá)到最小的情況。
由于這個(gè)函數(shù)在(-∞, 20)區(qū)間內(nèi)遞減系羞，所以我們肯定需要增大θ值來(lái)逼近最優(yōu)解郭计，于是我們對(duì)θ求偏導(dǎo)，然后再用θ減去偏導(dǎo)值椒振，這樣就可以使θ增大
(因?yàn)榇藭r(shí)斜率k為負(fù)昭伸，所以減去之后肯定是增大的，相應(yīng)的澎迎，如果這個(gè)點(diǎn)跑到對(duì)稱(chēng)軸右邊庐杨，斜率變成正數(shù)选调，那么θ則會(huì)減少，仍然向最優(yōu)解方向逼近)
說(shuō)了這么多灵份，我們?nèi)蓚€(gè)點(diǎn)畫(huà)個(gè)圖來(lái)看看:
k -
梯度下降迭代到最后仁堪，斜率越來(lái)越接近0，代價(jià)函數(shù)也越來(lái)越接近最優(yōu)解填渠，此時(shí)我們得到的便是擬合度最高的系數(shù)θ
上兩圖的代碼
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()
- 梯度下降代碼編寫(xiě),還是在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
- 解釋一下幾個(gè)參數(shù):
alpha是公式中的??弦聂，我們通常稱(chēng)之為步長(zhǎng)，它決定了θ逼近最優(yōu)解的速度, 過(guò)大會(huì)導(dǎo)致函數(shù)無(wú)法收斂得不到最優(yōu)解氛什，過(guò)小則會(huì)使逼近速度過(guò)慢横浑。
iters是迭代次數(shù)，足夠大的情況下得到的θ才有說(shuō)服力
theta則是初始化迭代時(shí)給的一組θ值屉更，最后得到擬合度最高的θ并在函數(shù)中返回
- 開(kāi)始回歸
- 開(kāi)始回歸之前，我們先設(shè)置一下學(xué)習(xí)的參數(shù):
theta = np.mat([0,0,0])
iters = 100000
alpha = 0.001
- 回歸結(jié)果及可視化
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') # 坐標(biāo)軸
Ax.set_ylabel('displacement')
Ax.set_xlabel('mpg')
plt.show()
finalTheta洒缀，也就是最終的θ:
[[ 15.47017446 2.99065096 -3.31870705]]
最小二乘法估計(jì)的結(jié)果
[[ 17.74518558]
[ -0.63629322]
[ -5.95952521]]
代價(jià)函數(shù)的值:
[ 124.67279846 124.37076615 124.06949161 ..., 2.77449658 2.77449481
2.77449304]
- 可以看到迭代十萬(wàn)次之后得到的θ是和之前估計(jì)的比較接近的瑰谜，如果你把迭代次數(shù)改為100萬(wàn)次，雖然計(jì)算時(shí)間長(zhǎng)一點(diǎn)但可以看到結(jié)果是基本一樣的
-
可視化:
回歸結(jié)果
比較密集的部分還是基本分布在平面上下的
- 梯度下降的過(guò)程可視化:
代碼：
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()
梯度下降
如圖树绩，隨著迭代次數(shù)的增加萨脑，代價(jià)函數(shù)越來(lái)越小，趨勢(shì)越來(lái)越平穩(wěn)饺饭，同時(shí)逼近最優(yōu)解渤早。
- 完整代碼:
- 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') # 創(chuàng)建一個(gè)三維的繪圖工程
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') # 坐標(biāo)軸
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') # 坐標(biāo)軸
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()
總結(jié)
這個(gè)系列算是學(xué)習(xí)吳恩達(dá)機(jī)器學(xué)習(xí)的筆記與作業(yè)，由于平時(shí)課程較多瘫俊，所以不定期更新鹊杖，下一篇大概是用Python實(shí)現(xiàn)邏輯回歸，希望不足之處大家可以討論一下扛芽。
