先调用需要使用到的包
import matplotlib.pyplot as plt import numpy as np import pandas as pd import scipy.io as sio import matplotlib import scipy.optimize as opt from sklearn.metrics import classification_report
输入数据的函数，transpose为是否需要转置，由于所给数据是需要转置的，所以这里需要对数据进行一个转置
def load_data(path,transpose=True):
data = sio.loadmat(path)
y = data.get('y')
y = y.reshape(y.shape[0])
X = data.get('X')
if transpose: # 原始数据需要进行转置才能得到正确的方向
X = np.array([im.reshape((20,20)).T for im in X])
X = np.array([im.reshape(400) for im in X])
return [X,y]
载入数据，数据是 20 ∗ 20 20*20 20∗20 的图像，一共有 5000 5000 5000 个sample
[X,y] = load_data('ex3data1.mat')
print(X.shape)
print(y.shape)
(5000, 400) (5000,)
def PlotAnImage(X,y):
pick_one = np.random.randint(0,5000)
image = X[pick_one,:]
plt.show()
[fig,ax] = plt.subplots(figsize=(3,3))
ax.matshow(image.reshape((20,20)),cmap=matplotlib.cm.binary)
plt.xticks(np.array([]))
plt.yticks(np.array([]))
plt.show()
print('this is {}'.format(y[pick_one]))
随机挑一个图像输出看一下，注意， 0 0 0 的 lable 是 10 10 10
PlotAnImage(X,y)
this is 9
def Plot100Image(X):
sample_index = np.random.choice(np.arange(5000),100)
sample_image = X[sample_index,:]
[fig,ax]=plt.subplots(nrows=10,ncols=10,sharex=True,sharey=True,figsize=(8,8))
for row in range(10):
for col in range(10):
ax[row,col].matshow(sample_image[10*row+col].reshape((20,20)),cmap=matplotlib.cm.binary)
plt.xticks(np.array([]))
plt.yticks(np.array([]))
plt.show()
我们也可以输出更多的图形进行一个观察数据的合理性
Plot100Image(X)准备数据
先插入一列全为
0
0
0对应
x
0
x_0
x0
insert(被插数组, 指定位置, 插入数据, 插入维度)
print(X.shape[0]) X = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1) print(X.shape) print(X)
5000 (5000, 401) [[1. 0. 0. ... 0. 0. 0.] [1. 0. 0. ... 0. 0. 0.] [1. 0. 0. ... 0. 0. 0.] ... [1. 0. 0. ... 0. 0. 0.] [1. 0. 0. ... 0. 0. 0.] [1. 0. 0. ... 0. 0. 0.]]
print(y.shape) print(y)
(5000,) [10 10 10 ... 9 9 9]
给y向量化标签，具体看图“向量化标签”
y_matrix = [y_matrix[-1]] + y_matrix[:-1] 在这里y_matrix[-1]加上一个[]是因为y_matrix[-1]为一维数组，y_matrix[:-1]是2维数组，相加不能得到我们想要得到的内容
a = [[10,20],[30,40]] b = [1,2] print(a+b) # [[10, 20], [30, 40], 1, 2]
python 的 list 的 + 就是拼接两个list
y_matrix = []
for k in range(1,11):
y_matrix.append((y == k).astype(int))
y_matrix = [y_matrix[-1]] + y_matrix[:-1]
y = np.array(y_matrix)
print(y.shape)
print(y)
(10, 5000) [[1 1 1 ... 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 ... 1 1 1]]训练一维模型
Logistic Regression :
Sigmoid Function (
h
Θ
(
x
)
h_{Theta}(x)
hΘ(x))
1
1
+
e
−
z
frac{1}{1+e^{-z}}
1+e−z1
Cost Function
1
m
∑
i
=
1
m
[
−
y
(
i
)
l
o
g
(
h
Θ
(
x
(
i
)
)
)
−
(
1
−
y
(
i
)
)
l
o
g
(
1
−
h
Θ
(
x
(
i
)
)
)
]
frac{1}{m}sum_{i=1}^{m}[-y^{(i)}log(h_{Theta}(x^{(i)})) - (1-y^{(i)})log(1-h_{Theta}(x^{(i)}))]
m1i=1∑m[−y(i)log(hΘ(x(i)))−(1−y(i))log(1−hΘ(x(i)))]
Regularized Term
λ
2
m
∑
j
=
1
n
θ
j
2
frac{lambda}{2m}sum_{j=1}^{n}theta ^ 2 _j
2mλj=1∑nθj2
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def cost(theta, X, y):
return np.mean(-y * np.log(sigmoid(X @ theta)) - (1 - y) * np.log(1 - sigmoid(X @ theta)))
def regularized_cost(theta, X, y, Lambda=1):
penalized_theta = theta[1:]
regularized_term = (Lambda/(2*len(X))) * np.power(penalized_theta,2).sum()
return cost(theta,X,y) + regularized_term
∂ J ( θ ) ∂ θ j = 1 m ∑ i = 1 m [ h Θ ( x ( i ) ) − y ( i ) ] x j ( i ) frac{partial J(theta)}{partial theta _j}=frac{1}{m}sum _{i=1}^{m}[h_{Theta}(x^{(i)}) - y^{(i)}]x^{(i)}_j ∂θj∂J(θ)=m1i=1∑m[hΘ(x(i))−y(i)]xj(i)
在一维模型训练中： θ = ( 401 , 1 ) theta=(401,1) θ=(401,1) , X = ( 5000 , 401 ) X = (5000,401) X=(5000,401) , Y = ( 5000 , 1 ) Y = (5000,1) Y=(5000,1)
x j ( i ) x^{(i)}_j xj(i) 表示每一个sample的第 j j j个特征， i i i range from 1 1 1 to 5000 5000 5000
向量化之后，先计算出列向量 h Θ ( x ( i ) ) − y ( i ) h_{Theta}(x^{(i)}) - y^{(i)} hΘ(x(i))−y(i) 然后左乘 X T X^T XT, X T X^T XT一行就是一个特征的 θ theta θ
def gradient(theta, X, y):
return (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)
正则化项求导
d
(
λ
2
m
∑
j
=
1
m
θ
j
2
)
d
θ
=
λ
m
∑
j
=
1
m
θ
j
frac{mathrm{d} (frac{lambda}{2m}sum_{j=1}^mtheta _j ^2)}{mathrm{d} theta} = frac{lambda}{m}sum_{j=1}^mtheta _j
d θd (2mλ∑j=1mθj2)=mλj=1∑mθj
np.concatenate([a,b,c],axis=..) 在维度axis连接a,b,c
不对 θ 0 theta_0 θ0进行惩罚，所以求导的时候不带上他，求导值就为 0 0 0
在更新值的时候就不会收到正则化项的影响
θ j = θ j − α ∂ J ( θ ) ∂ θ j theta _j=theta _j - alpha frac{partial J(theta)}{partial theta_j} θj=θj−α∂θj∂J(θ) 其中 ∂ J ( θ ) ∂ θ j frac{partial J(theta)}{partial theta_j} ∂θj∂J(θ) 不带有 θ 0 theta_0 θ0的正则化项
def regularized_gradient(theta, X, y, Lambda):
penalized_theta = theta[1:]
regularized_theta = (Lambda / len(X)) * penalized_theta
regularized_term = np.concatenate([np.array([0]),regularized_theta])
return gradient(theta, X, y) + regularized_term
scipy.optimize.minimize
fun : 最小化的损失函数的名字
x0 : θ theta θ 的初始值，从传入这个 θ theta θ 进行 minimize
args : 优化的附加参数，传递给 fun, jac, hess 函数
method : 优化参数的方式，详情可以点进去看
options : disp 可以设置打印收敛消息，最后返回的值是一个结构体，可以调用结构体内的message查看，maxiter 可以设置最大的迭代次数
jac : 计算梯度的方法，也就是 θ theta θ 减去多少，也可以不填写，自动计算梯度，但是自动计算梯度速度很慢，并且正确率稍低，自动计算梯度花费56.3s，使用函数计算梯度花费0.5s，正确率都在99%以上
def logistic_regression(X, y, Lambda = 1):
theta = np.zeros(X.shape[1])
res = opt.minimize(
fun = regularized_cost,
x0 = theta,
args = (X, y, Lambda),
method = 'TNC',
jac = regularized_gradient,
options = {'disp':True}
)
final_theta = res.x
return final_theta
def predict(x, theta):
probability = sigmoid(x @ theta)
return (probability >= 0.5).astype(int)
theta0 = logistic_regression(X,y[0])
print(theta0.shape)
y_predict = predict(X,theta0)
print('Accuracy = {}'.format(np.mean(y[0] == y_predict)))
(401,) Accuracy = 0.9974训练k维模型
theta = np.array([logistic_regression(X,y[k]) for k in range(10)]) print(theta.shape) print(X.shape) print(y.shape)
(10, 401) (5000, 401) (10, 5000)
进行预测 N e w y = X θ T Newy = Xtheta^T Newy=XθT
probability_matrix = sigmoid(X @ theta.T) np.set_printoptions(suppress=True) print(probability_matrix) print(probability_matrix.shape)
[[0.99577193 0. 0.00053462 ... 0.00006468 0.00003912 0.00172204] [0.99834626 0.0000001 0.0000561 ... 0.0000968 0.0000029 0.00008493] [0.99139522 0. 0.00056788 ... 0.00000655 0.02652965 0.00197572] ... [0.00000068 0.04144051 0.00321396 ... 0.00012721 0.00297496 0.70764763] [0.00001844 0.00000013 0.00000009 ... 0.00164783 0.06820371 0.86121653] [0.02880042 0. 0.00012975 ... 0.36613909 0.0049755 0.14829976]] (5000, 10)
argmax是找到沿轴最大的值的索引，在这里就是找到哪一个概率最大，这里的索引正好对应哪一个数字
y_predict = np.argmax(probability_matrix,axis=1) print(y_predict) print(y_predict.shape)
[0 0 0 ... 9 9 7] (5000,)
[raw_x,raw_y] = load_data('ex3data1.mat')
y_answer = raw_y
y_answer[y_answer == 10] = 0
print(y_answer)
[0 0 0 ... 9 9 9]
print(classification_report(y_answer,y_predict))
precision recall f1-score support
0 0.97 0.99 0.98 500
1 0.95 0.99 0.97 500
2 0.95 0.92 0.93 500
3 0.95 0.91 0.93 500
4 0.95 0.95 0.95 500
5 0.92 0.92 0.92 500
6 0.97 0.98 0.97 500
7 0.95 0.95 0.95 500
8 0.93 0.92 0.92 500
9 0.92 0.92 0.92 500
accuracy 0.94 5000
macro avg 0.94 0.94 0.94 5000
weighted avg 0.94 0.94 0.94 5000
