基于Logistic回归和Sigmoid函数的分类

对于二分类问题，需要找到这样一个函数，他对任意给定的输入输出0或1，海维赛德阶跃函数符合这一特性，但其在跳跃点处很难处理，因此我们用其他函数来代替，本章所用的函数为sigmoid函数该函数的表达式及图像如下所示

可以看到，当x趋于正无穷时，函数值趋近于1；当x趋于负无穷时，函数值趋近于0；

我们将输出大与0.5的归入1类，将输出小于0.5时归入0类。

假设输入数据的特征是(x0, x1, x2, …, xn)，我们在每个特征上乘以一个回归系数 (w0, w1, w2, … , wn)

利用线性代数的知识，我们可以把式子简化为：

​

再将Z输入sigmoid函数就可以得到最终的值啦~

可以看到，式子中最重要的是要确定w(又称权重矩阵)，那么应该怎么确定这个权重矩阵呢？

梯度上升法

梯度上升法的思想是：要找到某函数的最大值，最好的办法就是沿着该函数的梯度方向探寻。

梯度上升算法的公式如下，这里的f(w)为损失函数，关于损失函数的知识，可以参照吴恩达的机器学习视频或者网上的梯度上升算法的推导。损失函数代表了真实值与预测值之间的差异，我们对其求偏导即可得到权重矩阵的梯度，这里的alpha是学习率，学习率的选择是一项很大的学问，本章暂时给定学习率，只需知道选择合适的学习率是一个很重要的步骤，学习率太小会导致开销大，学习率太大会导致越来越偏离真实值。

以下是梯度上升法的具体例子

OK，有了理论支撑我们就可以开始写代码了

简单练习

作为练习，我们首先使用一个灰常简单的数据集

这是test.txt中的一小部分数据，可以看到数据有两个特征，最后有一个0,1的标签

先看一下我们需要导入的包，因为是从底层开始写，所以我们需要的特别简单（环境为python2.7）

from numpy import *		##实际不推荐这种写法
from math import exp

整个文档共有100个这样的数据，我们将其导入

def loaddataset():
    datamat = []; labelmat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
            linearr = line.strip().split()
            datamat.append([1.0,float(linearr[0]),float(linearr[1])])
            labelmat.append(int(linearr[2]))
    return datamat,labelmat

定义sigmoid函数

def sigmoid0(inX):
    return 1.0/(1+exp(-inX))

然后写出梯度上升算法

"""未优化梯度上升算法"""
def gradascent(datamatin,classlabels):
    datamatrix = mat(datamatin)
    labelmat = mat(classlabels).transpose()
    m,n = shape(datamatrix)
    alpha = 0.001
    maxcycles = 500
    weights = ones((n,1))
    for k in range(maxcycles):
        h = sigmoid(datamatrix*weights)
        error = (labelmat - h)
        weights = weights + alpha*datamatrix.transpose()*error
    plotbestfit(weights)
    return weights

在运行之前，我们先对sigmoid函数做一个改写。可以看到传入sigmoid函数的参数dataamatrixweights是一个m1的矩阵，但是Python3.5并不支持用exp ()函数直接对矩阵进行操作。这一点书中并没有指出，算是一个勘误。以下给出风车自己的方法，这个方法有一点蠢，既然不支持直接对矩阵中的每个元素进行操作，那我们就先把它转换为numpy中的数组，然后用数组遍历的方法来进行每个元素的sigmoid操作，最后别忘记再把它转成矩阵。修改后的sigmoid()函数如下所示

"""求sigmoid"""
def sigmoid(data):
    s = data.getA()
    datatemp = [[ 1.0/(1+exp(-s[i][j])) for j in range(len(s[i]))]for i in range(len(s))]
    datafinal = mat(datatemp)
    return datafinal

好了，这就可以开始测试了

"""测试代码"""
dataarr,labelmat = loaddataset
print(gradascent(dataarr,labelmat))

输出如下：

可以看到，我们的方法确实是有效的，但是写了这么久看到一行黑白的输出，是不是有些小失落，别急，我们来让这些数据可视化

绘制图形

这里直接给出代码，用到的只是一些基础的功能，搜索一下相关函数的功能很容易理解。风车也还在学进一步的图形绘制功能，一起加油～

