機(jī)器學(xué)習(xí)算法Python實(shí)現(xiàn)

目錄

• 機(jī)器學(xué)習(xí)算法Python實(shí)現(xiàn)
  • 邏輯回歸_手寫(xiě)數(shù)字識(shí)別_OneVsAll
  • 六、PCA主成分分析（降維）
    • 3、主成分分析PCA與線(xiàn)性回歸的區(qū)別
    • 6、主成分個(gè)數(shù)的選擇（即要降的維度）
    • 9、使用scikit-learn庫(kù)中的PCA實(shí)現(xiàn)降維
  • 七、異常檢測(cè) Anomaly Detection
    • 1、高斯分布（正態(tài)分布）
    • 3、評(píng)價(jià)的好壞，以及的選取
    • 4、選擇使用什么樣的feature（單元高斯分布）
    • 6、單元和多元高斯分布特點(diǎn)

一、 線(xiàn)性回歸

• 全部代碼

1、代價(jià)函數(shù)

•  

• 其中：

• 下面就是要求出theta，使代價(jià)最小，即代表我們擬合出來(lái)的方程距離真實(shí)值最近

• 共有m條數(shù)據(jù)，其中 代表我們要擬合出來(lái)的方程到真實(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的更新可以寫(xiě)為：
• 其中 為學(xué)習(xí)速率，控制梯度下降的速度，一般取 0.01,0.03,0.1,0.3.....
• 為什么梯度下降可以逐步減小代價(jià)函數(shù)
• 假設(shè)函數(shù) f(x)
• 泰勒展開(kāi)： f(x+△x)=f(x)+f'(x)*△x+o(△x)
• 令： △x=-α*f'(x) ,即負(fù)梯度方向乘以一個(gè)很小的步長(zhǎng) α
• 將 △x 代入泰勒展開(kāi)式中： 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ù)測(cè)的時(shí)候也需要均值歸一化數(shù)據(jù)

4、最終運(yùn)行結(jié)果

• 代價(jià)隨迭代次數(shù)的變化

  

5、 使用scikit-learn庫(kù)中的線(xiàn)性模型實(shí)現(xiàn)

• 導(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]))

• 線(xiàn)性模型擬合

# 線(xiàn)性模型擬合
    model = linear_model.LinearRegression()
    model.fit(x_train, y)

• 預(yù)測(cè)

#預(yù)測(cè)結(jié)果
    result = model.predict(x_test)

二、 邏輯回歸

• 全部代碼

1、代價(jià)函數(shù)

•  
• 可以綜合起來(lái)為： 其中：
• 為什么不用線(xiàn)性回歸的代價(jià)函數(shù)表示，因?yàn)榫�(xiàn)性回歸的代價(jià)函數(shù)可能是非凸的，對(duì)于分類(lèi)問(wèn)題，使用梯度下降很難得到最小值，上面的代價(jià)函數(shù)是凸函數(shù)
• 的圖像如下，即 y=1 時(shí)：

可以看出，當(dāng) 趨于 1 ， y=1 ,與預(yù)測(cè)值一致，此時(shí)付出的代價(jià) cost 趨于 0 ，若 趨于 0 ， y=1 ,此時(shí)的代價(jià) cost 值非常大，我們最終的目的是最小化代價(jià)值

• 同理 的圖像如下（ y=0 ）：

  

2、梯度

• 同樣對(duì)代價(jià)函數(shù)求偏導(dǎo)：
  可以看出與線(xiàn)性回歸的偏導(dǎo)數(shù)一致
• 推到過(guò)程

3、正則化

• 目的是為了防止過(guò)擬合
• 在代價(jià)函數(shù)中加上一項(xiàng)
• 注意j是重1開(kāi)始的，因?yàn)閠heta(0)為一個(gè)常數(shù)項(xiàng)，X中最前面一列會(huì)加上1列1，所以乘積還是theta(0),feature沒(méi)有關(guān)系，沒(méi)有必要正則化
• 正則化后的代價(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開(kāi)始，不包含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的長(zhǎng)度一置
    
    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是其余測(cè)參數(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)確度
  

8、 使用scikit-learn庫(kù)中的邏輯回歸模型實(shí)現(xià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)練集和測(cè)試集

# 劃分為訓(xùn)練集和測(cè)試集
    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)

