人生苦短，就用 Python。
在 Kaggle 最新發布的全球數據科學/機器學習現狀報告中，來自 50 多個國家的 16000 多位從業者紛紛向新手們推薦 Python 語言，用以學習機器學習。
那麼，用Python實現出來的機器學習算法都是什麼樣子呢？營長剛好在 GitHub 上發現了東南大學研究生「Lawlite」的一個項目——機器學習算法的Python實現，下面從線性回歸到反向傳播算法、從SVM到K-means聚類算法，咱們一一來分析其中的Python代碼。
目錄
一、線性回歸
二、邏輯回歸
邏輯回歸_手寫數字識別_OneVsAll
三、BP神經網絡
1、神經網絡model
2、代價函數
3、正則化
4、反向傳播BP
5、BP可以求梯度的原因
6、梯度檢查
7、權重的隨機初始化
8、預測
9、輸出結果
四、SVM支持向量機
1、代價函數
2、Large Margin
3、SVM Kernel（核函數）
4、使用中的模型代碼
5、運行結果
五、K-Means聚類算法
六、PCA主成分分析（降維）
七、異常檢測 Anomaly Detection
正文
一、線性回歸
https://github.com/lawlite19/MachineLearning_Python/tree/master/LinearRegression
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/LinearRegression/LinearRegression.py
1、代價函數
其中：
下面就是要求出theta，使代價最小，即代表我們擬合出來的方程距離真實值最近
共有m條數據，其中代表我們要擬合出來的方程到真實值距離的平方，平方的原因是因為可能有負值，正負可能會抵消
前面有係數2的原因是下面求梯度是對每個變量求偏導，2可以消去
實現代碼：
# 計算代價函數
def computerCost(X,y,theta):
m = len(y)
J = 0
J = (np.transpose(X*theta-y))*(X*theta-y)/(2*m) #計算代價J
return J
注意這裡的X是真實數據前加了一列1，因為有theta(0)
2、梯度下降算法
代價函數對求偏導得到：
所以對theta的更新可以寫為：
其中為學習速率，控制梯度下降的速度，一般取0.01,0.03,0.1,0.3
為什麼梯度下降可以逐步減小代價函數？
假設函數f(x)
泰勒展開：f(x+△x)=f(x)+f'(x)*△x+o(△x)，
令：△x=-α*f'(x) ,即負梯度方向乘以一個很小的步長α
將△x代入泰勒展開式中：f(x+x)=f(x)-α*[f'(x)]²+o(△x)
可以看出，α是取得很小的正數，[f'(x)]²也是正數，所以可以得出：f(x+△x)<=f(x)
所以沿著負梯度放下，函數在減小，多維情況一樣。
# 梯度下降算法
def gradientDescent(X,y,theta,alpha,num_iters):
m = len(y)
n = len(theta)
temp = np.matrix(np.zeros((n,num_iters))) # 暫存每次迭代計算的theta，轉化為矩陣形式
J_history = np.zeros((num_iters,1)) #記錄每次迭代計算的代價值
for i in range(num_iters): # 遍歷迭代次數
h = np.dot(X,theta) # 計算內積，matrix可以直接乘
temp[:,i] = theta - ((alpha/m)*(np.dot(np.transpose(X),h-y))) #梯度的計算
theta = temp[:,i]
J_history[i] = computerCost(X,y,theta) #調用計算代價函數
print '.',
return theta,J_history
3、均值歸一化
目的是使數據都縮放到一個範圍內，便於使用梯度下降算法
其中 為所有此feture數據的平均值
可以是最大值-最小值，也可以是這個feature對應的數據的標準差
實現代碼：
# 歸一化feature
def featureNormaliza(X):
X_norm = np.array(X) #將X轉化為numpy數組對象，才可以進行矩陣的運算
#定義所需變量
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) # 求每一列的標準差
for i in range(X.shape[1]): # 遍歷列
X_norm[:,i] = (X_norm[:,i]-mu[i])/sigma[i] # 歸一化
return X_norm,mu,sigma
注意預測的時候也需要均值歸一化數據
4、最終運行結果
代價隨迭代次數的變化
5、使用scikit-learn庫中的線性模型實現
https://github.com/lawlite19/MachineLearning_Python/blob/master/LinearRegression/LinearRegression_scikit-learn.py
導入包
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)
預測
#預測結果
result = model.predict(x_test)
二、邏輯回歸
https://github.com/lawlite19/MachineLearning_Python/tree/master/LogisticRegression
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/LogisticRegression/LogisticRegression.py
1、代價函數
可以綜合起來為：
其中：
為什麼不用線性回歸的代價函數表示，因為線性回歸的代價函數可能是非凸的，對於分類問題，使用梯度下降很難得到最小值，上面的代價函數是凸函數
的圖像如下，即y=1時：
可以看出，當趨於1，y=1,與預測值一致，此時付出的代價cost趨於0，若趨於0，y=1,此時的代價cost值非常大，我們最終的目的是最小化代價值
同理的圖像如下（y=0）：
2、梯度
同樣對代價函數求偏導：
可以看出與線性回歸的偏導數一致
推導過程
3、正則化
目的是為了防止過擬合
在代價函數中加上一項
注意j是重1開始的，因為theta(0)為一個常數項，X中最前面一列會加上1列1，所以乘積還是theta(0),feature沒有關係，沒有必要正則化
正則化後的代價：
# 代價函數
def costFunction(initial_theta,X,y,inital_lambda):
m = len(y)
J = 0
h = sigmoid(np.dot(X,initial_theta)) # 計算h(z)
theta1 = initial_theta.copy() # 因為正則化j=1從1開始，不包含0，所以複製一份，前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 # 正則化的代價方程
return 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))# 計算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型函數（即）
實現代碼：
# S型函數
def sigmoid(z):
h = np.zeros((len(z),1)) # 初始化，與z的長度一置
h = 1.0/(1.0+np.exp(-z)) return h
5、映射為多項式
因為數據的feture可能很少，導致偏差大，所以創造出一些feture結合
eg:映射為2次方的形式:
實現代碼：
# 映射為多項式
def mapFeature(X1,X2):
degree = 3; # 映射的最高次方
out = np.ones((X1.shape[0],1)) # 映射後的結果數組（取代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) #矩陣直接乘相當於matlab中的點乘.*
out = np.hstack((out, temp.reshape(-1,1)))
return out
6、使用scipy的優化方法
梯度下降使用scipy中optimize中的fmin_bfgs函數
調用scipy中的優化算法fmin_bfgs（擬牛頓法Broyden-Fletcher-Goldfarb-Shanno
costFunction是自己實現的一個求代價的函數，
initial_theta表示初始化的值,
fprime指定costFunction的梯度
args是其餘測參數，以元組的形式傳入，最後會將最小化costFunction的theta返回
result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, args=(X,y,initial_lambda))
7、運行結果
data1決策邊界和準確度
data2決策邊界和準確度
8、使用scikit-learn庫中的邏輯回歸模型實現
https://github.com/lawlite19/MachineLearning_Python/blob/master/LogisticRegression/LogisticRegression_scikit-learn.py
導入包
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.cross_validation import train_test_split
import numpy as np
劃分訓練集和測試集
# 劃分為訓練集和測試集
x_train,x_test,y_train,y_test = train_test_split(X,y,test_size=0.2)
歸一化
# 歸一化
scaler = StandardScaler()
scaler.fit(x_train)
x_train = scaler.fit_transform(x_train)
x_test = scaler.fit_transform(x_test)
邏輯回歸
#邏輯回歸
model = LogisticRegression()
model.fit(x_train,y_train)
預測
# 預測
predict = model.predict(x_test)
right = sum(predict == y_test)
predict = np.hstack((predict.reshape(-1,1),y_test.reshape(-1,1))) # 將預測值和真實值放在一塊，好觀察
print predict
print ('測試集準確率：%f%%'%(right*100.0/predict.shape[0])) #計算在測試集上的準確度
邏輯回歸_手寫數字識別_OneVsAll
https://github.com/lawlite19/MachineLearning_Python/blob/master/LogisticRegression
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/LogisticRegression/LogisticRegression_OneVsAll.py
1、隨機顯示100個數字
我沒有使用scikit-learn中的數據集，像素是20*20px，彩色圖如下
灰度圖：
實現代碼：
# 顯示100個數字
def display_data(imgData):
sum = 0
'''
顯示100個數（若是一個一個繪製將會非常慢，可以將要畫的數字整理好，放到一個矩陣中，顯示這個矩陣即可）
- 初始化一個二維數組
- 將每行的數據調整成圖像的矩陣，放進二維數組
- 顯示即可
'''
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指定以列優先，在matlab中是這樣的，python中需要指定，默認以行
sum += 1
plt.imshow(display_array,cmap='gray') #顯示灰度圖像
plt.axis('off')
plt.show()
2、OneVsAll
如何利用邏輯回歸解決多分類的問題，OneVsAll就是把當前某一類看成一類，其他所有類別看作一類，這樣有成了二分類的問題了
如下圖，把途中的數據分成三類，先把紅色的看成一類，把其他的看作另外一類，進行邏輯回歸，然後把藍色的看成一類，其他的再看成一類，以此類推...
可以看出大於2類的情況下，有多少類就要進行多少次的邏輯回歸分類
3、手寫數字識別
共有0-9，10個數字，需要10次分類
由於數據集y給出的是0,1,2...9的數字，而進行邏輯回歸需要0/1的label標記，所以需要對y處理
說一下數據集，前500個是0,500-1000是1,...,所以如下圖，處理後的y，前500行的第一列是1，其餘都是0,500-1000行第二列是1，其餘都是0....
然後調用梯度下降算法求解theta
實現代碼：
# 求每個分類的theta，最後返回所有的all_theta
def oneVsAll(X,y,num_labels,Lambda):
# 初始化變量
m,n = X.shape
all_theta = np.zeros((n+1,num_labels)) # 每一列對應相應分類的theta,共10列
X = np.hstack((np.ones((m,1)),X)) # X前補上一列1的偏置bias
class_y = np.zeros((m,num_labels)) # 數據的y對應0-9，需要映射為0/1的關係
initial_theta = np.zeros((n+1,1)) # 初始化一個分類的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=',')
'''遍歷每個分類，計算對應的theta值'''
for i in range(num_labels):
result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, args=(X,class_y[:,i],Lambda)) # 調用梯度下降的優化方法
all_theta[:,i] = result.reshape(1,-1) # 放入all_theta中
all_theta = np.transpose(all_theta)
return all_theta
4、預測
之前說過，預測的結果是一個概率值，利用學習出來的theta代入預測的S型函數中，每行的最大值就是是某個數字的最大概率，所在的列號就是預測的數字的真實值,因為在分類時，所有為0的將y映射在第一列，為1的映射在第二列，依次類推
實現代碼：
# 預測
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))) #預測
'''
返回h中每一行最大值所在的列號
- np.max(h, axis=1)返回h中每一行的最大值（是某個數字的最大概率）
- 最後where找到的最大概率所在的列號（列號即是對應的數字）
'''
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、運行結果
10次分類，在訓練集上的準確度：
6、使用scikit-learn庫中的邏輯回歸模型實現
https://github.com/lawlite19/MachineLearning_Python/blob/master/LogisticRegression/LogisticRegression_OneVsAll_scikit-learn.py
1、導入包
from scipy import io as spio
import numpy as np
from sklearn import svm
from sklearn.linear_model import LogisticRegression
2、加載數據
data = loadmat_data("data_digits.mat")
X = data['X'] # 獲取X數據，每一行對應一個數字20x20px
y = data['y'] # 這裡讀取mat文件y的shape=(5000, 1)
y = np.ravel(y) # 調用sklearn需要轉化成一維的(5000,)
3、擬合模型
model = LogisticRegression()
model.fit(X, y) # 擬合
4、預測
predict = model.predict(X) #預測
print u"預測準確度為：%f%%"%np.mean(np.float64(predict == y)*100)
5、輸出結果（在訓練集上的準確度）
三、BP神經網絡
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/NeuralNetwok/NeuralNetwork.py
1、神經網絡model
先介紹個三層的神經網絡，如下圖所示
輸入層（input layer）有三個units（為補上的bias，通常設為1）
表示第j層的第i個激勵，也稱為為單元unit
為第j層到第j+1層映射的權重矩陣，就是每條邊的權重
所以可以得到：
隱含層：
輸出層，
其中，S型函數，也成為激勵函數
可以看出 為3x4的矩陣，為1x4的矩陣
==》j+1的單元數x（j層的單元數+1）
2、代價函數
假設最後輸出的，即代表輸出層有K個單元
，
其中，代表第i個單元輸出與邏輯回歸的代價函數
差不多，就是累加上每個輸出（共有K個輸出）
3、正則化
L-->所有層的個數
-->第l層unit的個數
正則化後的代價函數為
共有L-1層，然後是累加對應每一層的theta矩陣，注意不包含加上偏置項對應的theta(0)
正則化後的代價函數實現代碼：
# 代價函數
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)) # 數據的y對應0-9，需要映射為0/1的關係
# 映射y
for i in range(num_labels):
class_y[:,i] = np.int32(y==i).reshape(1,-1) # 注意reshape(1,-1)才可以賦值
'''去掉theta1和theta2的第一列，因為正則化時從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))))
'''正向傳播,每次需要補上一列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)
'''代價'''
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(θ),使用梯度下降法還需要求它的梯度
BP反向傳播的目的就是求代價函數的梯度
假設4層的神經網絡,記為-->l層第j個單元的誤差
《===》（向量化）
沒有，因為對於輸入沒有誤差
因為S型函數的倒數為：
，
所以上面的和可以在前向傳播中計算出來
反向傳播計算梯度的過程為：
（是大寫的）
for i=1-m:-
-正向傳播計算（l=2,3,4...L）
-反向計算、...；
-
-
最後，即得到代價函數的梯度
實現代碼：
# 梯度
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)
Theta2 = nn_params[hidden_layer_size*(input_layer_size+1):length].reshape(num_labels,hidden_layer_size+1)
m = X.shape[0]
class_y = np.zeros((m,num_labels)) # 數據的y對應0-9，需要映射為0/1的關係
# 映射y
for i in range(num_labels):
class_y[:,i] = np.int32(y==i).reshape(1,-1) # 注意reshape(1,-1)才可以賦值
'''去掉theta1和theta2的第一列，因為正則化時從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)) #第一層到第二層的權重
Theta2_grad = np.zeros((Theta2.shape)) #第二層到第三層的權重
Theta1[:,0] = 0;
Theta2[:,0] = 0;
'''正向傳播，每次需要補上一列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,:]
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))
'''梯度'''
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可以求梯度的原因
實際是利用了鏈式求導法則
因為下一層的單元利用上一層的單元作為輸入進行計算
大體的推導過程如下，最終我們是想預測函數與已知的y非常接近，求均方差的梯度沿著此梯度方向可使代價函數最小化。可對照上面求梯度的過程。
求誤差更詳細的推導過程：
6、梯度檢查
檢查利用BP求的梯度是否正確
利用導數的定義驗證：
求出來的數值梯度應該與BP求出的梯度非常接近
驗證BP正確後就不需要再執行驗證梯度的算法了
實現代碼：
# 檢驗梯度是否計算正確
# 檢驗梯度是否計算正確
def checkGradient(Lambda = 0):
'''構造一個小型的神經網絡驗證，因為數值法計算梯度很浪費時間，而且驗證正確後之後就不再需要驗證了'''
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)
'''使用數值法計算梯度'''
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、權重的隨機初始化
神經網絡不能像邏輯回歸那樣初始化theta為0,因為若是每條邊的權重都為0，每個神經元都是相同的輸出，在反向傳播中也會得到同樣的梯度，最終只會預測一種結果。
所以應該初始化為接近0的數
實現代碼
# 隨機初始化權重theta
def randInitializeWeights(L_in,L_out):
W = np.zeros((L_out,1+L_in)) # 對應theta的權重
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)產生L_out*(1+L_in)大小的隨機矩陣
return W
8、預測
正向傳播預測結果
實現代碼
# 預測
def predict(Theta1,Theta2,X):
m = X.shape[0]
num_labels = Theta2.shape[0]
#p = np.zeros((m,1))
'''正向傳播，預測結果'''
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中每一行最大值所在的列號
- np.max(h, axis=1)返回h中每一行的最大值（是某個數字的最大概率）
- 最後where找到的最大概率所在的列號（列號即是對應的數字）
'''
#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、輸出結果
梯度檢查：
隨機顯示100個手寫數字
顯示theta1權重
訓練集預測準確度
歸一化後訓練集預測準確度
四、SVM支持向量機
1、代價函數
在邏輯回歸中，我們的代價為：
，
其中：
，
如圖所示，如果y=1，cost代價函數如圖所示
我們想讓，即z>>0，這樣的話cost代價函數才會趨於最小（這是我們想要的），所以用途中紅色的函數代替邏輯回歸中的cost
當y=0時同樣，用代替
最終得到的代價函數為：
最後我們想要
之前我們邏輯回歸中的代價函數為：
可以認為這裡的，只是表達形式問題，這裡C的值越大，SVM的決策邊界的margin也越大，下面會說明
2、Large Margin
如下圖所示,SVM分類會使用最大的margin將其分開
先說一下向量內積
，
表示u的歐幾裡得範數（歐式範數），
向量V在向量u上的投影的長度記為p，則：向量內積：
根據向量夾角公式推導一下即可，
前面說過，當C越大時，margin也就越大，我們的目的是最小化代價函數J(θ),當margin最大時，C的乘積項
要很小，所以近似為：
，
我們最後的目的就是求使代價最小的θ
由
可以得到：
，
p即為x在θ上的投影
如下圖所示，假設決策邊界如圖，找其中的一個點，到θ上的投影為p,則或者，若是p很小，則需要很大，這與我們要求的θ使最小相違背，所以最後求的是large margin
3、SVM Kernel（核函數）
對於線性可分的問題，使用線性核函數即可
對於線性不可分的問題，在邏輯回歸中，我們是將feature映射為使用多項式的形式，SVM中也有多項式核函數，但是更常用的是高斯核函數，也稱為RBF核
高斯核函數為：
假設如圖幾個點，
令：
，
. . .
可以看出，若是x與距離較近，==》，（即相似度較大），若是x與距離較遠，==》，（即相似度較低）
高斯核函數的σ越小，f下降的越快
如何選擇初始的
訓練集：
選擇：
對於給出的x，計算f,令：，
所以：
最小化J求出θ，
如果，==》預測y=1
4、使用scikit-learn中的SVM模型代碼
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/SVM/SVM_scikit-learn.py
線性可分的,指定核函數為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) # 指定核函數為線性核函數
非線性可分的，默認核函數為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為核函數的係數，值越大擬合的越好
5、運行結果
線性可分的決策邊界：
線性不可分的決策邊界：
五、K-Means聚類算法
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/K-Means/K-Menas.py
1、聚類過程
聚類屬於無監督學習，不知道y的標記分為K類
K-Means算法分為兩個步驟
第一步：簇分配，隨機選K個點作為中心，計算到這K個點的距離，分為K個簇
第二步：移動聚類中心：重新計算每個簇的中心，移動中心，重複以上步驟。
如下圖所示：
隨機分配的聚類中心
重新計算聚類中心，移動一次
最後10步之後的聚類中心
計算每條數據到哪個中心最近實現代碼：
# 找到每條數據距離哪個類中心最近
def findClosestCentroids(X,initial_centroids):
m = X.shape[0] # 數據條數
K = initial_centroids.shape[0] # 類的總數
dis = np.zeros((m,K)) # 存儲計算每個點分別到K個類的距離
idx = np.zeros((m,1)) # 要返回的每條數據屬於哪個類
'''計算每個點到每個類中心的距離'''
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每一行的最小值對應的列號，即為對應的類別
- np.min(dis, axis=1)返回每一行的最小值
- np.where(dis == np.min(dis, axis=1).reshape(-1,1)) 返回對應最小值的坐標
- 注意：可能最小值對應的坐標有多個，where都會找出來，所以返回時返回前m個需要的即可（因為對於多個最小值，屬於哪個類別都可以）
'''
dummy,idx = np.where(dis == np.min(dis, axis=1).reshape(-1,1))
return idx[0:dis.shape[0]] # 注意截取一下
計算類中心實現代碼：
# 計算類中心
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一次找出屬於哪一類的，然後計算均值
return centroids
2、目標函數
也叫做失真代價函數
最後我們想得到：
其中表示第i條數據距離哪個類中心最近，其中即為聚類的中心
3、聚類中心的選擇
隨機初始化，從給定的數據中隨機抽取K個作為聚類中心
隨機一次的結果可能不好，可以隨機多次，最後取使代價函數最小的作為中心
實現代碼：(這裡隨機一次)
# 初始化類中心--隨機取K個點作為聚類中心
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個
centroids = X[rand_indices,:]
return centroids
4、聚類個數K的選擇
聚類是不知道y的label的，所以不知道真正的聚類個數
肘部法則（Elbow method）
作代價函數J和K的圖，若是出現一個拐點，如下圖所示，K就取拐點處的值，下圖此時K=3
若是很平滑就不明確，人為選擇。
第二種就是人為觀察選擇
5、應用——圖片壓縮
將圖片的像素分為若干類，然後用這個類代替原來的像素值
執行聚類的算法代碼：
# 聚類算法
def runKMeans(X,initial_centroids,max_iters,plot_process):
m,n = X.shape # 數據條數和維度
K = initial_centroids.shape[0] # 類數
centroids = initial_centroids # 記錄當前類中心
previous_centroids = centroids # 記錄上一次類中心
idx = np.zeros((m,1)) # 每條數據屬於哪個類
for i in range(max_iters): # 迭代次數
print u'迭代計算次數：%d'%(i+1)
idx = findClosestCentroids(X, centroids)
if plot_process: # 如果繪製圖像
plt = plotProcessKMeans(X,centroids,previous_centroids) # 畫聚類中心的移動過程
previous_centroids = centroids # 重置
centroids = computerCentroids(X, idx, K) # 重新計算類中心
if plot_process: # 顯示最終的繪製結果
plt.show()
return centroids,idx # 返回聚類中心和數據屬於哪個類
6、使用scikit-learn庫中的線性模型實現聚類
https://github.com/lawlite19/MachineLearning_Python/blob/master/K-Means/K-Means_scikit-learn.py
導入包
from sklearn.cluster import KMeans
使用模型擬合數據
model = KMeans(n_clusters=3).fit(X) # n_clusters指定3類，擬合數據
聚類中心
centroids = model.cluster_centers_ # 聚類中心
7、運行結果
二維數據類中心的移動
圖片壓縮
六、PCA主成分分析（降維）
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/PCA/PCA.py
1、用處
數據壓縮（Data Compression）,使程序運行更快
可視化數據，例如3D-->2D等
.
2、2D-->1D，nD-->kD
如下圖所示，所有數據點可以投影到一條直線，是投影距離的平方和（投影誤差）最小
注意數據需要歸一化處理
思路是找1個向量u,所有數據投影到上面使投影距離最小
那麼nD-->kD就是找k個向量，
所有數據投影到上面使投影誤差最小
eg:3D-->2D,2個向量
就代表一個平面了，所有點投影到這個平面的投影誤差最小即可
3、主成分分析PCA與線性回歸的區別
線性回歸是找x與y的關係，然後用於預測y
PCA是找一個投影面，最小化data到這個投影面的投影誤差
4、PCA降維過程
數據預處理（均值歸一化）
公式：
就是減去對應feature的均值，然後除以對應特徵的標準差（也可以是最大值-最小值）
實現代碼：
# 歸一化數據
def featureNormalize(X):
'''（每一個數據-當前列的均值）/當前列的標準差'''
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
計算協方差矩陣Σ（Covariance Matrix）：
注意這裡的Σ和求和符號不同
協方差矩陣對稱正定（不理解正定的看看線代）
大小為nxn,n為feature的維度
實現代碼：
Sigma = np.dot(np.transpose(X_norm),X_norm)/m # 求Sigma
計算Σ的特徵值和特徵向量
可以是用svd奇異值分解函數：U,S,V = svd(Σ)
返回的是與Σ同樣大小的對角陣S（由Σ的特徵值組成）[注意：matlab中函數返回的是對角陣，在python中返回的是一個向量，節省空間]
還有兩個酉矩陣U和V，且
注意：svd函數求出的S是按特徵值降序排列的，若不是使用svd,需要按特徵值大小重新排列U
降維
選取U中的前K列（假設要降為K維）
Z就是對應降維之後的數據
實現代碼：
# 映射數據
def projectData(X_norm,U,K):
Z = np.zeros((X_norm.shape[0],K))
U_reduce = U[:,0:K] # 取前K個
Z = np.dot(X_norm,U_reduce)
return Z
過程總結：
Sigma = X'*X/m
U,S,V = svd(Sigma)
Ureduce = U[:,0:k]
Z = Ureduce'*x
5、數據恢復
因為：，
所以： （注意這裡是X的近似值）
又因為Ureduce為正定矩陣，【正定矩陣滿足：，所以：】，
所以這裡：
實現代碼：
# 恢復數據
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)) # 還原數據（近似）
return X_rec
6、主成分個數的選擇（即要降的維度）
如何選擇
投影誤差（project error）：
總變差（total variation）:
若誤差率（error ratio）：，則稱99%保留差異性
誤差率一般取1%，5%，10%等
如何實現
若是一個個試的話代價太大
之前U,S,V = svd(Sigma),我們得到了S，這裡誤差率error ratio:
可以一點點增加K嘗試。
7、使用建議
不要使用PCA去解決過擬合問題Overfitting，還是使用正則化的方法（如果保留了很高的差異性還是可以的）
只有在原數據上有好的結果，但是運行很慢，才考慮使用PCA
8、運行結果
2維數據降為1維
要投影的方向
2D降為1D及對應關係
人臉數據降維
原始數據
可視化部分U矩陣信息
恢復數據
9、使用scikit-learn庫中的PCA實現降維
https://github.com/lawlite19/MachineLearning_Python/blob/master/PCA/PCA.py_scikit-learn.py
導入需要的包：
#-*- 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
歸一化數據
'''歸一化數據並作圖'''
scaler = StandardScaler()
scaler.fit(X)
x_train = scaler.transform(X)
使用PCA模型擬合數據，並降維
n_components對應要將的維度
'''擬合數據'''
K=1 # 要降的維度
model = pca.PCA(n_components=K).fit(x_train) # 擬合數據，n_components定義要降的維度
Z = model.transform(x_train) # transform就會執行降維操作
數據恢復
model.components_會得到降維使用的U矩陣
'''數據恢復並作圖'''
Ureduce = model.components_ # 得到降維用的Ureduce
x_rec = np.dot(Z,Ureduce) # 數據恢復
七、異常檢測 Anomaly Detection
全部代碼
https://github.com/lawlite19/MachineLearning_Python/blob/master/AnomalyDetection/AnomalyDetection.py
1、高斯分布（正態分布）Gaussian distribution
分布函數：
其中，u為數據的均值，σ為數據的標準差
σ越小，對應的圖像越尖
參數估計（parameter estimation）
2、異常檢測算法
例子
訓練集：,其中
假設
相互獨立，建立model模型：
過程
選擇具有代表異常的feature:xi
參數估計：
計算p(x),若是P(x)<ε則認為異常，其中ε為我們要求的概率的臨界值threshold
這裡只是單元高斯分布，假設了feature之間是獨立的，下面會講到多元高斯分布，會自動捕捉到feature之間的關係
參數估計實現代碼
# 參數估計函數（就是求均值和方差）
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(x)的好壞，以及ε的選取
對偏斜數據的錯誤度量
因為數據可能是非常偏斜的（就是y=1的個數非常少，(y=1表示異常)），所以可以使用Precision/Recall，計算F1Score(在CV交叉驗證集上)
例如：預測癌症，假設模型可以得到99%能夠預測正確，1%的錯誤率，但是實際癌症的概率很小，只有0.5%，那麼我們始終預測沒有癌症y=0反而可以得到更小的錯誤率。使用error rate來評估就不科學了。
如下圖記錄：
，即：正確預測正樣本/所有預測正樣本
，即：正確預測正樣本/真實值為正樣本
總是讓y=1(較少的類)，計算Precision和Recall
還是以癌症預測為例，假設預測都是no-cancer，TN=199，FN=1，TP=0，FP=0，所以：Precision=0/0，Recall=0/1=0，儘管accuracy=199/200=99.5%，但是不可信。
ε的選取
嘗試多個ε值，使F1Score的值高
實現代碼
# 選擇最優的epsilon，即：使F1Score最大
def selectThreshold(yval,pval):
'''初始化所需變量'''
bestEpsilon = 0.
bestF1 = 0.
F1 = 0.
step = (np.max(pval)-np.min(pval))/1000
'''計算'''
for epsilon in np.arange(np.min(pval),np.max(pval),step):
cvPrecision = pval<epsilon
tp = np.sum((cvPrecision == 1) & (yval == 1)).astype(float) # sum求和是int型的，需要轉為float
fp = np.sum((cvPrecision == 1) & (yval == 0)).astype(float)
fn = np.sum((cvPrecision == 1) & (yval == 0)).astype(float)
precision = tp/(tp+fp) # 精準度
recision = tp/(tp+fn) # 召回率
F1 = (2*precision*recision)/(precision+recision) # F1Score計算公式
if F1 > bestF1: # 修改最優的F1 Score
bestF1 = F1
bestEpsilon = epsilon
return bestEpsilon,bestF1
4、選擇使用什麼樣的feature（單元高斯分布）
如果一些數據不是滿足高斯分布的，可以變化一下數據，例如log(x+C),x^(1/2)等
如果p(x)的值無論異常與否都很大，可以嘗試組合多個feature,(因為feature之間可能是有關係的)
5、多元高斯分布
單元高斯分布存在的問題
如下圖，紅色的點為異常點，其他的都是正常點（比如CPU和memory的變化）
x1對應的高斯分布如下：
x2對應的高斯分布如下：
可以看出對應的p(x1)和p(x2)的值變化並不大，就不會認為異常
因為我們認為feature之間是相互獨立的，所以如上圖是以正圓的方式擴展
多元高斯分布
，並不是建立p(x1),p(x2)...p(xn)，而是統一建立p(x)
其中參數：,Σ為協方差矩陣
同樣，|Σ|越小，p(x)越尖
例如：
，
表示x1,x2正相關，即x1越大，x2也就越大，如下圖，也就可以將紅色的異常點檢查出了
若：
表示x1,x2負相關
實現代碼：
# 多元高斯分布函數
def multivariateGaussian(X,mu,Sigma2):
k = len(mu)
if (Sigma2.shape[0]>1):
Sigma2 = np.diag(Sigma2)
'''多元高斯分布函數'''
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、單元和多元高斯分布特點
單元高斯分布
人為可以捕捉到feature之間的關係時可以使用
計算量小
多元高斯分布
自動捕捉到相關的feature
計算量大，因為：
m>n或Σ可逆時可以使用。（若不可逆，可能有冗餘的x，因為線性相關，不可逆，或者就是m<n）
7、程序運行結果
顯示數據
等高線
異常點標註
原文地址
https://github.com/lawlite19/MachineLearning_Python#
作者 | Lawlite
西普教育一直深耕於高校教育領域，依託自身十五年的高校教育經驗和強大的研發創新能力，基於大數據、人工智慧等實驗教學產品，為高校提供教學、科研相關專業的實驗教學解決方案，旨在為大數據人工智慧領域提供實戰、複合型人才。
西普教育
微信號 : xipujiaoyu
新浪微博：@西普教育
現已入駐：今日頭條、百度百家、搜狐號、一點資訊，歡迎訂閱