"""画出决策边界"""
def plotbestfit(wei):
    import matplotlib.pyplot as plt
    weights = wei
    datamat,labelmat = loaddataset()
    dataarr = array(datamat)
    n = shape(dataarr)[0]
    xcord1 = []; ycord1 = []
    xcord2 = []; ycord2 = []
    for i in range(n):
        if int(labelmat[i]) == 1:
            xcord1.append(dataarr[i,1]);ycord1.append(dataarr[i,2])
        else:
            xcord2.append(dataarr[i,1]);ycord2.append(dataarr[i,2])
    fig = plt.figure()
    ax =fig.add_subplot(111)
    ax.scatter(xcord1,ycord1,s=30,c='red',marker='s')
    ax.scatter(xcord2,ycord2,s=30,c='green')
    x = arange(-3.0,3.0,0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('x1');plt.ylabel('x2');
    plt.show()

我们在程序末尾加上

plotbestfit(gradascent(dataarr,labelmat))

运行得到如下图像

可以看到结果对数据的划分还是很令人满意的

随机梯度上升算法

梯度上升法给出的结果很棒，但其本身仍有有一定的不足。梯度上升算法在每次更新权重矩阵时都要便利整个数据集，我们给出的test.txt只有100个样例，但若是有数十亿个样本，那么这样的方式复杂度就太高了。因此我们有了随机梯度算法。随机梯度算法，顾名思义，每次随机的取一个值（为了方便，下面的样例按顺序取一个值）来对我们的权重矩阵进行更新。

随机梯度上升算法的代码如下

"""随机梯度上升算法"""
def stogradascent(datamatrix,classlabels,num=150):
    m,n = shape(datamatrix)
    weights = ones(n)
    for j in range(num):
        dataindex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001
            randindex = int(random.uniform(0,len(dataindex)))
            """此处代码为自己修改"""
            h = sigmoid0(sum(datamatrix[randindex]*weights))
            """sigmoid0()出现OverflowError:math range error"""
            error = classlabels[randindex] - h
            weights = weights + alpha * error * datamatrix[randindex]
            del(dataindex[randindex])
    return weights

这里的sigmoid0函数就是开头给出的

def sigmoid0(inX):
    return 1.0/(1+exp(-inX))

这一次我们传入的值sum(datamatrix[randindex]weights)并不是一个矩阵，sum()函数内两个矩阵的大小分别为1n和n1，因此得到的是一个11的矩阵，我们用这个值来更新我们的权重矩阵weights。

最后输入测试代码

dataarr,labelmat = loaddataset()
weights = stogradascent(array(dataarr),labelmat)
plotbestfit(weights)

可以看到如下结果

可以看到结果并没有梯度上升算法好，但这样的比较是不公平的，因为后者是在整个数据集上迭代了500次得到的，当然提升随机梯度算法的迭代次数num也能的到更好的结果。

小测试也做完了，这回来个实际应用吧！

从疝气病症预测病马的死亡率

给出的horseColicTraining.txt中每个样本有20个特征，并有一个标签，0代表未能存活，1代表可以存活。horseColicTest.txt作为测试集来判断我们训练出的模型的好坏。

话不多说，直接上代码

def classifyvector(inX,weights):
    prob = sigmoid0(sum(inX*weights))
    if prob > 0.5:return 1.0
    else:return 0.0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = [];trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stogradascent(array(trainingSet),trainingLabels,1000)
    errorCount = 0; numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr=[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyvector(array(lineArr),trainWeights)) != int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)
    print("the error rata of this test is:%f" %errorRate)
    return errorRate
def multiTest():
    numTests = 10; errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rata is:%f" %(numTests,errorSum/float(numTests)))

"""测试代码"""
multiTest()

运行结果如下：

tensorflow实现

数据集

本次我们所使用的数据集，是MNIST数据集。MNIST是一个手写数字数据库，它有60000个训练样本集和10000个测试样本集。它是NIST数据库的一个子集。其中每张图片固定大小为28281，最后的数字1表示一个通道，即灰度图。

我们首先导入一些必要的模块

