手撸机器学习系列之 - 多元逻辑回归不加正则项

这个例子是逻辑回归不加正则项的编写笔记

首先导入依赖

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>

png

打出一张图片是像上面这样

下面随机打出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)))

png

上面这些数字显示好像有点不对,可能数据需要进一步处理,这里我们将数据转置即可. 没有为什么,人家数据就是转置之后存的,我们用的时候要转回来.

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)

png

使用 .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()

png

添加参数 ,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>]

png

从图中来看,是正确的呢

接下里,我们来定义 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))

png

上图中,图片左上角的即是图片的预测值

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值