机器学习（一）逻辑回归与softmax回归及代码示例

本文适合已经对机器学习、人工智能有过一定了解，但是还没有自己写过代码，或者一直在使用现有框架的同学。不用框架自己写一次代码的过程还是很有必要的，能让你真正地理解原理与机器学习中各个步骤的实现过程，而不是停留在“好像懂了”、只会调库的阶段。

目录

 

一、logistics回归简介（仅理论）

二、softmax回归（softmax_model.py文件）

2.1 代价函数

2.2 批量梯度下降（反向传播）

2.3 正则化（权重衰减）

2.4 mnist数据集的提取

2.5模型训练与保存

2.6 准确度评估

三、模型运用（Application.py文件)

3.1模型的加载

3.2 鼠标写数字并保存为图像

3.3 模型的使用

四、模型缺陷

---

一、logistics回归简介（仅理论）

线性回归输出的值可以很大，很难确定一个阈值作为分类线。因此使用logistics回归，将模型的输出变量映射到0和1之间，就能以 0.5 作为分类的阈值。

逻辑回归模型的假设是： 其中： ​X 代表特征向量 , g 代表逻辑函数. 一个常用的逻辑函数为S形函数（Sigmoid function），公式为：

simoid函数的图像为：

其中 (也可以用更复杂的组合)

代码实现

import numpy as np
  ​
  def sigmoid(z):
      return 1/(1+np.exp(-z))

构造一个代价函数 J(θ)，衡量模型的输出值与真实值的差别大小。

有以下两种情况

• 当类标签 y 为1时， h(x)越接近1， J(θ)越小； h(x)越接近0， J(θ)越大

• 当类标签 y 为1时， h(x)越接近1， J(θ)越小； h(x)越接近0， J(θ)越大

这其实就是个数学建模中的 0-1 型整数规划问题

可根据这种特性构造代价函数：

 

二、softmax回归（softmax_model.py文件）

softmax回归模型是logistics回归模型在多分类问题上的推广 ，类标签 可以取两个以上的值

2.1 代价函数

由于 y 有k个类别，故 变成一个 维的向量，来表示这 个类别估计的概率值, 并归一化（概率和为1）

def softmax(z):
      [m,k]= z.shape
      p = np.zeros([m,k])
      for i in range(m):
          p[i,:]=np.exp(z[i,:])/np.sum(np.exp(z[i,:]))
      return p

将logistics回归的 0-1 整数规划思想（二维），拓展为 k维 , 引入一个指示函数 1{·}

其取值规则为： 1{表达式为真} = 1

利用指示函数，构造出了代价函数

 

def cost(theta, x,y):  #x(m行n列），y（m行k列），theta（k行n列）
      [k,m]=y.shape
      theta = np.matrix(theta)
      sum = 0
      p = softmax(np.dot(x,theta.T))  ## [m,k]     P(i)为 [1,k]
      p = p.T.reshape([k*m,1])
      y = y.reshape([k*m,1])
      temp_p=np.mat(np.log(p))
      #punish = np.sum(np.power(theta,2))   #惩罚项，在梯度下降时求导，这里不必加上去
      cost = -1/m*np.dot(y.T,temp_p) #+ punish
      return cost         #输出m行k列，代表m个样本，k个类别各自概率

2.2 批量梯度下降（反向传播）

利用批量梯度下降（BGD）的思想更新 θ：

其中，代价函数的梯度 经过公式推导为：

​

  def gradientDescent(x,y,theta,iters,alpha,regulization_rate):
      #theta:权重系数    iters:迭代次数    alpha:学习率
      COST = np.zeros((iters,1))    #存放每次迭代后，cost值的变化
      thetaNums = int(theta.shape[0])  #维数，即j的取值个数
      print(thetaNums)
      for i in range(iters):
          #bb = x*theta.T
          p = softmax(np.dot(x,theta.T));
          grad = (-1/m*  np.dot(x.T ,(y.T-p)) ).T   #[3,784]
          #更新theta
          theta = theta - alpha*grad #- regulization_rate * theta
          COST[i] = cost(theta,x,y)
          #每训练一次，输出当前训练步数与损失值
          print("训练次数： ",i)
          print(COST[i])
          print("\n")
      #返回迭代后的theta值，和每次迭代的代价函数值
      return theta,COST

 

