照斯坦福的機器學習課程用Octave敲了一下邏輯迴歸的演算法  然後搬到python上各種出問題...感覺自己對向量化 還是理解不夠徹底 又重新推導了一遍

用的資料集是課程的作業  training data包含100條資料 兩個特徵（兩門課程的分數） 一個類別（能否入學）

#載入資料集
data = np.loadtxt(path, dtype = float, encoding = 'utf-8', delimiter = ',')
X = data[:, 0:2]#二維陣列 100x2
y = data[:, 2]#一維陣列 1x100

需要把y轉換為100x1的列向量 方法有很多 

#y = np.c_[data[:, 2]]#100x1
y = data[:, 2].reshape(-1, 1)#100x1
 

邏輯迴歸假設函式：

#S形函式(h(x;theta))
def sigmoid(z):
	return 1 / (1 + np.exp(-z))
#100x1

X = [ ]100x3(新增x0 ==1), theta = [ ]1x3, X * theta'= [ ]100x1

代價函式：

向量化：

#代價函式
def costFunction(theta, X, y):
	m = y.size#100
	h = sigmoid(X.dot(theta.reshape(-1, 1)))#100x1

	J = -1 * (1 / m) * (np.log(h).T.dot(y) + np.log(1 - h).T.dot(1 - y))#1x1

	if np.isnan(J[0]):
		return(np.inf)
	return J[0]
 

相當於兩矩陣對應行點乘之後 求和

梯度 也就是代價函式對theta(j)的偏導:

向量化：

#計算梯度
def compute_grad(theta, X, y):
	m = y.size
	h = sigmoid(X.dot(theta.reshape(-1, 1)))#100x1

	grad = (1 / m) * ((h - y).T.dot(X))#應該是1x3的行向量 有三個theta值
	#但是打印出的結果為100x3的矩陣
	#===因為剛開始定義的y 是1x100的行向量
	return grad#1x3

梯度演算法 目的是求costFunction的最小值和最小值點的theta

用scipy.optimize中的minimize函式

from scipy.optimize import minimize#無約束最小化損失函式
res = minimize(costFunction, initial_theta, args = (X, y),
	jac = compute_grad, options = {'maxiter': 400})

#test1LogisticRegression.py:43: RuntimeWarning: 
#divide by zero encountered in log
#  J = -1 * (1 / m) * (np.log(h).T.dot(y) + np.log(1 - h).T.dot(1 - y))#1x1
#test1LogisticRegression.py:43: RuntimeWarning: divide by zero encountered in log
#  J = -1 * (1 / m) * (np.log(h).T.dot(y) + np.log(1 - h).T.dot(1 - y))#1x1
#      fun: 0.20349770158950983
# hess_inv: array([[ 2.85339493e+03, -2.32908823e+01, -2.27416470e+01],
#       [-2.32908823e+01,  2.04489131e-01,  1.72969525e-01],
#       [-2.27416470e+01,  1.72969525e-01,  1.96170322e-01]])
#      jac: array([-2.68557634e-09,  4.36433479e-07, -1.39671758e-06])
#  message: 'Optimization terminated successfully.'
#     nfev: 34
#      nit: 25
#     njev: 30
#   status: 0
#  success: True
#        x: array([-25.16131634,   0.2062316 ,   0.20147143])

[-25.16131634,   0.2062316 ,   0.20147143]

預測函式

def predict(theta, X, threshold = 0.5):
	p = sigmoid(X.dot(theta.T)) >= threshold
	return(p.astype('int'))

對分數為[30, 43]的同學進行預測

pre = sigmoid(np.array([1, 30, 43]).dot(res.x.T))

3.326083521580104e-05
#近似於零

參考：