• 預(yù)測(cè)

# 預(yù)測(cè)
    predict = model.predict(x_test)
    right = sum(predict == y_test)
    
    predict = np.hstack((predict.reshape(-1,1),y_test.reshape(-1,1)))   # 將預(yù)測(cè)值和真實(shí)值放在一塊，好觀(guān)察
    print predict
    print ('測(cè)試集準(zhǔn)確率：%f%%'%(right*100.0/predict.shape[0]))          #計(jì)算在測(cè)試集上的準(zhǔn)確度

邏輯回歸_手寫(xiě)數(shù)字識(shí)別_OneVsAll

• 全部代碼

1、隨機(jī)顯示100個(gè)數(shù)字

• 我沒(méi)有使用scikit-learn中的數(shù)據(jù)集，像素是20*20px，彩色圖如下 灰度圖：
• 實(shí)現(xiàn)代碼：

# 顯示100個(gè)數(shù)字
def display_data(imgData):
    sum = 0
    '''
    顯示100個(gè)數(shù)（若是一個(gè)一個(gè)繪制將會(huì)非常慢，可以將要畫(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

• 如何利用邏輯回歸解決多分類(lèi)的問(wèn)題，OneVsAll就是把當(dāng)前某一類(lèi)看成一類(lèi)，其他所有類(lèi)別看作一類(lèi)，這樣有成了二分類(lèi)的問(wèn)題了
• 如下圖，把途中的數(shù)據(jù)分成三類(lèi)，先把紅色的看成一類(lèi)，把其他的看作另外一類(lèi)，進(jìn)行邏輯回歸，然后把藍(lán)色的看成一類(lèi)，其他的再看成一類(lèi)，以此類(lèi)推...
• 可以看出大于2類(lèi)的情況下，有多少類(lèi)就要進(jìn)行多少次的邏輯回歸分類(lèi)

3、手寫(xiě)數(shù)字識(shí)別

• 共有0-9，10個(gè)數(shù)字，需要10次分類(lèi)
• 由于 數(shù)據(jù)集y 給出的是 0,1,2...9 的數(shù)字，而進(jìn)行邏輯回歸需要 0/1 的label標(biāo)記，所以需要對(duì)y處理
• 說(shuō)一下數(shù)據(jù)集，前 500 個(gè)是 0 , 500-1000 是 1 , ... ,所以如下圖，處理后的 y ， 前500行的第一列是1，其余都是0,500-1000行第二列是1，其余都是0....
• 然后調(diào)用 梯度下降算法 求解 theta
• 實(shí)現(xiàn)代碼：

# 求每個(gè)分類(lèi)的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)分類(lèi)的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è)分類(lèi)的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è)分類(lèi)，計(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ù)測(cè)

• 之前說(shuō)過(guò)，預(yù)測(cè)的結(jié)果是一個(gè) 概率值 ，利用學(xué)習(xí)出來(lái)的 theta 代入預(yù)測(cè)的 S型函數(shù) 中，每行的最大值就是是某個(gè)數(shù)字的最大概率，所在的 列號(hào) 就是預(yù)測(cè)的數(shù)字的真實(shí)值,因?yàn)樵诜诸?lèi)時(shí)，所有為 0 的將 y 映射在第一列，為1的映射在第二列，依次類(lèi)推
• 實(shí)現(xiàn)代碼：

# 預(yù)測(cè)
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ù)測(cè)

    '''
    返回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次分類(lèi)，在訓(xùn)練集上的準(zhǔn)確度：

  

6、 使用scikit-learn庫(kù)中的邏輯回歸模型實(shí)現(xià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ù)測(cè)

predict = model.predict(X) #預(yù)測(cè)
    
    print u"預(yù)測(cè)準(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ì)，也稱(chēng)為為單元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的中長(zhǎng)度
    # 還原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開(kāi)始'''    
    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è)單元的誤差