2.3 正则化（权重衰减）

在代价函数中添加一个权重衰减项 ，使代价函数变成凸函数，以期能够收敛到最优解

现在代价函数变为：

对应代码更新：

  punish = np.sum(np.power(theta,2))   
  cost = -1/m*np.dot(y.T,temp_p) + punish

 

对应地，梯度公式变为：

 

也就是说，更新 θ 时，需要再减去一个 λθ ，公式变为：

对应代码更新：

  theta = theta - alpha*grad - regulization_rate * theta

 

2.4 mnist数据集的提取

读取图像数据方法：

  def loadImageSet(filename):
      print("load image set", filename)
      binfile = open(filename, 'rb')
      buffers = binfile.read()
  ​
      head = struct.unpack_from('>IIII', buffers, 0)
      print("head,", head)
  ​
      offset = struct.calcsize('>IIII')
      imgNum = head[1]
      width = head[2]
      height = head[3]
      # [60000]*28*28
      bits = imgNum * width * height
      bitsString = '>' + str(bits) + 'B'  # like '>47040000B'
  ​
      imgs = struct.unpack_from(bitsString, buffers, offset)
  ​
      binfile.close()
      imgs = np.reshape(imgs, [imgNum, 1, width * height])
      print("load imgs finished")
      return imgs

读取标签方法：

  def loadLabelSet(filename):
      print("load label set", filename)
      binfile = open(filename, 'rb')
      buffers = binfile.read()
  ​
      head = struct.unpack_from('>II', buffers, 0)
      print("head,", head)
      imgNum = head[1]
  ​
      offset = struct.calcsize('>II')
      numString = '>' + str(imgNum) + "B"
      labels = struct.unpack_from(numString, buffers, offset)
      binfile.close()
      labels = np.reshape(labels, [imgNum, 1])
  ​
      print('load label finished')
      return labels

加载mnist数据集：

  def data_Load():
      imgs_train = loadImageSet("train-images.idx3-ubyte")
      imgs_train = imgs_train.reshape(-1, 784)
      labels_train = loadLabelSet("train-labels.idx1-ubyte")
      imgs_test = loadImageSet("t10k-images.idx3-ubyte")
      imgs_test = imgs_test.reshape(-1, 784)
      labels_test = loadLabelSet("t10k-labels.idx1-ubyte")
      return imgs_train,labels_train,imgs_test,labels_test

并在主函数中调用

[imgs_train, labels_train, imgs_test, labels_test] = data_Load()  # 读取训练集与测试集
x = imgs_train[0:num_train, :]   #num_train为读取的样本的数量
#将灰度图转化为二值图
x[x<=40]=0;
x[x>40]=1;
y = np.zeros((num_train,10))    
#mnist每个样本的标签是一个数字，以独热码的形式将其扩充为10维的向量
for t in np.arange(0,num_train):
    y[t,labels_train[t]]=1;
y=y.T

 

2.5模型训练与保存

训练模型

首先进行参数的设置（没来得及写参数调优可视化的代码了，多次尝试后设定成如下数值）：

iters = 1000;       #训练次数
alpha = 0.5;        #学习率
regulization_rate = 0.1;    #惩罚系数
k = 10;             #标签类别数
num_train = 2000    #样本数量

调用2.2写好的​gradientDescent() 方法完成训练，并得到训练后的 θ 与损失值 cost

[finalTheta,finalCost] = gradientDescent(x,y,theta,iters,alpha,regulization_rate)

保存模型

import pickle
  ​
  f = open('saved_model/model.pickle', 'wb') #保存到文件夹saved_model中
      pickle.dump(finalTheta, f)
      f.close()

 