import tensorflow as tf
import matplotlib.pyplot as plt
import input_data
mnist=input_data.read_data_sets("MNIST_data/",one_hot=True)

这里input_data.read_data_sets函数生成的类会自动将MNIST数据集划分为train,validation和test三个数据集，其中train中有55000张图片，validation集合内有5000张图片，这两个集合合起来就是MNIST的训练数据集。剩下的test集合中有10000张图片，这是MNIST的测试数据集。我们的每张图片将被处理成一个长度784（28281）一维数组。

解释一下one_hot,我们最后的分类有10个数字，我们采用二分类的模型，one_hot=True就是让其中的一个元素为1，其余元素为0，这样就将10分类简化为2分类啦~

我们可以打印MNIST数据集中的一张图片来看看效果

img = mnist.test.images[0]
image = img.reshape(28,28)
plt.imshow(image)
plt.show()

这里我们打印test数据集中的第一张图片，打印之前我们先通过reshape将其还原为一个28*28的像素矩阵，然后用matplotlib把他打印出来

大概就是这个样子，但是其实际为一个灰度图，只有黑白两种颜色

实现SLogistic Regression

"""定义出一些参数"""
training_num = 20 #总共训练几轮
learning_rate = 0.01 #学习率
batch_size = 100 #每次取一部分图片，随机梯度上升

"""mnist Graph input"""
x  = tf.placeholder("float",shape=[None,784])
y_ = tf.placeholder("float",shape=[None,10])

"""定义权重矩阵W和偏置b"""
W = tf.Variable(tf.zeros([784,10]))
b = tf.Variable(tf.zeros([10]))

"""softmax计算输出"""
y = tf.nn.softmax(tf.matmul(x,W) + b)

"""用交叉熵作为损失函数"""
cost = tf.reduce_mean(-tf.reduce_sum(y_*tf.log(y),reduction_indices=1))

"""梯度下降优化器"""
train_step = tf.train.GradientDescentOptimizer(learning_rate).minimize(cost)

我们定义了一些参数，这些参数在后面的过程中会用到，写在前面是为了在环境改变时能快速的做相应的调整。这里的reduce_sum根据函数名很容易知道是求和，而外面的的reduce_mean自然是求平均数，这里解释一下reduction_indices。这个参数表示函数处理的维度，若没有写则默认为None，即把输入的tensor降到0维，也就是一个数。

placeholder是一个占位符，在后面计算时需要给他相应的值。

创建Session

注意一定要有Session！！！不然上面所有的工作都只是创建了一个graph，实际并不做任何计算。我们需要用Session来计算图或图的一部分。

在开始之前还要对所有的变量进行初始化。

这里的batch_size每次取训练集中的一小部分，这样可以一定程度上增大最后到达全局最优的可能性，也相当于一个随机的梯度下降。我们用feed_dict将数据“喂给”前文所提到的占位符。

"""初始化所有变量"""
init = tf.initialize_all_variables()
with tf.Session() as sess:
    sess.run(init)
    #训练20轮
    for lun in range(training_num):
        avg_cost = 0;
        num = int(mnist.train.num_examples/batch_size)
        for i in range(num):
            batch = mnist.train.next_batch(batch_size)
            c=sess.run([train_step,cost],feed_dict={
    x:batch[0],y_:batch[1]})
            avg_cost += c[1]/num
        print("after %d training,the cost is %.9f" %(lun+1,avg_cost))

    print("training over!")
    """测试其准确率"""
    correct_prediction = tf.equal(tf.argmax(y, 1), tf.argmax(y_, 1))
    accuracy = tf.reduce_mean(tf.cast(correct_prediction, "float"))
    print("Accuracy:",sess.run(accuracy, feed_dict={
    x: mnist.test.images, y_: mnist.test.labels}))

准确率

可以看到我们最后的准确率是91%左右，这个结果怎么样呢？很抱歉，91的准确率并不是很理想，不过也别灰心，这是我们的第一个算法，也是灰常灰常简单的一种，后面还有很多美妙神奇的算法等着我们去学习，加油吧，骚年！

风车是个渣渣，如果文章有什么错误，欢迎指出~

参考文献

周志华《机器学习》

李锐《机器学习实战》