機器學習筆記（三）：線性迴歸大解剖（程式碼部分）
阿新 • • 發佈：2018-12-12
這裡，讓我手把手教你如何用邏輯迴歸分析資料
根據學生分數預測是否錄取：
#必備3個庫
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
讓我們讀入資料：
import os path = "data" + os.sep + "LogiReg_data.txt"; # os.sep 表示的是檔案分割符 “ / ” 因為win和linux分隔符不同所以 為了跨平臺就不直接寫 pdData = pd.read_csv(path,header=None,names=['Exam1','Exam2','Admitted']); #header預設使用第一行做標題行，這裡改為自己設定 pdData.shape; #資料的緯度 (100 , 3) pdData.head() #不加引數表示預設顯示5行
前兩列資料表示分數，最後一列表示是否錄取：
讓我們簡單看一下，錄取和沒錄取在不同分數上的顯示情況
positive = pdData[pdData['Admitted' ] == 1] # 錄取的集合 nagative = pdData[pdData['Admitted' ] == 0] # 沒錄取集合 fig, ax=plt.subplots(figsize=(10,5)) # ax前面必須有fig figsize設定子圖的寬和高 #分1 做x 分2做y ax.scatter(positive['Exam1'],positive['Exam2'],s=30,c='blue',marker='o',label='Admitted') #s=點大小 c='顏色' marker=標記 ax.scatter(nagative['Exam1'],nagative['Exam2'],s=30,c='r',marker='x',label='Not Admitted') ax.legend() ax.set_xlabel('Exam1 score') ax.set_ylabel('Exam2 score' plt.show()
可以看到有明顯的分界線：
接下來就是正式的分析了 我們的目標是：建立分類器 ，求出3個引數（分數1的引數，分數2的引數，誤差） 設定閾值，根據閾值判斷是否錄取（這裡>=0.5就當作錄取）
要完成的模組： #1.sigmoid：對映到概率（0-1）的函式 g(z) = 1/(1+e^-z) #2.model: 返回預測結果值 #3.cost: 根據引數計算損失 #4.gradient: 計算每個引數的梯度方向 #5.descent: 進行引數更新 #6.accuracy: 計算精度
1.sigmoid 函式
值域（0-1） g(-∞) = 0 g(+∞) = 1 g(0)=0.5
def sigmoid(z): return 1/(1+np.exp(-z))
2.model
就是用資料作為sigmoid的輸入，從而構造出預測函 資料 (1,x1,x2) * (θ0 , θ1 , θ2)的轉置 = θ0+θ1X1+θ2X2 然後傳到 sigmoid中得到結果
def model(X,theta):
return sigmoid( np.dot(X,theta.T) )
#結果中 θ0需要和1相乘，資料中缺少1這1列，所以要補進
pdData.insert(0,'Ones',1); #在第0列（第一列之前，索引之後 插入列名為‘ones’，所有資料為1的一列
data = pdData.as_matrix() #將pdData從DataFrame 轉為矩陣，方便進行矩陣乘法
cols = data.shape[1] #資料有幾列 4 ones exam1 exam2 admitte
X = data[:,0:cols-1] #資料矩陣X 1，分數1，分數2
y = data[:,cols-1:cols]; # 結果（是否選上
theta=np.zeros([1,3]) # 3個引數 引數個數=資料個數+1
3.損失函式
這裡就是對數似然函式 * -1/m [加-是因為要把梯度上什轉為梯度下降，1/m是為了求平均] 似然函式 = 每個 x 對應的 P(y|x;θ) 的連乘 ， 然後取對數得到對數似然函式
損失函式（loss function）是用來衡量模型的預測值f(x)與真實值Y的不一致程度，它是一個非負實值函式，損失函式越小，模型越優（還需考慮過擬合等問題）
這裡選擇用 似然函式*-1/m 作為損失函式，是不是因為 似然函式越大越好 ， 加一個 負號就越小越好 ， 所以只要平均似然函式的值越來越小了，就說明越來越好了?
損失函式用： J(θ) = -1/m * sum( (y*log(hθ(X)) + (1-y)log(1-hθ(x)) ) )
def cost(X,y,theta):
tt = y-model(X,theta)
left = np.multiply(tt,tt)
return np.sum(left)/(len(X))
用 1/m * sum( (y - hθ(X) )^2 ) 做損失函式 效果相同
def cost(X,y,theta):
left = np.multiply(y,np.log(model(X,theta)
right = np.multiply(1-y,np.log(1-model(X,theta)))
return np.sum(left+right)/(-len(X)
4.求梯度
多維影象的梯度，也就是多維影象在某點各個方向的斜率中斜率最大的那個
演算法就是對J(θ) 求偏導 dJ(θ)/dθ = -1/m * sum( (yi - hθ(xi))xij )
def gradient(X,y,theta):
grad = np.zeros(theta.shape) #有幾個θ，就找幾個梯度
# dJ(θ)/dθ = -1/m * sum( (yi - hθ(xi))xij )
error = (model(X,theta) - y).ravel(); # yi - hθ(xi) 將-號提取進來 ravel將多維降為一維矩陣 ,因為原error是列向量，所以作用和專T類似
#求每個θi的梯度
for j in range (theta.shape[1]):
term = np.multiply(error,X[:,j]) #Xij 表示的是 X的每一個樣本，針對每一個引數θj 都只取和θj有關的那一列即X（ ， j）行不限
grad[0,j] = np.sum(term) / len(X)
return grad
4 （附屬） 停止策略
梯度更新需要有個頭，有個結束點
#設計了3種不同的停止策略
STOP_ITER = 0 #按照迭代次數停止 (一次迭代即更新一次梯度值)
STOP_COST = 1 #目標函式兩次迭代之間的損失（差異）很小
STOP_GRAD = 2 #兩次迭代之間梯度沒什麼變化
def stopCriterion(type,value,threshold):
#設定3種不同的停止策略
if type == STOP_ITER: return value > threshold
elif type == STOP_COST: return abs(value[-1]-value[-2]) < threshold
elif type == STOP_GRAD: return np.linalg.norm(value) < threshold; #求2範數 sqrt(x1^2+x2^2+.....)
4 （附屬） 對資料進行洗牌
打亂（防止自己收集資料的時候按照某種規律收集）
def shuffleData(data):
np.random.shuffle(data) #洗牌資料
#重新獲取資料
cols = data.shape[1] #資料有幾列 4 ones exam1 exam2 admitted
X = data[:,0:cols-1] #資料矩陣X 1，分數1，分數2
y = data[:,cols-1:cols]; # 結果（是否選上）
return X,y
5，6重點！
檢視不同梯度下降策略對時間的影響
import time
def descent(data, theta, batchSize, stopType ,thresh, alpha): #batchSize = 1 隨機梯度下降 =樣本個數 梯度下降 =1-樣本個數 小批量樣本夏婧
init_time = time.time();#開始時候的時間 #thresh 閾值 alpha學習率
i = 0;#迭代次數
k = 0;#batch ，即當前已經消耗的樣本數量
X,y = shuffleData(data)
grad = np.zeros(theta.shape) #梯度
costs = [cost(X,y,theta)] #損失值
#***5.descent: 進行引數更新***
while True:
grad = gradient(X[k:k+batchSize],y[k:k+batchSize],theta)
#更新引數 ， 讓θ 不停的靠近梯度為0的點
theta = theta - alpha*grad;
costs.append(cost(X,y,theta)) #用新的 θ 計算出新的損失值
i += 1; #迭代次數 +1
k+=batchSize
if k >= len(X):
k=0;
X,y = shuffleData(data); #重新洗牌
value = 0;
if stopType == STOP_ITER: value = i;
elif stopType == STOP_COST: value = costs;
elif stopType == STOP_GRAD: value = grad;
if stopCriterion(stopType,value,thresh):break; #是否跳出
return theta,i-1,costs,grad,time.time()-init_time
功能函式，畫圖
#功能函式，畫圖
def runExpe(data, theta, batchSize, stopType, thresh, alpha):
#import pdb; pdb.set_trace();
theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, alpha)
name = "Original" if (data[:,1]>2).sum() > 1 else "Scaled"
name += " data - learning rate: {} - ".format(alpha)
if batchSize==n: strDescType = "Gradient"
elif batchSize==1: strDescType = "Stochastic"
else: strDescType = "Mini-batch ({})".format(batchSize)
name += strDescType + " descent - Stop: "
if stopType == STOP_ITER: strStop = "{} iterations".format(thresh)
elif stopType == STOP_COST: strStop = "costs change < {}".format(thresh)
else: strStop = "gradient norm < {}".format(thresh)
name += strStop
print ("***{}\nTheta: {} - Iter: {} - Last cost: {:03.2f} - Duration: {:03.2f}s".format(
name, theta, iter, costs[-1], dur))
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(len(costs)), costs, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title(name.upper() + ' - Error vs. Iteration')
plt.show()
return theta
資料標準化
（（a-均值）/方差）
from sklearn import preprocessing as pp
scaled_data = data.copy()
scaled_data[:, 1:3] = pp.scale(data[:, 1:3])
n=100
#損失值 < 0.000001時停
theta = runExpe(scaled_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001)
6.accuracy: 計算精度
def predict(X,theda): #設定閾值 預測 大於0.5就是可以
return [1 if x>=0.5 else 0 for x in model(X,theta)]
scaled_X = scaled_data[:, :3]
y = scaled_data[:, 3]
predictions = predict(scaled_X, theta)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
accuracy = (sum(map(int, correct)) % len(correct))
print ('accuracy = {0}%'.format(accuracy))
