import pandas as pd import numpy as np import matplotlib.pyplot as plt from scipy.io import loadmat
数据集以matlab格式储存，需要Scipy工具
加载数据data=loadmat('ex3data1.mat')
data
{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
'__version__': '1.0',
'__globals__': [],
'X': array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]),
'y': array([[10],
[10],
[10],
...,
[ 9],
[ 9],
[ 9]], dtype=uint8)}
data['X'].shape,data['y'].shape
((5000, 400), (5000, 1))数据可视化
x=data['X'] s = np.random.permutation(x) a = s[:100,:] #选前50行,(50, 400)
#定义可视化函数
def displayData(x):
plt.figure()
n = np.round(np.sqrt(x.shape[0])).astype(int)
fig, ax = plt.subplots(nrows=n, ncols=n, sharex=True, sharey=True, figsize=(6, 6))
#在每个子画布中画出一个数字
for row in range(n):
for column in range(n):
ax[row, column].imshow(x[10 * row + column].reshape(20,20).T, cmap='gray')
plt.xticks([]) #去掉坐标轴
plt.yticks([])
plt.show()
displayData(a)
矢量化逻辑回归
def sigmoid(z):
return 1.0/(1.0+np.exp(-z))
sigmoid(0)
0.5
def cost(theta, X, y, learningRate):
# STEP1：将theta, X, y转换为numpy类型的矩阵
theta =np.matrix(theta)
X = np.matrix(X)
y = np.matrix(y)
# STEP2：根据公式计算损失函数（不含正则化）
cross_cost =np.multiply(-y, np.log(sigmoid(X * theta.T)))-np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
# STEP3：根据公式计算损失函数中的正则化部分
reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[1:], 2))
# STEP4：把上两步当中的结果加起来得到整体损失函数
# your code here (appro ~ 1 lines)
whole_cost=np.sum(cross_cost)/len(X)+reg #sum不要忘掉 ，否则最后准确率差别很大
return whole_cost
def gradient(theta, X, y, learningRate):
# STEP1：将theta, X, y转换为numpy类型的矩阵
theta = np.matrix(theta)
X = np.matrix(X)
y = np.matrix(y)
# STEP2：将theta矩阵拉直（转换为一个向量）
parameters =int(theta.ravel().shape[1]) #theta数量
# STEP3：计算预测的误差
error = sigmoid(X * theta.T) - y
# STEP4：根据上面的公式计算梯度
grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
# STEP5：由于j=0时不需要正则化，所以这里重置一下
grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
return np.array(grad).ravel()
from scipy.optimize import minimize
def one_vs_all(X, y, num_labels, learning_rate):
rows = X.shape[0]
params = X.shape[1]
# k X (n + 1) array for the parameters of each of the k classifiers
all_theta = np.zeros((num_labels, params + 1))
# insert a column of ones at the beginning for the intercept term
X = np.insert(X, 0, values=np.ones(rows), axis=1) #插了一列1
# labels are 1-indexed instead of 0-indexed
for i in range(1, num_labels + 1):
theta = np.zeros(params + 1)
y_i = np.array([1 if label == i else 0 for label in y])
y_i = np.reshape(y_i, (rows, 1))
# minimize the objective function
fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient) # 参数位置保证正确
all_theta[i-1,:] = fmin.x
return all_theta
rows = data['X'].shape[0] params = data['X'].shape[1] all_theta = np.zeros((10, params + 1)) theta = np.zeros(params + 1) X_0=np.insert(data['X'],0,1.0,axis=1) y_0 = np.array([1 if label == 0 else 0 for label in data['y']]) y_0 = np.reshape(y_0, (rows, 1)) X_0.shape, y_0.shape, theta.shape, all_theta.shape
((5000, 401), (5000, 1), (401,), (10, 401))
#看y有几类标签 np.unique(data['y'])
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=uint8)
all_theta = one_vs_all(data['X'], data['y'], 10, 1) all_theta
array([[-2.37935917e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.29766577e-03, 8.08552795e-10, 0.00000000e+00],
[-3.17481434e+00, 0.00000000e+00, 0.00000000e+00, ...,
4.51583692e-03, -5.14550467e-04, 0.00000000e+00],
[-4.79789731e+00, 0.00000000e+00, 0.00000000e+00, ...,
-2.87494906e-05, -2.47839422e-07, 0.00000000e+00],
...,
[-7.98362233e+00, 0.00000000e+00, 0.00000000e+00, ...,
-8.85576929e-05, 7.12346918e-06, 0.00000000e+00],
[-4.58713014e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.39122413e-03, 1.02609812e-04, 0.00000000e+00],
[-5.31727020e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.14812908e-04, 8.49987514e-06, 0.00000000e+00]])
def predict_all(X, all_theta):
# STEP1：获取矩阵的维度信息
rows = X.shape[0]
params = X.shape[1]
num_labels = all_theta.shape[0]
# STEP2：把矩阵X加入一行零元素
X = np.insert(X, 0, values=np.ones(rows), axis=1)
# STEP3：把矩阵X和all_theta转换为numpy型矩阵
X = np.matrix(X)
all_theta = np.matrix(all_theta)
# STEP4：计算样本属于每一类的概率
h = sigmoid(X * all_theta.T)
# STEP5：找到每个样本中预测概率最大的值
h_argmax = np.argmax(h, axis=1)
# STEP6：因为我们的数组是零索引的，所以我们需要为真正的标签+1
h_argmax = h_argmax + 1
return h_argmax
y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))
accuracy = 94.46%