2.6 准确度评估

y_predict = np.argmax(np.dot(imgs_test, finalTheta.T), axis=1)
      for i in range(len(labels_test)):
          if (y_predict[i]==labels_test[i]):
              num = num+1
      print("准确率为 ：",num/len(labels_test))

运行结果：

准确率为 0.8751 （对于mnist自带的测试集）

训练耗时 151秒

三、模型运用（Application.py文件)

3.1模型的加载

def model_load(filepath):
      f = open(filepath,'rb')
      w = pickle.load(f)
      f.close()
      return w
  ##主函数中
  theta = model_load('saved_model/model.pickle')

filepath 为之前保存模型的路径'saved_model/model.pickle'

3.2 鼠标写数字并保存为图像

鼠标绘图方法

drawing = False  # true if mouse is pressed
  ix, iy = -1, -1
  # 鼠标手写数字
  def draw_circle(event, x, y, s,a):
      global ix, iy, drawing, mode
      if event == cv2.EVENT_LBUTTONDOWN:
          drawing = True
          ix, iy = x, y
      elif event == cv2.EVENT_MOUSEMOVE:
          if drawing == True:
              #圆半径（笔画粗细）设置为25是为了和mnist数据集中的数字尽可能粗细相似
              cv2.circle(image, (x, y), 25, (255, 255, 255), -1)    #画笔颜色为白色
      elif event == cv2.EVENT_LBUTTONUP:
          drawing = False
          cv2.circle(image, (x, y), 25, (255, 255, 255), -1)
  image = np.zeros((512, 512, 3), np.uint8)       #初始化画布，黑色
  cv2.namedWindow('image')
  cv2.setMouseCallback('image', draw_circle)

 

3.3 模型的使用

思路

  1.展现一张512x512大小的黑色画布

  2.利用opencv的鼠标回调方法，鼠标在画布上按下时绘制一个小圆，随着鼠标的拖动，完成数字的书写

  3.当书写完成，按下键盘'p'键，程序将当前画布上的图像保存到本地

  4.将保存的图像进行预处理，转化成mnist数据集的格式

  5.把处理好的图像矩阵交给模型进行识别

图片预处理

def pre_picture(picName):
      img = cv2.resize(picName,(28,28))    #将图像等比例缩小至28x28大小
      im_arr = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)   #转化为灰度图
      x_input = im_arr.reshape([1,784])               #把图像矩阵转化为一维向量
      x_input[x_input < 40] = 0                       #设置灰度阈值，将向量二值化
      x_input[x_input>=40]=1                          ##
      return x_input                                  #返回图像处理后的生成的特征向量

---

鼠标书写数字并识别

while(1):
          cv2.imshow('image', image)
          k = cv2.waitKey(1) & 0xFF
          #按下 p 键，识别当前画布上的数字
          if k == ord('p'):
              #img_new = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
              cv2.imwrite("REC.jpg", image)
              x = pre_picture(image)
              labels_pro=np.dot(x,theta.T)
              y = np.argmax(np.dot(x,theta.T),axis=1)
              #输出模型识别的结果
              print("手写数字识别结果 ： ",y)
              print("各标签概率  = ",labels_pro)
          #按下 c 键，清空当前图像，还原成黑色画布
          elif k == ord('c'):
              image = np.zeros((512, 512, 3), np.uint8)#清空
          #按下 “ESC” 键，退出程序
          elif k == 27:
              break

运行截图

 

四、模型缺陷

一个很大的问题是mnist中数字的字体风格是西方式的，和中国人常用手写体有很大的区别。因此当使用Application.py文件进行手写数字识别时，必须要刻意去模仿mnist的字体风格，才能提高识别准确率。

 

嘛，直接用全像素点训练的后果就是这样啦。有两个解决方向：（1）选用表现力好的特征 （2）使用更复杂的模型

如果用CNN等神经网络做的话显然效果会好很多。以后有时间我会更新。