• 《===》 （向量化）

•  

•  

• 沒(méi)有 ，因?yàn)閷?duì)于輸入沒(méi)有誤差

• 因?yàn)镾型函數(shù) 的導(dǎo)數(shù)為： ，所以上面的 和 可以在前向傳播中計(jì)算出來(lái)

• 反向傳播計(jì)算梯度的過(guò)程為：

• （ 是大寫(xiě)的 ）

• 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開(kāi)始'''
    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)橄乱粚拥膯卧�利用上一層的單元作為輸入進(jìn)行計(jì)算
• 大體的推導(dǎo)過(guò)程如下，最終我們是想預(yù)測(cè)函數(shù)與已知的 y 非常接近，求均方差的梯度沿著此梯度方向可使代價(jià)函數(shù)最小化�？蓪�(duì)照上面求梯度的過(guò)程。
• 求誤差更詳細(xì)的推導(dǎo)過(guò)程：

6、梯度檢查

• 檢查利用 BP 求的梯度是否正確
• 利用導(dǎo)數(shù)的定義驗(yàn)證：
• 求出來(lái)的數(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)))  #展開(kāi)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ù)測(cè)一種結(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ù)測(cè)

• 正向傳播預(yù)測(cè)結(jié)果
• 實(shí)現(xiàn)代碼

# 預(yù)測(cè)
def predict(Theta1,Theta2,X):
    m = X.shape[0]
    num_labels = Theta2.shape[0]
    #p = np.zeros((m,1))
    '''正向傳播，預(yù)測(cè)結(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è)手寫(xiě)數(shù)字
  
• 顯示theta1權(quán)重
  
• 訓(xùn)練集預(yù)測(cè)準(zhǔn)確度
  
• 歸一化后訓(xùn)練集預(yù)測(cè)準(zhǔn)確度
  

四、SVM支持向量機(jī)

1、代價(jià)函數(shù)

• 在邏輯回歸中，我們的代價(jià)為：
  ，
  其中： ，
• 如圖所示，如果 y=1 ， cost 代價(jià)函數(shù)如圖所示
  
  我們想讓 ，即 z>>0 ，這樣的話(huà) cost 代價(jià)函數(shù)才會(huì)趨于最�。ㄟ@是我們想要的），所以用途中 紅色 的函數(shù) 代替邏輯回歸中的cost
• 當(dāng) y=0 時(shí)同樣，用 代替
• 最終得到的代價(jià)函數(shù)為：
  
  最后我們想要
• 之前我們邏輯回歸中的代價(jià)函數(shù)為：
  
  可以認(rèn)為這里的 ，只是表達(dá)形式問(wèn)題，這里 C 的值越大，SVM的決策邊界的 margin 也越大，下面會(huì)說(shuō)明

2、Large Margin

• 如下圖所示,SVM分類(lèi)會(huì)使用最大的 margin 將其分開(kāi)

  

• 先說(shuō)一下向量?jī)?nèi)積

• ，

• 表示 u 的 歐幾里得范數(shù) （歐式范數(shù)），

• 向量V 在 向量u 上的投影的長(zhǎng)度記為 p ，則：向量?jī)?nèi)積：

   

  

  根據(jù)向量夾角公式推導(dǎo)一下即可，

• 前面說(shuō)過(guò)，當(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ì)于線(xiàn)性可分的問(wèn)題，使用 線(xiàn)性核函數(shù) 即可

• 對(duì)于線(xiàn)性不可分的問(wèn)題，在邏輯回歸中，我們是將 feature 映射為使用多項(xiàng)式的形式 ， SVM 中也有 多項(xiàng)式核函數(shù) ，但是更常用的是 高斯核函數(shù) ，也稱(chēng)為 RBF核

• 高斯核函數(shù)為：

  假設(shè)如圖幾個(gè)點(diǎn)， 令：

   

  . . .

