决策算法python_GitHub - nxety/MachineLearning_Python: 机器学习算法python实现
機(jī)器學(xué)習(xí)算法Python實(shí)現(xiàn)
目錄
1、代價(jià)函數(shù)
其中:
下面就是要求出theta,使代價(jià)最小,即代表我們擬合出來的方程距離真實(shí)值最近
共有m條數(shù)據(jù),其中代表我們要擬合出來的方程到真實(shí)值距離的平方,平方的原因是因?yàn)榭赡苡胸?fù)值,正負(fù)可能會(huì)抵消
前面有系數(shù)2的原因是下面求梯度是對(duì)每個(gè)變量求偏導(dǎo),2可以消去
實(shí)現(xiàn)代碼:
# 計(jì)算代價(jià)函數(shù)
def computerCost(X,y,theta):
m = len(y)
J = 0
J = (np.transpose(X*theta-y))*(X*theta-y)/(2*m) #計(jì)算代價(jià)J
return J
注意這里的X是真實(shí)數(shù)據(jù)前加了一列1,因?yàn)橛衪heta(0)
2、梯度下降算法
代價(jià)函數(shù)對(duì)求偏導(dǎo)得到:
所以對(duì)theta的更新可以寫為:
其中為學(xué)習(xí)速率,控制梯度下降的速度,一般取0.01,0.03,0.1,0.3.....
為什么梯度下降可以逐步減小代價(jià)函數(shù)
假設(shè)函數(shù)f(x)
泰勒展開:f(x+△x)=f(x)+f'(x)*△x+o(△x)
令:△x=-α*f'(x) ,即負(fù)梯度方向乘以一個(gè)很小的步長α
將△x代入泰勒展開式中:f(x+△x)=f(x)-α*[f'(x)]2+o(△x)
可以看出,α是取得很小的正數(shù),[f'(x)]2也是正數(shù),所以可以得出:f(x+△x)<=f(x)
所以沿著負(fù)梯度放下,函數(shù)在減小,多維情況一樣。
實(shí)現(xiàn)代碼
# 梯度下降算法
def gradientDescent(X,y,theta,alpha,num_iters):
m = len(y)
n = len(theta)
temp = np.matrix(np.zeros((n,num_iters))) # 暫存每次迭代計(jì)算的theta,轉(zhuǎn)化為矩陣形式
J_history = np.zeros((num_iters,1)) #記錄每次迭代計(jì)算的代價(jià)值
for i in range(num_iters): # 遍歷迭代次數(shù)
h = np.dot(X,theta) # 計(jì)算內(nèi)積,matrix可以直接乘
temp[:,i] = theta - ((alpha/m)*(np.dot(np.transpose(X),h-y))) #梯度的計(jì)算
theta = temp[:,i]
J_history[i] = computerCost(X,y,theta) #調(diào)用計(jì)算代價(jià)函數(shù)
print '.',
return theta,J_history
3、均值歸一化
目的是使數(shù)據(jù)都縮放到一個(gè)范圍內(nèi),便于使用梯度下降算法
其中 為所有此feture數(shù)據(jù)的平均值
可以是最大值-最小值,也可以是這個(gè)feature對(duì)應(yīng)的數(shù)據(jù)的標(biāo)準(zhǔn)差
實(shí)現(xiàn)代碼:
# 歸一化feature
def featureNormaliza(X):
X_norm = np.array(X) #將X轉(zhuǎn)化為numpy數(shù)組對(duì)象,才可以進(jìn)行矩陣的運(yùn)算
#定義所需變量
mu = np.zeros((1,X.shape[1]))
sigma = np.zeros((1,X.shape[1]))
mu = np.mean(X_norm,0) # 求每一列的平均值(0指定為列,1代表行)
sigma = np.std(X_norm,0) # 求每一列的標(biāo)準(zhǔn)差
for i in range(X.shape[1]): # 遍歷列
X_norm[:,i] = (X_norm[:,i]-mu[i])/sigma[i] # 歸一化
return X_norm,mu,sigma
注意預(yù)測的時(shí)候也需要均值歸一化數(shù)據(jù)
4、最終運(yùn)行結(jié)果
代價(jià)隨迭代次數(shù)的變化
導(dǎo)入包
from sklearn import linear_model
from sklearn.preprocessing import StandardScaler #引入縮放的包
歸一化
# 歸一化操作
scaler = StandardScaler()
scaler.fit(X)
x_train = scaler.transform(X)
x_test = scaler.transform(np.array([1650,3]))
線性模型擬合
# 線性模型擬合
model = linear_model.LinearRegression()
model.fit(x_train, y)
預(yù)測
#預(yù)測結(jié)果
result = model.predict(x_test)
1、代價(jià)函數(shù)
可以綜合起來為:
其中:
為什么不用線性回歸的代價(jià)函數(shù)表示,因?yàn)榫€性回歸的代價(jià)函數(shù)可能是非凸的,對(duì)于分類問題,使用梯度下降很難得到最小值,上面的代價(jià)函數(shù)是凸函數(shù)
的圖像如下,即y=1時(shí):
可以看出,當(dāng)趨于1,y=1,與預(yù)測值一致,此時(shí)付出的代價(jià)cost趨于0,若趨于0,y=1,此時(shí)的代價(jià)cost值非常大,我們最終的目的是最小化代價(jià)值
同理的圖像如下(y=0):
2、梯度
同樣對(duì)代價(jià)函數(shù)求偏導(dǎo):
可以看出與線性回歸的偏導(dǎo)數(shù)一致
推到過程
3、正則化
目的是為了防止過擬合
在代價(jià)函數(shù)中加上一項(xiàng)
注意j是重1開始的,因?yàn)閠heta(0)為一個(gè)常數(shù)項(xiàng),X中最前面一列會(huì)加上1列1,所以乘積還是theta(0),feature沒有關(guān)系,沒有必要正則化
正則化后的代價(jià):
# 代價(jià)函數(shù)
def costFunction(initial_theta,X,y,inital_lambda):
m = len(y)
J = 0
h = sigmoid(np.dot(X,initial_theta)) # 計(jì)算h(z)
theta1 = initial_theta.copy() # 因?yàn)檎齽t化j=1從1開始,不包含0,所以復(fù)制一份,前theta(0)值為0
theta1[0] = 0
temp = np.dot(np.transpose(theta1),theta1)
J = (-np.dot(np.transpose(y),np.log(h))-np.dot(np.transpose(1-y),np.log(1-h))+temp*inital_lambda/2)/m # 正則化的代價(jià)方程
return J
正則化后的代價(jià)的梯度
# 計(jì)算梯度
def gradient(initial_theta,X,y,inital_lambda):
m = len(y)
grad = np.zeros((initial_theta.shape[0]))
h = sigmoid(np.dot(X,initial_theta))# 計(jì)算h(z)
theta1 = initial_theta.copy()
theta1[0] = 0
grad = np.dot(np.transpose(X),h-y)/m+inital_lambda/m*theta1 #正則化的梯度
return grad
4、S型函數(shù)(即)
實(shí)現(xiàn)代碼:
# S型函數(shù)
def sigmoid(z):
h = np.zeros((len(z),1)) # 初始化,與z的長度一置
h = 1.0/(1.0+np.exp(-z))
return h
5、映射為多項(xiàng)式
因?yàn)閿?shù)據(jù)的feture可能很少,導(dǎo)致偏差大,所以創(chuàng)造出一些feture結(jié)合
eg:映射為2次方的形式:
實(shí)現(xiàn)代碼:
# 映射為多項(xiàng)式
def mapFeature(X1,X2):
degree = 3; # 映射的最高次方
out = np.ones((X1.shape[0],1)) # 映射后的結(jié)果數(shù)組(取代X)
'''
這里以degree=2為例,映射為1,x1,x2,x1^2,x1,x2,x2^2
'''
for i in np.arange(1,degree+1):
for j in range(i+1):
temp = X1**(i-j)*(X2**j) #矩陣直接乘相當(dāng)于matlab中的點(diǎn)乘.*
out = np.hstack((out, temp.reshape(-1,1)))
return out
6、使用scipy的優(yōu)化方法
梯度下降使用scipy中optimize中的fmin_bfgs函數(shù)
調(diào)用scipy中的優(yōu)化算法fmin_bfgs(擬牛頓法Broyden-Fletcher-Goldfarb-Shanno
costFunction是自己實(shí)現(xiàn)的一個(gè)求代價(jià)的函數(shù),
initial_theta表示初始化的值,
fprime指定costFunction的梯度
args是其余測參數(shù),以元組的形式傳入,最后會(huì)將最小化costFunction的theta返回
result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, args=(X,y,initial_lambda))
7、運(yùn)行結(jié)果
data1決策邊界和準(zhǔn)確度
data2決策邊界和準(zhǔn)確度
導(dǎo)入包
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.cross_validation import train_test_split
import numpy as np
劃分訓(xùn)練集和測試集
# 劃分為訓(xùn)練集和測試集
x_train,x_test,y_train,y_test = train_test_split(X,y,test_size=0.2)
歸一化
# 歸一化
scaler = StandardScaler()
x_train = scaler.fit_transform(x_train)
x_test = scaler.fit_transform(x_test)
邏輯回歸
#邏輯回歸
model = LogisticRegression()
model.fit(x_train,y_train)
預(yù)測
# 預(yù)測
predict = model.predict(x_test)
right = sum(predict == y_test)
predict = np.hstack((predict.reshape(-1,1),y_test.reshape(-1,1))) # 將預(yù)測值和真實(shí)值放在一塊,好觀察
print predict
print ('測試集準(zhǔn)確率:%f%%'%(right*100.0/predict.shape[0])) #計(jì)算在測試集上的準(zhǔn)確度
1、隨機(jī)顯示100個(gè)數(shù)字
我沒有使用scikit-learn中的數(shù)據(jù)集,像素是20*20px,彩色圖如下
灰度圖:
實(shí)現(xiàn)代碼:
# 顯示100個(gè)數(shù)字
def display_data(imgData):
sum = 0
'''
顯示100個(gè)數(shù)(若是一個(gè)一個(gè)繪制將會(huì)非常慢,可以將要畫的數(shù)字整理好,放到一個(gè)矩陣中,顯示這個(gè)矩陣即可)
- 初始化一個(gè)二維數(shù)組
- 將每行的數(shù)據(jù)調(diào)整成圖像的矩陣,放進(jìn)二維數(shù)組
- 顯示即可
'''
pad = 1
display_array = -np.ones((pad+10*(20+pad),pad+10*(20+pad)))
for i in range(10):
for j in range(10):
display_array[pad+i*(20+pad):pad+i*(20+pad)+20,pad+j*(20+pad):pad+j*(20+pad)+20] = (imgData[sum,:].reshape(20,20,order="F")) # order=F指定以列優(yōu)先,在matlab中是這樣的,python中需要指定,默認(rèn)以行
sum += 1
plt.imshow(display_array,cmap='gray') #顯示灰度圖像
plt.axis('off')
plt.show()
2、OneVsAll
如何利用邏輯回歸解決多分類的問題,OneVsAll就是把當(dāng)前某一類看成一類,其他所有類別看作一類,這樣有成了二分類的問題了
如下圖,把途中的數(shù)據(jù)分成三類,先把紅色的看成一類,把其他的看作另外一類,進(jìn)行邏輯回歸,然后把藍(lán)色的看成一類,其他的再看成一類,以此類推...
可以看出大于2類的情況下,有多少類就要進(jìn)行多少次的邏輯回歸分類
3、手寫數(shù)字識(shí)別
共有0-9,10個(gè)數(shù)字,需要10次分類
由于數(shù)據(jù)集y給出的是0,1,2...9的數(shù)字,而進(jìn)行邏輯回歸需要0/1的label標(biāo)記,所以需要對(duì)y處理
說一下數(shù)據(jù)集,前500個(gè)是0,500-1000是1,...,所以如下圖,處理后的y,前500行的第一列是1,其余都是0,500-1000行第二列是1,其余都是0....
然后調(diào)用梯度下降算法求解theta
實(shí)現(xiàn)代碼:
# 求每個(gè)分類的theta,最后返回所有的all_theta
def oneVsAll(X,y,num_labels,Lambda):
# 初始化變量
m,n = X.shape
all_theta = np.zeros((n+1,num_labels)) # 每一列對(duì)應(yīng)相應(yīng)分類的theta,共10列
X = np.hstack((np.ones((m,1)),X)) # X前補(bǔ)上一列1的偏置bias
class_y = np.zeros((m,num_labels)) # 數(shù)據(jù)的y對(duì)應(yīng)0-9,需要映射為0/1的關(guān)系
initial_theta = np.zeros((n+1,1)) # 初始化一個(gè)分類的theta
# 映射y
for i in range(num_labels):
class_y[:,i] = np.int32(y==i).reshape(1,-1) # 注意reshape(1,-1)才可以賦值
#np.savetxt("class_y.csv", class_y[0:600,:], delimiter=',')
'''遍歷每個(gè)分類,計(jì)算對(duì)應(yīng)的theta值'''
for i in range(num_labels):
result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, args=(X,class_y[:,i],Lambda)) # 調(diào)用梯度下降的優(yōu)化方法
all_theta[:,i] = result.reshape(1,-1) # 放入all_theta中
all_theta = np.transpose(all_theta)
return all_theta
4、預(yù)測
之前說過,預(yù)測的結(jié)果是一個(gè)概率值,利用學(xué)習(xí)出來的theta代入預(yù)測的S型函數(shù)中,每行的最大值就是是某個(gè)數(shù)字的最大概率,所在的列號(hào)就是預(yù)測的數(shù)字的真實(shí)值,因?yàn)樵诜诸悤r(shí),所有為0的將y映射在第一列,為1的映射在第二列,依次類推
實(shí)現(xiàn)代碼:
# 預(yù)測
def predict_oneVsAll(all_theta,X):
m = X.shape[0]
num_labels = all_theta.shape[0]
p = np.zeros((m,1))
X = np.hstack((np.ones((m,1)),X)) #在X最前面加一列1
h = sigmoid(np.dot(X,np.transpose(all_theta))) #預(yù)測
'''
返回h中每一行最大值所在的列號(hào)
- np.max(h, axis=1)返回h中每一行的最大值(是某個(gè)數(shù)字的最大概率)
- 最后where找到的最大概率所在的列號(hào)(列號(hào)即是對(duì)應(yīng)的數(shù)字)
'''
p = np.array(np.where(h[0,:] == np.max(h, axis=1)[0]))
for i in np.arange(1, m):
t = np.array(np.where(h[i,:] == np.max(h, axis=1)[i]))
p = np.vstack((p,t))
return p
5、運(yùn)行結(jié)果
10次分類,在訓(xùn)練集上的準(zhǔn)確度:
1、導(dǎo)入包
from scipy import io as spio
import numpy as np
from sklearn import svm
from sklearn.linear_model import LogisticRegression
2、加載數(shù)據(jù)
data = loadmat_data("data_digits.mat")
X = data['X'] # 獲取X數(shù)據(jù),每一行對(duì)應(yīng)一個(gè)數(shù)字20x20px
y = data['y'] # 這里讀取mat文件y的shape=(5000, 1)
y = np.ravel(y) # 調(diào)用sklearn需要轉(zhuǎn)化成一維的(5000,)
3、擬合模型
model = LogisticRegression()
model.fit(X, y) # 擬合
4、預(yù)測
predict = model.predict(X) #預(yù)測
print u"預(yù)測準(zhǔn)確度為:%f%%"%np.mean(np.float64(predict == y)*100)
5、輸出結(jié)果(在訓(xùn)練集上的準(zhǔn)確度)
三、BP神經(jīng)網(wǎng)絡(luò)
1、神經(jīng)網(wǎng)絡(luò)model
先介紹個(gè)三層的神經(jīng)網(wǎng)絡(luò),如下圖所示
輸入層(input layer)有三個(gè)units(為補(bǔ)上的bias,通常設(shè)為1)
表示第j層的第i個(gè)激勵(lì),也稱為為單元unit
為第j層到第j+1層映射的權(quán)重矩陣,就是每條邊的權(quán)重
所以可以得到:
隱含層:
輸出層
其中,S型函數(shù),也成為激勵(lì)函數(shù)
可以看出 為3x4的矩陣,為1x4的矩陣
==》j+1的單元數(shù)x(j層的單元數(shù)+1)
2、代價(jià)函數(shù)
假設(shè)最后輸出的,即代表輸出層有K個(gè)單元
其中,代表第i個(gè)單元輸出
與邏輯回歸的代價(jià)函數(shù)差不多,就是累加上每個(gè)輸出(共有K個(gè)輸出)
3、正則化
L-->所有層的個(gè)數(shù)
-->第l層unit的個(gè)數(shù)
正則化后的代價(jià)函數(shù)為
共有L-1層,
然后是累加對(duì)應(yīng)每一層的theta矩陣,注意不包含加上偏置項(xiàng)對(duì)應(yīng)的theta(0)
正則化后的代價(jià)函數(shù)實(shí)現(xiàn)代碼:
# 代價(jià)函數(shù)
def nnCostFunction(nn_params,input_layer_size,hidden_layer_size,num_labels,X,y,Lambda):
length = nn_params.shape[0] # theta的中長度
# 還原theta1和theta2
Theta1 = nn_params[0:hidden_layer_size*(input_layer_size+1)].reshape(hidden_layer_size,input_layer_size+1)
Theta2 = nn_params[hidden_layer_size*(input_layer_size+1):length].reshape(num_labels,hidden_layer_size+1)
# np.savetxt("Theta1.csv",Theta1,delimiter=',')
m = X.shape[0]
class_y = np.zeros((m,num_labels)) # 數(shù)據(jù)的y對(duì)應(yīng)0-9,需要映射為0/1的關(guān)系
# 映射y
for i in range(num_labels):
class_y[:,i] = np.int32(y==i).reshape(1,-1) # 注意reshape(1,-1)才可以賦值
'''去掉theta1和theta2的第一列,因?yàn)檎齽t化時(shí)從1開始'''
Theta1_colCount = Theta1.shape[1]
Theta1_x = Theta1[:,1:Theta1_colCount]
Theta2_colCount = Theta2.shape[1]
Theta2_x = Theta2[:,1:Theta2_colCount]
# 正則化向theta^2
term = np.dot(np.transpose(np.vstack((Theta1_x.reshape(-1,1),Theta2_x.reshape(-1,1)))),np.vstack((Theta1_x.reshape(-1,1),Theta2_x.reshape(-1,1))))
'''正向傳播,每次需要補(bǔ)上一列1的偏置bias'''
a1 = np.hstack((np.ones((m,1)),X))
z2 = np.dot(a1,np.transpose(Theta1))
a2 = sigmoid(z2)
a2 = np.hstack((np.ones((m,1)),a2))
z3 = np.dot(a2,np.transpose(Theta2))
h = sigmoid(z3)
'''代價(jià)'''
J = -(np.dot(np.transpose(class_y.reshape(-1,1)),np.log(h.reshape(-1,1)))+np.dot(np.transpose(1-class_y.reshape(-1,1)),np.log(1-h.reshape(-1,1)))-Lambda*term/2)/m
return np.ravel(J)
4、反向傳播BP
上面正向傳播可以計(jì)算得到J(θ),使用梯度下降法還需要求它的梯度
BP反向傳播的目的就是求代價(jià)函數(shù)的梯度
假設(shè)4層的神經(jīng)網(wǎng)絡(luò),記為-->l層第j個(gè)單元的誤差
《===》(向量化)
沒有,因?yàn)閷?duì)于輸入沒有誤差
因?yàn)镾型函數(shù)的導(dǎo)數(shù)為:,所以上面的和可以在前向傳播中計(jì)算出來
反向傳播計(jì)算梯度的過程為:
(是大寫的)
for i=1-m:
-
-正向傳播計(jì)算(l=2,3,4...L)
-反向計(jì)算、...;
-
-
最后,即得到代價(jià)函數(shù)的梯度
實(shí)現(xiàn)代碼:
# 梯度
def nnGradient(nn_params,input_layer_size,hidden_layer_size,num_labels,X,y,Lambda):
length = nn_params.shape[0]
Theta1 = nn_params[0:hidden_layer_size*(input_layer_size+1)].reshape(hidden_layer_size,input_layer_size+1).copy() # 這里使用copy函數(shù),否則下面修改Theta的值,nn_params也會(huì)一起修改
Theta2 = nn_params[hidden_layer_size*(input_layer_size+1):length].reshape(num_labels,hidden_layer_size+1).copy()
m = X.shape[0]
class_y = np.zeros((m,num_labels)) # 數(shù)據(jù)的y對(duì)應(yīng)0-9,需要映射為0/1的關(guān)系
# 映射y
for i in range(num_labels):
class_y[:,i] = np.int32(y==i).reshape(1,-1) # 注意reshape(1,-1)才可以賦值
'''去掉theta1和theta2的第一列,因?yàn)檎齽t化時(shí)從1開始'''
Theta1_colCount = Theta1.shape[1]
Theta1_x = Theta1[:,1:Theta1_colCount]
Theta2_colCount = Theta2.shape[1]
Theta2_x = Theta2[:,1:Theta2_colCount]
Theta1_grad = np.zeros((Theta1.shape)) #第一層到第二層的權(quán)重
Theta2_grad = np.zeros((Theta2.shape)) #第二層到第三層的權(quán)重
'''正向傳播,每次需要補(bǔ)上一列1的偏置bias'''
a1 = np.hstack((np.ones((m,1)),X))
z2 = np.dot(a1,np.transpose(Theta1))
a2 = sigmoid(z2)
a2 = np.hstack((np.ones((m,1)),a2))
z3 = np.dot(a2,np.transpose(Theta2))
h = sigmoid(z3)
'''反向傳播,delta為誤差,'''
delta3 = np.zeros((m,num_labels))
delta2 = np.zeros((m,hidden_layer_size))
for i in range(m):
#delta3[i,:] = (h[i,:]-class_y[i,:])*sigmoidGradient(z3[i,:]) # 均方誤差的誤差率
delta3[i,:] = h[i,:]-class_y[i,:] # 交叉熵誤差率
Theta2_grad = Theta2_grad+np.dot(np.transpose(delta3[i,:].reshape(1,-1)),a2[i,:].reshape(1,-1))
delta2[i,:] = np.dot(delta3[i,:].reshape(1,-1),Theta2_x)*sigmoidGradient(z2[i,:])
Theta1_grad = Theta1_grad+np.dot(np.transpose(delta2[i,:].reshape(1,-1)),a1[i,:].reshape(1,-1))
Theta1[:,0] = 0
Theta2[:,0] = 0
'''梯度'''
grad = (np.vstack((Theta1_grad.reshape(-1,1),Theta2_grad.reshape(-1,1)))+Lambda*np.vstack((Theta1.reshape(-1,1),Theta2.reshape(-1,1))))/m
return np.ravel(grad)
5、BP可以求梯度的原因
實(shí)際是利用了鏈?zhǔn)角髮?dǎo)法則
因?yàn)橄乱粚拥膯卧蒙弦粚拥膯卧鳛檩斎脒M(jìn)行計(jì)算
大體的推導(dǎo)過程如下,最終我們是想預(yù)測函數(shù)與已知的y非常接近,求均方差的梯度沿著此梯度方向可使代價(jià)函數(shù)最小化。可對(duì)照上面求梯度的過程。
求誤差更詳細(xì)的推導(dǎo)過程:
6、梯度檢查
檢查利用BP求的梯度是否正確
利用導(dǎo)數(shù)的定義驗(yàn)證:
求出來的數(shù)值梯度應(yīng)該與BP求出的梯度非常接近
驗(yàn)證BP正確后就不需要再執(zhí)行驗(yàn)證梯度的算法了
實(shí)現(xiàn)代碼:
# 檢驗(yàn)梯度是否計(jì)算正確
# 檢驗(yàn)梯度是否計(jì)算正確
def checkGradient(Lambda = 0):
'''構(gòu)造一個(gè)小型的神經(jīng)網(wǎng)絡(luò)驗(yàn)證,因?yàn)閿?shù)值法計(jì)算梯度很浪費(fèi)時(shí)間,而且驗(yàn)證正確后之后就不再需要驗(yàn)證了'''
input_layer_size = 3
hidden_layer_size = 5
num_labels = 3
m = 5
initial_Theta1 = debugInitializeWeights(input_layer_size,hidden_layer_size);
initial_Theta2 = debugInitializeWeights(hidden_layer_size,num_labels)
X = debugInitializeWeights(input_layer_size-1,m)
y = 1+np.transpose(np.mod(np.arange(1,m+1), num_labels))# 初始化y
y = y.reshape(-1,1)
nn_params = np.vstack((initial_Theta1.reshape(-1,1),initial_Theta2.reshape(-1,1))) #展開theta
'''BP求出梯度'''
grad = nnGradient(nn_params, input_layer_size, hidden_layer_size,
num_labels, X, y, Lambda)
'''使用數(shù)值法計(jì)算梯度'''
num_grad = np.zeros((nn_params.shape[0]))
step = np.zeros((nn_params.shape[0]))
e = 1e-4
for i in range(nn_params.shape[0]):
step[i] = e
loss1 = nnCostFunction(nn_params-step.reshape(-1,1), input_layer_size, hidden_layer_size,
num_labels, X, y,
Lambda)
loss2 = nnCostFunction(nn_params+step.reshape(-1,1), input_layer_size, hidden_layer_size,
num_labels, X, y,
Lambda)
num_grad[i] = (loss2-loss1)/(2*e)
step[i]=0
# 顯示兩列比較
res = np.hstack((num_grad.reshape(-1,1),grad.reshape(-1,1)))
print res
7、權(quán)重的隨機(jī)初始化
神經(jīng)網(wǎng)絡(luò)不能像邏輯回歸那樣初始化theta為0,因?yàn)槿羰敲織l邊的權(quán)重都為0,每個(gè)神經(jīng)元都是相同的輸出,在反向傳播中也會(huì)得到同樣的梯度,最終只會(huì)預(yù)測一種結(jié)果。
所以應(yīng)該初始化為接近0的數(shù)
實(shí)現(xiàn)代碼
# 隨機(jī)初始化權(quán)重theta
def randInitializeWeights(L_in,L_out):
W = np.zeros((L_out,1+L_in)) # 對(duì)應(yīng)theta的權(quán)重
epsilon_init = (6.0/(L_out+L_in))**0.5
W = np.random.rand(L_out,1+L_in)*2*epsilon_init-epsilon_init # np.random.rand(L_out,1+L_in)產(chǎn)生L_out*(1+L_in)大小的隨機(jī)矩陣
return W
8、預(yù)測
正向傳播預(yù)測結(jié)果
實(shí)現(xiàn)代碼
# 預(yù)測
def predict(Theta1,Theta2,X):
m = X.shape[0]
num_labels = Theta2.shape[0]
#p = np.zeros((m,1))
'''正向傳播,預(yù)測結(jié)果'''
X = np.hstack((np.ones((m,1)),X))
h1 = sigmoid(np.dot(X,np.transpose(Theta1)))
h1 = np.hstack((np.ones((m,1)),h1))
h2 = sigmoid(np.dot(h1,np.transpose(Theta2)))
'''
返回h中每一行最大值所在的列號(hào)
- np.max(h, axis=1)返回h中每一行的最大值(是某個(gè)數(shù)字的最大概率)
- 最后where找到的最大概率所在的列號(hào)(列號(hào)即是對(duì)應(yīng)的數(shù)字)
'''
#np.savetxt("h2.csv",h2,delimiter=',')
p = np.array(np.where(h2[0,:] == np.max(h2, axis=1)[0]))
for i in np.arange(1, m):
t = np.array(np.where(h2[i,:] == np.max(h2, axis=1)[i]))
p = np.vstack((p,t))
return p
9、輸出結(jié)果
梯度檢查:
隨機(jī)顯示100個(gè)手寫數(shù)字
顯示theta1權(quán)重
訓(xùn)練集預(yù)測準(zhǔn)確度
歸一化后訓(xùn)練集預(yù)測準(zhǔn)確度
四、SVM支持向量機(jī)
1、代價(jià)函數(shù)
在邏輯回歸中,我們的代價(jià)為:
,
其中:,
如圖所示,如果y=1,cost代價(jià)函數(shù)如圖所示
我們想讓,即z>>0,這樣的話cost代價(jià)函數(shù)才會(huì)趨于最小(這是我們想要的),所以用途中紅色的函數(shù)代替邏輯回歸中的cost
當(dāng)y=0時(shí)同樣,用代替
最終得到的代價(jià)函數(shù)為:
最后我們想要
之前我們邏輯回歸中的代價(jià)函數(shù)為:
可以認(rèn)為這里的,只是表達(dá)形式問題,這里C的值越大,SVM的決策邊界的margin也越大,下面會(huì)說明
2、Large Margin
如下圖所示,SVM分類會(huì)使用最大的margin將其分開
先說一下向量內(nèi)積
,
表示u的歐幾里得范數(shù)(歐式范數(shù)),
向量V在向量u上的投影的長度記為p,則:向量內(nèi)積:
根據(jù)向量夾角公式推導(dǎo)一下即可,
前面說過,當(dāng)C越大時(shí),margin也就越大,我們的目的是最小化代價(jià)函數(shù)J(θ),當(dāng)margin最大時(shí),C的乘積項(xiàng)要很小,所以近似為:
,
我們最后的目的就是求使代價(jià)最小的θ
由
可以得到:
,p即為x在θ上的投影
如下圖所示,假設(shè)決策邊界如圖,找其中的一個(gè)點(diǎn),到θ上的投影為p,則或者,若是p很小,則需要很大,這與我們要求的θ使最小相違背,所以最后求的是large margin
3、SVM Kernel(核函數(shù))
對(duì)于線性可分的問題,使用線性核函數(shù)即可
對(duì)于線性不可分的問題,在邏輯回歸中,我們是將feature映射為使用多項(xiàng)式的形式,SVM中也有多項(xiàng)式核函數(shù),但是更常用的是高斯核函數(shù),也稱為RBF核
高斯核函數(shù)為:
假設(shè)如圖幾個(gè)點(diǎn),
令:
.
.
.
可以看出,若是x與距離較近,==》,(即相似度較大)
若是x與距離較遠(yuǎn),==》,(即相似度較低)
高斯核函數(shù)的σ越小,f下降的越快
如何選擇初始的
訓(xùn)練集:
選擇:
對(duì)于給出的x,計(jì)算f,令:所以:
最小化J求出θ,
如果,==》預(yù)測y=1
4、使用scikit-learn中的SVM模型代碼
線性可分的,指定核函數(shù)為linear:
'''data1——線性分類'''
data1 = spio.loadmat('data1.mat')
X = data1['X']
y = data1['y']
y = np.ravel(y)
plot_data(X,y)
model = svm.SVC(C=1.0,kernel='linear').fit(X,y) # 指定核函數(shù)為線性核函數(shù)
非線性可分的,默認(rèn)核函數(shù)為rbf
'''data2——非線性分類'''
data2 = spio.loadmat('data2.mat')
X = data2['X']
y = data2['y']
y = np.ravel(y)
plt = plot_data(X,y)
plt.show()
model = svm.SVC(gamma=100).fit(X,y) # gamma為核函數(shù)的系數(shù),值越大擬合的越好
5、運(yùn)行結(jié)果
線性可分的決策邊界:
線性不可分的決策邊界:
五、K-Means聚類算法
1、聚類過程
聚類屬于無監(jiān)督學(xué)習(xí),不知道y的標(biāo)記分為K類
K-Means算法分為兩個(gè)步驟
第一步:簇分配,隨機(jī)選K個(gè)點(diǎn)作為中心,計(jì)算到這K個(gè)點(diǎn)的距離,分為K個(gè)簇
第二步:移動(dòng)聚類中心:重新計(jì)算每個(gè)簇的中心,移動(dòng)中心,重復(fù)以上步驟。
如下圖所示:
隨機(jī)分配的聚類中心
重新計(jì)算聚類中心,移動(dòng)一次
最后10步之后的聚類中心
計(jì)算每條數(shù)據(jù)到哪個(gè)中心最近實(shí)現(xiàn)代碼:
# 找到每條數(shù)據(jù)距離哪個(gè)類中心最近
def findClosestCentroids(X,initial_centroids):
m = X.shape[0] # 數(shù)據(jù)條數(shù)
K = initial_centroids.shape[0] # 類的總數(shù)
dis = np.zeros((m,K)) # 存儲(chǔ)計(jì)算每個(gè)點(diǎn)分別到K個(gè)類的距離
idx = np.zeros((m,1)) # 要返回的每條數(shù)據(jù)屬于哪個(gè)類
'''計(jì)算每個(gè)點(diǎn)到每個(gè)類中心的距離'''
for i in range(m):
for j in range(K):
dis[i,j] = np.dot((X[i,:]-initial_centroids[j,:]).reshape(1,-1),(X[i,:]-initial_centroids[j,:]).reshape(-1,1))
'''返回dis每一行的最小值對(duì)應(yīng)的列號(hào),即為對(duì)應(yīng)的類別
- np.min(dis, axis=1)返回每一行的最小值
- np.where(dis == np.min(dis, axis=1).reshape(-1,1)) 返回對(duì)應(yīng)最小值的坐標(biāo)
- 注意:可能最小值對(duì)應(yīng)的坐標(biāo)有多個(gè),where都會(huì)找出來,所以返回時(shí)返回前m個(gè)需要的即可(因?yàn)閷?duì)于多個(gè)最小值,屬于哪個(gè)類別都可以)
'''
dummy,idx = np.where(dis == np.min(dis, axis=1).reshape(-1,1))
return idx[0:dis.shape[0]] # 注意截取一下
計(jì)算類中心實(shí)現(xiàn)代碼:
# 計(jì)算類中心
def computerCentroids(X,idx,K):
n = X.shape[1]
centroids = np.zeros((K,n))
for i in range(K):
centroids[i,:] = np.mean(X[np.ravel(idx==i),:], axis=0).reshape(1,-1) # 索引要是一維的,axis=0為每一列,idx==i一次找出屬于哪一類的,然后計(jì)算均值
return centroids
2、目標(biāo)函數(shù)
也叫做失真代價(jià)函數(shù)
最后我們想得到:
其中表示第i條數(shù)據(jù)距離哪個(gè)類中心最近,
其中即為聚類的中心
3、聚類中心的選擇
隨機(jī)初始化,從給定的數(shù)據(jù)中隨機(jī)抽取K個(gè)作為聚類中心
隨機(jī)一次的結(jié)果可能不好,可以隨機(jī)多次,最后取使代價(jià)函數(shù)最小的作為中心
實(shí)現(xiàn)代碼:(這里隨機(jī)一次)
# 初始化類中心--隨機(jī)取K個(gè)點(diǎn)作為聚類中心
def kMeansInitCentroids(X,K):
m = X.shape[0]
m_arr = np.arange(0,m) # 生成0-m-1
centroids = np.zeros((K,X.shape[1]))
np.random.shuffle(m_arr) # 打亂m_arr順序
rand_indices = m_arr[:K] # 取前K個(gè)
centroids = X[rand_indices,:]
return centroids
4、聚類個(gè)數(shù)K的選擇
聚類是不知道y的label的,所以不知道真正的聚類個(gè)數(shù)
肘部法則(Elbow method)
作代價(jià)函數(shù)J和K的圖,若是出現(xiàn)一個(gè)拐點(diǎn),如下圖所示,K就取拐點(diǎn)處的值,下圖此時(shí)K=3
若是很平滑就不明確,人為選擇。
第二種就是人為觀察選擇
5、應(yīng)用——圖片壓縮
將圖片的像素分為若干類,然后用這個(gè)類代替原來的像素值
執(zhí)行聚類的算法代碼:
# 聚類算法
def runKMeans(X,initial_centroids,max_iters,plot_process):
m,n = X.shape # 數(shù)據(jù)條數(shù)和維度
K = initial_centroids.shape[0] # 類數(shù)
centroids = initial_centroids # 記錄當(dāng)前類中心
previous_centroids = centroids # 記錄上一次類中心
idx = np.zeros((m,1)) # 每條數(shù)據(jù)屬于哪個(gè)類
for i in range(max_iters): # 迭代次數(shù)
print u'迭代計(jì)算次數(shù):%d'%(i+1)
idx = findClosestCentroids(X, centroids)
if plot_process: # 如果繪制圖像
plt = plotProcessKMeans(X,centroids,previous_centroids) # 畫聚類中心的移動(dòng)過程
previous_centroids = centroids # 重置
centroids = computerCentroids(X, idx, K) # 重新計(jì)算類中心
if plot_process: # 顯示最終的繪制結(jié)果
plt.show()
return centroids,idx # 返回聚類中心和數(shù)據(jù)屬于哪個(gè)類
導(dǎo)入包
from sklearn.cluster import KMeans
使用模型擬合數(shù)據(jù)
model = KMeans(n_clusters=3).fit(X) # n_clusters指定3類,擬合數(shù)據(jù)
聚類中心
centroids = model.cluster_centers_ # 聚類中心
7、運(yùn)行結(jié)果
二維數(shù)據(jù)類中心的移動(dòng)
圖片壓縮
六、PCA主成分分析(降維)
1、用處
數(shù)據(jù)壓縮(Data Compression),使程序運(yùn)行更快
可視化數(shù)據(jù),例如3D-->2D等
......
2、2D-->1D,nD-->kD
如下圖所示,所有數(shù)據(jù)點(diǎn)可以投影到一條直線,是投影距離的平方和(投影誤差)最小
注意數(shù)據(jù)需要?dú)w一化處理
思路是找1個(gè)向量u,所有數(shù)據(jù)投影到上面使投影距離最小
那么nD-->kD就是找k個(gè)向量,所有數(shù)據(jù)投影到上面使投影誤差最小
eg:3D-->2D,2個(gè)向量就代表一個(gè)平面了,所有點(diǎn)投影到這個(gè)平面的投影誤差最小即可
3、主成分分析PCA與線性回歸的區(qū)別
線性回歸是找x與y的關(guān)系,然后用于預(yù)測y
PCA是找一個(gè)投影面,最小化data到這個(gè)投影面的投影誤差
4、PCA降維過程
數(shù)據(jù)預(yù)處理(均值歸一化)
公式:
就是減去對(duì)應(yīng)feature的均值,然后除以對(duì)應(yīng)特征的標(biāo)準(zhǔn)差(也可以是最大值-最小值)
實(shí)現(xiàn)代碼:
# 歸一化數(shù)據(jù)
def featureNormalize(X):
'''(每一個(gè)數(shù)據(jù)-當(dāng)前列的均值)/當(dāng)前列的標(biāo)準(zhǔn)差'''
n = X.shape[1]
mu = np.zeros((1,n));
sigma = np.zeros((1,n))
mu = np.mean(X,axis=0)
sigma = np.std(X,axis=0)
for i in range(n):
X[:,i] = (X[:,i]-mu[i])/sigma[i]
return X,mu,sigma
計(jì)算協(xié)方差矩陣Σ(Covariance Matrix):
注意這里的Σ和求和符號(hào)不同
協(xié)方差矩陣對(duì)稱正定(不理解正定的看看線代)
大小為nxn,n為feature的維度
實(shí)現(xiàn)代碼:
Sigma = np.dot(np.transpose(X_norm),X_norm)/m # 求Sigma
計(jì)算Σ的特征值和特征向量
可以是用svd奇異值分解函數(shù):U,S,V = svd(Σ)
返回的是與Σ同樣大小的對(duì)角陣S(由Σ的特征值組成)[注意:matlab中函數(shù)返回的是對(duì)角陣,在python中返回的是一個(gè)向量,節(jié)省空間]
還有兩個(gè)酉矩陣U和V,且
注意:svd函數(shù)求出的S是按特征值降序排列的,若不是使用svd,需要按特征值大小重新排列U
降維
選取U中的前K列(假設(shè)要降為K維)
Z就是對(duì)應(yīng)降維之后的數(shù)據(jù)
實(shí)現(xiàn)代碼:
# 映射數(shù)據(jù)
def projectData(X_norm,U,K):
Z = np.zeros((X_norm.shape[0],K))
U_reduce = U[:,0:K] # 取前K個(gè)
Z = np.dot(X_norm,U_reduce)
return Z
過程總結(jié):
Sigma = X'*X/m
U,S,V = svd(Sigma)
Ureduce = U[:,0:k]
Z = Ureduce'*x
5、數(shù)據(jù)恢復(fù)
因?yàn)?#xff1a;
所以: (注意這里是X的近似值)
又因?yàn)閁reduce為正定矩陣,【正定矩陣滿足:,所以:】,所以這里:
實(shí)現(xiàn)代碼:
# 恢復(fù)數(shù)據(jù)
def recoverData(Z,U,K):
X_rec = np.zeros((Z.shape[0],U.shape[0]))
U_recude = U[:,0:K]
X_rec = np.dot(Z,np.transpose(U_recude)) # 還原數(shù)據(jù)(近似)
return X_rec
6、主成分個(gè)數(shù)的選擇(即要降的維度)
如何選擇
投影誤差(project error):
總變差(total variation):
若誤差率(error ratio):,則稱99%保留差異性
誤差率一般取1%,5%,10%等
如何實(shí)現(xiàn)
若是一個(gè)個(gè)試的話代價(jià)太大
之前U,S,V = svd(Sigma),我們得到了S,這里誤差率error ratio:
可以一點(diǎn)點(diǎn)增加K嘗試。
7、使用建議
不要使用PCA去解決過擬合問題Overfitting,還是使用正則化的方法(如果保留了很高的差異性還是可以的)
只有在原數(shù)據(jù)上有好的結(jié)果,但是運(yùn)行很慢,才考慮使用PCA
8、運(yùn)行結(jié)果
2維數(shù)據(jù)降為1維
要投影的方向
2D降為1D及對(duì)應(yīng)關(guān)系
人臉數(shù)據(jù)降維
原始數(shù)據(jù)
可視化部分U矩陣信息
恢復(fù)數(shù)據(jù)
導(dǎo)入需要的包:
#-*- coding: utf-8 -*-
# Author:bob
# Date:2016.12.22
import numpy as np
from matplotlib import pyplot as plt
from scipy import io as spio
from sklearn.decomposition import pca
from sklearn.preprocessing import StandardScaler
歸一化數(shù)據(jù)
'''歸一化數(shù)據(jù)并作圖'''
scaler = StandardScaler()
scaler.fit(X)
x_train = scaler.transform(X)
使用PCA模型擬合數(shù)據(jù),并降維
n_components對(duì)應(yīng)要將的維度
'''擬合數(shù)據(jù)'''
K=1 # 要降的維度
model = pca.PCA(n_components=K).fit(x_train) # 擬合數(shù)據(jù),n_components定義要降的維度
Z = model.transform(x_train) # transform就會(huì)執(zhí)行降維操作
數(shù)據(jù)恢復(fù)
model.components_會(huì)得到降維使用的U矩陣
'''數(shù)據(jù)恢復(fù)并作圖'''
Ureduce = model.components_ # 得到降維用的Ureduce
x_rec = np.dot(Z,Ureduce) # 數(shù)據(jù)恢復(fù)
七、異常檢測 Anomaly Detection
1、高斯分布(正態(tài)分布)Gaussian distribution
分布函數(shù):
其中,u為數(shù)據(jù)的均值,σ為數(shù)據(jù)的標(biāo)準(zhǔn)差
σ越小,對(duì)應(yīng)的圖像越尖
參數(shù)估計(jì)(parameter estimation)
2、異常檢測算法
例子
訓(xùn)練集:,其中
假設(shè)相互獨(dú)立,建立model模型:
過程
選擇具有代表異常的feature:xi
參數(shù)估計(jì):
計(jì)算p(x),若是P(x)
這里只是單元高斯分布,假設(shè)了feature之間是獨(dú)立的,下面會(huì)講到多元高斯分布,會(huì)自動(dòng)捕捉到feature之間的關(guān)系
參數(shù)估計(jì)實(shí)現(xiàn)代碼
# 參數(shù)估計(jì)函數(shù)(就是求均值和方差)
def estimateGaussian(X):
m,n = X.shape
mu = np.zeros((n,1))
sigma2 = np.zeros((n,1))
mu = np.mean(X, axis=0) # axis=0表示列,每列的均值
sigma2 = np.var(X,axis=0) # 求每列的方差
return mu,sigma2
3、評(píng)價(jià)p(x)的好壞,以及ε的選取
對(duì)偏斜數(shù)據(jù)的錯(cuò)誤度量
因?yàn)閿?shù)據(jù)可能是非常偏斜的(就是y=1的個(gè)數(shù)非常少,(y=1表示異常)),所以可以使用Precision/Recall,計(jì)算F1Score(在CV交叉驗(yàn)證集上)
例如:預(yù)測癌癥,假設(shè)模型可以得到99%能夠預(yù)測正確,1%的錯(cuò)誤率,但是實(shí)際癌癥的概率很小,只有0.5%,那么我們始終預(yù)測沒有癌癥y=0反而可以得到更小的錯(cuò)誤率。使用error rate來評(píng)估就不科學(xué)了。
如下圖記錄:
,即:正確預(yù)測正樣本/所有預(yù)測正樣本
,即:正確預(yù)測正樣本/真實(shí)值為正樣本
總是讓y=1(較少的類),計(jì)算Precision和Recall
還是以癌癥預(yù)測為例,假設(shè)預(yù)測都是no-cancer,TN=199,FN=1,TP=0,FP=0,所以:Precision=0/0,Recall=0/1=0,盡管accuracy=199/200=99.5%,但是不可信。
ε的選取
嘗試多個(gè)ε值,使F1Score的值高
實(shí)現(xiàn)代碼
# 選擇最優(yōu)的epsilon,即:使F1Score最大
def selectThreshold(yval,pval):
'''初始化所需變量'''
bestEpsilon = 0.
bestF1 = 0.
F1 = 0.
step = (np.max(pval)-np.min(pval))/1000
'''計(jì)算'''
for epsilon in np.arange(np.min(pval),np.max(pval),step):
cvPrecision = pval
tp = np.sum((cvPrecision == 1) & (yval == 1).ravel()).astype(float) # sum求和是int型的,需要轉(zhuǎn)為float
fp = np.sum((cvPrecision == 1) & (yval == 0).ravel()).astype(float)
fn = np.sum((cvPrecision == 0) & (yval == 1).ravel()).astype(float)
precision = tp/(tp+fp) # 精準(zhǔn)度
recision = tp/(tp+fn) # 召回率
F1 = (2*precision*recision)/(precision+recision) # F1Score計(jì)算公式
if F1 > bestF1: # 修改最優(yōu)的F1 Score
bestF1 = F1
bestEpsilon = epsilon
return bestEpsilon,bestF1
4、選擇使用什么樣的feature(單元高斯分布)
如果一些數(shù)據(jù)不是滿足高斯分布的,可以變化一下數(shù)據(jù),例如log(x+C),x^(1/2)等
如果p(x)的值無論異常與否都很大,可以嘗試組合多個(gè)feature,(因?yàn)閒eature之間可能是有關(guān)系的)
5、多元高斯分布
單元高斯分布存在的問題
如下圖,紅色的點(diǎn)為異常點(diǎn),其他的都是正常點(diǎn)(比如CPU和memory的變化)
x1對(duì)應(yīng)的高斯分布如下:
x2對(duì)應(yīng)的高斯分布如下:
可以看出對(duì)應(yīng)的p(x1)和p(x2)的值變化并不大,就不會(huì)認(rèn)為異常
因?yàn)槲覀冋J(rèn)為feature之間是相互獨(dú)立的,所以如上圖是以正圓的方式擴(kuò)展
多元高斯分布
,并不是建立p(x1),p(x2)...p(xn),而是統(tǒng)一建立p(x)
其中參數(shù):,Σ為協(xié)方差矩陣
同樣,|Σ|越小,p(x)越尖
例如:
,
表示x1,x2正相關(guān),即x1越大,x2也就越大,如下圖,也就可以將紅色的異常點(diǎn)檢查出了
若:
,
表示x1,x2負(fù)相關(guān)
實(shí)現(xiàn)代碼:
# 多元高斯分布函數(shù)
def multivariateGaussian(X,mu,Sigma2):
k = len(mu)
if (Sigma2.shape[0]>1):
Sigma2 = np.diag(Sigma2)
'''多元高斯分布函數(shù)'''
X = X-mu
argu = (2*np.pi)**(-k/2)*np.linalg.det(Sigma2)**(-0.5)
p = argu*np.exp(-0.5*np.sum(np.dot(X,np.linalg.inv(Sigma2))*X,axis=1)) # axis表示每行
return p
6、單元和多元高斯分布特點(diǎn)
單元高斯分布
人為可以捕捉到feature之間的關(guān)系時(shí)可以使用
計(jì)算量小
多元高斯分布
自動(dòng)捕捉到相關(guān)的feature
計(jì)算量大,因?yàn)?#xff1a;
m>n或Σ可逆時(shí)可以使用。(若不可逆,可能有冗余的x,因?yàn)榫€性相關(guān),不可逆,或者就是m
7、程序運(yùn)行結(jié)果
顯示數(shù)據(jù)
等高線
異常點(diǎn)標(biāo)注
總結(jié)
以上是生活随笔為你收集整理的决策算法python_GitHub - nxety/MachineLearning_Python: 机器学习算法python实现的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 软件测试要求太高,软件“故障门”频现 对
- 下一篇: matlab如何实现降维,matlab怎