这个例子是逻辑回归不加正则项的编写笔记
首先导入依赖
import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
加载数据, 加载的数据是一个字典
data = sio.loadmat('ex3data1.mat')
data,type(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)},
dict)
raw_X=data['X']
raw_X 是 5000 张图片,400维度,像 [0., 0., 0., …, 0., 0., 0.] 这样
raw_X.shape,raw_X
((5000, 400), 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.]]))
random_index=np.random.randint(5000)
fig,ax=plt.subplots(figsize=(1,1))
ax.imshow(raw_X[random_index].reshape((20,20)))
<matplotlib.image.AxesImage at 0x272e1bddcc8>
打出一张图片是像上面这样
下面随机打出100 张图片
type(raw_X),raw_X.shape
(numpy.ndarray, (5000, 400))
img_indexes=np.random.choice(len(raw_X),100)
imges=raw_X[img_indexes,:]
fig,ax=plt.subplots(figsize=(8,8),ncols=10,nrows=10)
for i in range(10):
for j in range(10):
img_index=i*10+j
ax[i,j].imshow(imges[img_index].reshape((20,20)))
上面这些数字显示好像有点不对,可能数据需要进一步处理,这里我们将数据转置即可. 没有为什么,人家数据就是转置之后存的,我们用的时候要转回来.
img_indexes=np.random.choice(len(raw_X),100)
imges=raw_X[img_indexes,:]
fig,ax=plt.subplots(figsize=(8,8),ncols=10,nrows=10)
for i in range(10):
for j in range(10):
img_index=i*10+j
ax[i,j].imshow(imges[img_index].reshape((20,20)).T)
使用 .T 即可获得其转置矩阵,然后再打出来,就好了. 另外,绘制矩阵图像用 imshow() 函数
这里有坐标轴比较难受,我们把坐标轴去掉
img_indexes=np.random.choice(len(raw_X),100)
imges=raw_X[img_indexes,:]
fig,ax=plt.subplots(figsize=(8,8),ncols=10,nrows=10,sharex=True,sharey=True)
for i in range(10):
for j in range(10):
img_index=i*10+j
ax[i,j].imshow(imges[img_index].reshape((20,20)).T)
plt.xticks([])
plt.yticks([])
plt.show()
添加参数 ,sharex=True,sharey=True
并 调用 plt.xticks([])
和 plt.yticks([])
就可以去除坐标轴了
接下来定义 sigmoid 函数
def sigmoid(z):
return 1/(1+np.exp(-z))
验证一下是否正确
fig,ax=plt.subplots(figsize=(10,10))
test_x=np.linspace(-10,10,100)
test_y=sigmoid(test_x)
ax.plot(test_x,test_y)
[<matplotlib.lines.Line2D at 0x272ea7fd388>]
从图中来看,是正确的呢
接下里,我们来定义 costFun
代价函数是: [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5yBgruXw-1600668076601)(attachment:%E5%9B%BE%E7%89%87.png)]
h θ ( x ) h\theta(x) hθ(x) 是: [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-I7qp7HwH-1600668076602)(attachment:%E5%9B%BE%E7%89%87.png)]
g(x) 是: [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mYHYzfne-1600668076603)(attachment:%E5%9B%BE%E7%89%87.png)]
def costFunc1(theta,X,y,alpha):
item1= -y*np.log(sigmoid(theta.T @ X))
item2= (1-y)*log(1-sigmoid(theta.T @ X))
cost=(1/lex(X))*np.sum(item1-item2)
上面这个函数就是比着公式理解得来的
接下来,我们在 X 的第一列插入 1, 也就是 x0
X=np.insert(raw_X,0,values=1,axis=1)
X.shape
(5000, 401)
接下来,咋们初始化theta
theta=np.zeros((len(X)))
theta.shape,theta
((5000,), array([0., 0., 0., ..., 0., 0., 0.]))
把 y 拿出来
y=data['y']
y.shape,y
((5000, 1), array([[10],
[10],
[10],
...,
[ 9],
[ 9],
[ 9]], dtype=uint8))
y现在是二维数组, 5000 行,1列,我们他把处理为一维(一个数组)
y=y.flatten()
y.shape,y
((5000,), array([10, 10, 10, ..., 9, 9, 9], dtype=uint8))
这就变成一维数组了
看下 theta
theta
array([0., 0., 0., ..., 0., 0., 0.])
计算下损失,这里用一个第三方库
from scipy.optimize import minimize
def one_vs_all(X,y,K):
# 参数个数,就本例来说,有401个参数
n = X.shape[1]
# 初始化 10 维矩阵, 10个模型, K 为模型个数
theta_all = np.zeros((K,n))
for i in range(1,K+1):
# 初始化一个一维矩阵
theta_i = np.zeros(n,)
# 优化
res = minimize(fun =costFunc,
x0 = theta_i,
args = (X, y == i),
method = 'TNC',
jac =gredient)
# 把优化后的值给 theta_all
theta_all[i-1,:] = res.x
return theta_all
# 定义损失函数
def costFunc(theta,X,y):
equal=X@theta
hx=sigmoid(equal)
item1=y*np.log(hx)
item2=(1-y)*np.log(1-hx)
return (1/(len(X)))*-np.sum(item1+item2)
# 定义梯度下降
def gredient(theta,X,y):
return (X.T@(sigmoid(X@theta) - y)) / len(X)
lamda = 1
K =10
X.shape,y.shape
theta_final=one_vs_all(X,y,K)
C:\Anaconda3\lib\site-packages\ipykernel_launcher.py:31: RuntimeWarning: divide by zero encountered in log
C:\Anaconda3\lib\site-packages\ipykernel_launcher.py:31: RuntimeWarning: invalid value encountered in multiply
C:\Anaconda3\lib\site-packages\ipykernel_launcher.py:30: RuntimeWarning: divide by zero encountered in log
C:\Anaconda3\lib\site-packages\ipykernel_launcher.py:30: RuntimeWarning: invalid value encountered in multiply
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-a4yY8l5p-1600668076604)(attachment:%E5%9B%BE%E7%89%87.png)]
预测一下
theta_final
array([[-6.05726801e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.84445122e-02, 2.81693595e-07, 0.00000000e+00],
[-5.84716907e+00, 0.00000000e+00, 0.00000000e+00, ...,
6.02797727e-02, -6.35474074e-03, 0.00000000e+00],
[-8.58433297e+00, 0.00000000e+00, 0.00000000e+00, ...,
-2.55725115e-04, -1.26557197e-06, 0.00000000e+00],
...,
[-1.32887557e+01, 0.00000000e+00, 0.00000000e+00, ...,
-6.69606820e+00, 7.73348511e-01, 0.00000000e+00],
[-8.68350387e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.91978559e-01, 9.67706194e-03, 0.00000000e+00],
[-1.25484541e+01, 0.00000000e+00, 0.00000000e+00, ...,
1.16414470e-03, 5.21182784e-05, 0.00000000e+00]])
X.shape,theta_final.shape
((5000, 401), (10, 401))
# 预测 X.shape,theta_final.shape
# ((5000, 401), (10, 401))
lineEqual=X@theta_final.T
hx=sigmoid(lineEqual)
# 从列的角度 (axis=1)找出 hx 中值最大的索引值. 索引值 +1 就得到了实际值. 因为在处理实际值时,
# 我们就是按这个套路来的
max_index = np.argmax(hx,axis=1)
pre_num=max_index+1
pre_num
array([10, 10, 10, ..., 9, 9, 9], dtype=int64)
y,pre_num
(array([10, 10, 10, ..., 9, 9, 9], dtype=uint8),
array([10, 10, 10, ..., 9, 9, 9], dtype=int64))
计算下准确率
correct=[1 if a==b else 0 for (a,b) in zip(y,pre_num)]
accuacy=np.sum(correct)/len(correct)
accuacy
0.974
我们将图片和预测可视化看看
fig,ax=plt.subplots(ncols=10,nrows=10,figsize=(10,10),sharex=True,sharey=True)
image_indexs=np.random.choice(len(raw_X),100)
images=raw_X[image_indexs,:]
# 配合 sharex=True,sharey=True 隐藏坐标
plt.xticks([])
plt.yticks([])
# 注意,0 = 10
for i in range(10):
for j in range(10):
f_index=i*10+j
ax[i,j].imshow(images[f_index].reshape((20,20)).T)
# pre_num[image_indexs[f_index]] 为图片的预测值
ax[i,j].text(2, 4, pre_num[image_indexs[f_index]] , horizontalalignment='center',bbox=dict(facecolor='red', alpha=0.5))
上图中,图片左上角的即是图片的预测值