• 可以看出，若是 x 與 距離較近，==》 ，（即相似度較大）

  若是 x 與 距離較遠(yuǎn)，==》 ，（即相似度較低）

• 高斯核函數(shù)的 σ 越小， f 下降的越快

  

• 如何選擇初始的

• 訓(xùn)練集：

• 選擇：

• 對(duì)于給出的 x ，計(jì)算 f ,令： 所以：

• 最小化 J 求出 θ ，

• 如果 ，==》預(yù)測(cè) y=1

4、使用 scikit-learn 中的 SVM 模型代碼

• 全部代碼
• 線(xiàn)性可分的,指定核函數(shù)為 linear ：

'''data1——線(xiàn)性分類(lèi)'''
    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ù)為線(xiàn)性核函數(shù)

• 非線(xiàn)性可分的，默認(rèn)核函數(shù)為 rbf

'''data2——非線(xiàn)性分類(lèi)'''
    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é)果

• 線(xiàn)性可分的決策邊界：
  
• 線(xiàn)性不可分的決策邊界：
  

五、K-Means聚類(lèi)算法

• 全部代碼

1、聚類(lèi)過(guò)程

• 聚類(lèi)屬于無(wú)監(jiān)督學(xué)習(xí)，不知道y的標(biāo)記分為K類(lèi)

• K-Means算法分為兩個(gè)步驟

• 第一步：簇分配，隨機(jī)選 K 個(gè)點(diǎn)作為中心，計(jì)算到這 K 個(gè)點(diǎn)的距離，分為 K 個(gè)簇

• 第二步：移動(dòng)聚類(lèi)中心：重新計(jì)算每個(gè) 簇 的中心，移動(dòng)中心，重復(fù)以上步驟。

• 如下圖所示：

• 隨機(jī)分配的聚類(lèi)中心

  

• 重新計(jì)算聚類(lèi)中心，移動(dòng)一次

  

• 最后 10 步之后的聚類(lèi)中心

  

• 計(jì)算每條數(shù)據(jù)到哪個(gè)中心最近實(shí)現(xiàn)代碼：

# 找到每條數(shù)據(jù)距離哪個(gè)類(lèi)中心最近    
def findClosestCentroids(X,initial_centroids):
    m = X.shape[0]                  # 數(shù)據(jù)條數(shù)
    K = initial_centroids.shape[0]  # 類(lèi)的總數(shù)
    dis = np.zeros((m,K))           # 存儲(chǔ)計(jì)算每個(gè)點(diǎn)分別到K個(gè)類(lèi)的距離
    idx = np.zeros((m,1))           # 要返回的每條數(shù)據(jù)屬于哪個(gè)類(lèi)
    
    '''計(jì)算每個(gè)點(diǎn)到每個(gè)類(lèi)中心的距離'''
    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)的類(lèi)別
    - 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ì)找出來(lái)，所以返回時(shí)返回前m個(gè)需要的即可（因?yàn)閷?duì)于多個(gè)最小值，屬于哪個(gè)類(lèi)別都可以）
    '''  
    dummy,idx = np.where(dis == np.min(dis, axis=1).reshape(-1,1))
    return idx[0:dis.shape[0]]  # 注意截取一下

• 計(jì)算類(lèi)中心實(shí)現(xiàn)代碼：

# 計(jì)算類(lèi)中心
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一次找出屬于哪一類(lèi)的，然后計(jì)算均值
    return centroids

2、目標(biāo)函數(shù)

• 也叫做 失真代價(jià)函數(shù)
•  
• 最后我們想得到：
  
• 其中 表示第 i 條數(shù)據(jù)距離哪個(gè)類(lèi)中心最近，
• 其中 即為聚類(lèi)的中心

3、聚類(lèi)中心的選擇

• 隨機(jī)初始化，從給定的數(shù)據(jù)中隨機(jī)抽取K個(gè)作為聚類(lèi)中心
• 隨機(jī)一次的結(jié)果可能不好，可以隨機(jī)多次，最后取使代價(jià)函數(shù)最小的作為中心
• 實(shí)現(xiàn)代碼：(這里隨機(jī)一次)

# 初始化類(lèi)中心--隨機(jī)取K個(gè)點(diǎn)作為聚類(lèi)中心
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、聚類(lèi)個(gè)數(shù)K的選擇

• 聚類(lèi)是不知道y的label的，所以不知道真正的聚類(lèi)個(gè)數(shù)
• 肘部法則（Elbow method）
• 作代價(jià)函數(shù) J 和 K 的圖，若是出現(xiàn)一個(gè)拐點(diǎn)，如下圖所示， K 就取拐點(diǎn)處的值，下圖此時(shí) K=3
• 若是很平滑就不明確，人為選擇。
• 第二種就是人為觀(guān)察選擇

5、應(yīng)用——圖片壓縮

• 將圖片的像素分為若干類(lèi)，然后用這個(gè)類(lèi)代替原來(lái)的像素值
• 執(zhí)行聚類(lèi)的算法代碼：

# 聚類(lèi)算法
def runKMeans(X,initial_centroids,max_iters,plot_process):
    m,n = X.shape                   # 數(shù)據(jù)條數(shù)和維度
    K = initial_centroids.shape[0]  # 類(lèi)數(shù)
    centroids = initial_centroids   # 記錄當(dāng)前類(lèi)中心
    previous_centroids = centroids  # 記錄上一次類(lèi)中心
    idx = np.zeros((m,1))           # 每條數(shù)據(jù)屬于哪個(gè)類(lèi)
    
    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) # 畫(huà)聚類(lèi)中心的移動(dòng)過(guò)程
            previous_centroids = centroids  # 重置
        centroids = computerCentroids(X, idx, K)    # 重新計(jì)算類(lèi)中心
    if plot_process:    # 顯示最終的繪制結(jié)果
        plt.show()
    return centroids,idx    # 返回聚類(lèi)中心和數(shù)據(jù)屬于哪個(gè)類(lèi)

6、 使用scikit-learn庫(kù)中的線(xiàn)性模型實(shí)現(xiàn)聚類(lèi)

• 導(dǎo)入包

from sklearn.cluster import KMeans

• 使用模型擬合數(shù)據(jù)

model = KMeans(n_clusters=3).fit(X) # n_clusters指定3類(lèi)，擬合數(shù)據(jù)

• 聚類(lèi)中心

centroids = model.cluster_centers_  # 聚類(lèi)中心

7、運(yùn)行結(jié)果

• 二維數(shù)據(jù)類(lèi)中心的移動(dòng)
  
• 圖片壓縮
  

六、PCA主成分分析（降維）

• 全部代碼

1、用處

• 數(shù)據(jù)壓縮（Data Compression）,使程序運(yùn)行更快
• 可視化數(shù)據(jù)，例如 3D-->2D 等
• ......

2、2D-->1D，nD-->kD

• 如下圖所示，所有數(shù)據(jù)點(diǎn)可以投影到一條直線(xiàn)，是 投影距離的平方和 （投影誤差）最小
• 注意數(shù)據(jù)需要 歸一化 處理
• 思路是找 1 個(gè) 向量u ,所有數(shù)據(jù)投影到上面使投影距離最小
• 那么 nD-->kD 就是找 k 個(gè)向量 ，所有數(shù)據(jù)投影到上面使投影誤差最小
• eg:3D-->2D,2個(gè)向量 就代表一個(gè)平面了，所有點(diǎn)投影到這個(gè)平面的投影誤差最小即可

3、主成分分析PCA與線(xiàn)性回歸的區(qū)別

• 線(xiàn)性回歸是找 x 與 y 的關(guān)系，然后用于預(yù)測(cè) y
• PCA 是找一個(gè)投影面，最小化data到這個(gè)投影面的投影誤差

4、PCA降維過(guò)程

• 數(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ì)稱(chēng)正定 （不理解正定的看看線(xiàn)代）
• 大小為 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

• 過(guò)程總結(jié)：
• Sigma = X'*X/m
• U,S,V = svd(Sigma)
• Ureduce = U[:,0:k]
• Z = Ureduce'*x

5、數(shù)據(jù)恢復(fù)

• 因?yàn)椋?/li>
• 所以： （注意這里是X的近似值）
• 又因?yàn)?Ureduce 為正定矩陣，【正定矩陣滿(mǎn)足： ，所以： 】，所以這里：
•  
• 實(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）： ，則稱(chēng) 99% 保留差異性
• 誤差率一般取 1%，5%，10% 等
• 如何實(shí)現(xiàn)
• 若是一個(gè)個(gè)試的話(huà)代價(jià)太大
• 之前 U,S,V = svd(Sigma) ,我們得到了 S ，這里誤差率error ratio:
• 可以一點(diǎn)點(diǎn)增加 K 嘗試。

7、使用建議

• 不要使用PCA去解決過(guò)擬合問(wèn)題 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ù)
  

9、 使用scikit-learn庫(kù)中的PCA實(shí)現(xiàn)降維

• 導(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ù)

七、異常檢測(cè) 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、異常檢測(cè)算法

• 例子
• 訓(xùn)練集： ,其中
• 假設(shè) 相互獨(dú)立，建立model模型：
• 過(guò)程
• 選擇具有代表異常的 feature :xi
• 參數(shù)估計(jì)：
• 計(jì)算 p(x) ,若是 P(x)<ε 則認(rèn)為異常，其中 ε 為我們要求的概率的臨界值 threshold
• 這里只是 單元高斯分布 ，假設(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ù)測(cè)癌癥，假設(shè)模型可以得到 99% 能夠預(yù)測(cè)正確， 1% 的錯(cuò)誤率，但是實(shí)際癌癥的概率很小，只有 0.5% ，那么我們始終預(yù)測(cè)沒(méi)有癌癥y=0反而可以得到更小的錯(cuò)誤率。使用 error rate 來(lái)評(píng)估就不科學(xué)了。

• 如下圖記錄：

  

• ，即： 正確預(yù)測(cè)正樣本/所有預(yù)測(cè)正樣本

• ，即： 正確預(yù)測(cè)正樣本/真實(shí)值為正樣本

• 總是讓 y=1 (較少的類(lèi))，計(jì)算 Precision 和 Recall

• 

• 還是以癌癥預(yù)測(cè)為例，假設(shè)預(yù)測(cè)都是no-cancer，TN=199，F(xiàn)N=1，TP=0，F(xiàn)P=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<epsilon
        tp = np.sum((cvPrecision == 1) & (yval == 1)).astype(float)  # sum求和是int型的，需要轉(zhuǎn)為float
        fp = np.sum((cvPrecision == 1) & (yval == 0)).astype(float)
        fn = np.sum((cvPrecision == 1) & (yval == 0)).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ù)不是滿(mǎn)足高斯分布的，可以變化一下數(shù)據(jù)，例如 log(x+C),x^(1/2) 等
• 如果 p(x) 的值無(wú)論異常與否都很大，可以嘗試組合多個(gè) feature ,(因?yàn)閒eature之間可能是有關(guān)系的)

5、多元高斯分布

• 單元高斯分布存在的問(wèn)題
• 如下圖，紅色的點(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)椋?/li>
• m>n 或 Σ 可逆時(shí)可以使用。（若不可逆，可能有冗余的x，因?yàn)榫�(xiàn)性相關(guān)，不可逆，或者就是m<n）

7、程序運(yùn)行結(jié)果

• 顯示數(shù)據(jù)
  
• 等高線(xiàn)
  
• 異常點(diǎn)標(biāo)注
  

 

 

標(biāo)簽： isp 代碼 網(wǎng)絡(luò)

版權(quán)申明：本站文章部分自網(wǎng)絡(luò)，如有侵權(quán)，請(qǐng)聯(lián)系：west999com@outlook.com
特別注意：本站所有轉(zhuǎn)載文章言論不代表本站觀(guān)點(diǎn)！
本站所提供的圖片等素材，版權(quán)歸原作者所有，如需使用，請(qǐng)與原作者聯(lián)系。