机器学习-Logistic回归

最新推荐文章于 2024-05-28 09:28:30 发布

joezarlove88

最新推荐文章于 2024-05-28 09:28:30 发布

阅读量861

点赞数

文章标签：算法机器学习 sklearn

本文链接：https://blog.csdn.net/joezarlove88/article/details/121460262

版权

前言：

逻辑回归是应用非常广泛的一个分类机器学习算法，它将数据拟合到一个logit函数(或者叫做logistic函数)中，从而能够完成对事件发生的概率进行预测。要说逻辑回归，我们得追溯到线性回归，想必大家对线性回归都有一定的了解，即对于多维空间中存在的样本点，我们用特征的线性组合去拟合空间中点的分布和轨迹。如下图所示：

线性回归能对连续值结果进行预测，而现实生活中常见的另外一类问题是，分类问题。最简单的情况是是与否的二分类问题。比如说医生需要判断病人是否生病，银行要判断一个人的信用程度是否达到可以给他发信用卡的程度，邮件收件箱要自动对邮件分类为正常邮件和垃圾邮件等等。很多实际的情况下，我们需要学习的分类数据并没有这么精准，比如说上述例子中突然有一个不按套路出牌的数据点出现.，而现实生活的分类问题的数据，会更为复杂，而这个时候我们借助于线性回归阈值的方式，已经很难完成一个鲁棒性很好的分类器了。

在这样的场景下，逻辑回归就诞生了。它的核心思想是，如果线性回归的结果输出是一个连续值，而值的范围是无法限定的，那我们有没有办法把这个结果值映射为可以帮助我们判断的结果呢。而如果输出结果是 (0,1) 的一个概率值，这个问题就很清楚了。我们在数学上找了一圈，还真就找着这样一个简单的函数了，就是很神奇的sigmoid函数(如下)：

函数图像如下：

从函数图上可以看出，函数y=g(z)在z=0的时候取值为1/2，而随着z逐渐变小，函数值趋于0，z逐渐变大的同时函数值逐渐趋于1，而这正是一个概率的范围。

所以我们定义线性回归的预测函数为Y=WTX，那么逻辑回归的输出Y= g(WTX)，其中y=g(z)函数正是上述sigmoid函数(或者简单叫做S形函数)。

1. 判定边界

在逻辑回归中，我们预测：

当hø大于等于0.5时，预测y=1
当hø小于0.5时，预测y=0

根据上面的预测，我们绘制出一条S形函数，如下：

又有：

所以

以上，为我们预知的逻辑回归的部分内容。好，现在假设我们有一个模型：并且参数ø是向量：[-3 1 1]。那么当-3+x1+x2大于等于0，即x1+x2大于等于3时，模型将预测 y=1。

我们可以绘制出来x1+x2=3，这条线便是我们模型的分界线，也称之为判定边界(Decision Boundary)，将预测为1的区域和预测为0的区域分隔开。

假设我们的数据呈现出如下图的分布情况，那么我们的模型是什么样才能适合这些数据呢？

如上图，函数图像为一个圆，圆点在原点且半径为1，这样一条曲线来分隔开了 y=1 和 y=0 的区域，所以我们需要的是一个二次方特征：

假设参数为 [-1 0 0 1 1]，则我们得到的判定边界恰好是圆点在原点并且半径为1的圆形。

我们可以使用非常复杂的模型来适应非常复杂形状的判定边界。

2.代价函数

逻辑回归模型的代价函数(Cost Function)

对于线性回归模型，我们定义的代价函数是所有模型误差的平方和。理论上讲，我们也可以沿用这个定义来对逻辑回归模型使用，但是问题在于，当我们将：

代入到这样定义的代价函数中时，我们得到的代价函数将会是一个非凸函数(Non-covex Function)。

这意味着，我们的代价函数将会有许多的局部最小值，这就会影响到梯度下降算法去找寻全局最小值。

因此，我们重新定义逻辑回归的代价函数为：

其中，Cost(hø(x(i), y(i))) 是我们定义的一个代价函数迭代形式，具体表示如下：

hø(x) 与 Cost(hø(x),y)之间的关系是如下图所示：

通过这样构建的Cost(hø(x), y)函数的特点是：

当实际的 y=1 且 hø=1 时，误差为0；当 y=1 但 hø != 1时，误差随hø的变小而变大；

当实际的 y=0 且 hø=0 时，误差代价为0；当 y=0 但 hø != 0 时，误差随hø的变大而变大。

将构建的Cost(hø(x), y) 进行一个简化，可以得到如下简化公式：

这个简化其实是对上面Cost(hø(x), y) 的两种表达式的一次性结合。

将简化代入到代价函数，得到：

这便是逻辑回归模型的代价函数了。

3.梯度下降

在得到这样的一个代价函数之后，我们便可以使用梯度下降算法(Gradient Descent)来求得能够使代价函数最小的参数了。

逻辑回归的梯度下降与线性回归的梯度下降是不一样的

梯度下降推导

注：虽然得到的梯度下降算法表面上看上去与线性回归的梯度下降算法一样，但是这里的 ℎ𝜃(𝑥)=𝑔(𝜃𝑇𝑋)与线性回归中不同，所以实际上是不一样的。另外，在运行梯度下降算法之前，进行特征缩放依旧是非常必要的。

一些梯度下降算法之外的选择：除了梯度下降算法以外，还有一些常被用来令代价函数最小的算法，这些算法更加复杂和优越，而且通常不需要人工选择学习率，通常比梯度下降算法要更加快速。这些算法有：
共轭梯度 Conjugate Gradient 局部优化法 (Broyden fletcher goldfarb shann,BFGS)和有限内存局部优化法 (LBFGS)
fminunc是 matlab和 octave 中都带的一个最小值优化函数，使用时我们需要提供代价函数和每个参数的求导.

4.代码与实现

1.数据集准备

2.数据集二维坐标分布

from numpy import loadtxt,where
from pylab import scatter, show, legend, xlabel, ylabel
 
#load the dataset
data = loadtxt("F:/PythonCode/LogisticRegression/data2.txt", delimiter=",")
#可以看出数据是一个二维数组，维度是100*3
print(data)
 
X = data[:,0:2]
#X存放的是数据的特征，维度是：100*2
# print(X.shape)
y = data[:, 2]
#y存放的是数据的标签，维度是：100*1
# print(y)
 
pos = where(y == 1)
#pos是y中数据等于1的下标索引
# print(pos)
neg = where(y==0)
#neg是y中数据等于0的下标索引
# print(neg)
 
#python中数据可视化函数scatter(数据的横坐标向量，数据的纵坐标向量，marker='0'数据以点的形式显示，c='b'数据点是blue颜色)
scatter(X[pos,0],X[pos, 1],marker='o', c='b')
scatter(X[neg,0],X[neg, 1],marker='x', c='r')
 
#说明二维坐标中o表示Pass,x表示Fail
legend(["y==1","y==0"])
show()

3. 下面我们写好计算sigmoid函数、代价函数、和梯度下降的程序：

from numpy import *
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize
 
filename = "F:/PythonCode/LogisticRegression/data2.txt"
 
def loadDataSet():
    # load the dataset
    data = loadtxt("F:/PythonCode/LogisticRegression/data2.txt", delimiter=",")
    # 拿到X和y
    y = np.c_[data[:, 2]]
    X = data[:, 0:2]
    return data,X,y
 
def map_feature(x1, x2):
    '''''
    Maps the two input features to polonomial features.
    Returns a new feature array with more features of
    X1, X2, X1 ** 2, X2 ** 2, X1*X2, X1*X2 ** 2, etc...
    '''
    x1.shape =(x1.size,1)
    x2.shape =(x2.size,1)
    degree =6
    mapped_fea = ones(shape=(x1[:,0].size,1))
    for i in range(1, degree +1):
        for j in range(i +1):
            r =(x1 **(i - j))*(x2 ** j)
            mapped_fea = append(mapped_fea, r, axis=1)
    return mapped_fea
 
#计算Sigmoid函数
def sigmoid(X):
    '''Compute sigmoid function'''
    den = 1.0 + exp(-1.0*X)
    gz = 1.0/den
    return gz
 
# 定义损失函数
def costFunctionReg(theta, X, y, l):
    m = y.size
    h = sigmoid(X.dot(theta))
 
    J = -1.0 * (1.0 / m) * (np.log(h).T.dot(y) + np.log(1 - h).T.dot(1 - y)) + (l / (2.0 * m)) * np.sum(np.square(theta[1:]))
 
    if np.isnan(J[0]):
        return (np.inf)
    return (J[0])
 
#计算梯度
def compute_grad(theta, X, y, l):
    m = y.size
    h = sigmoid(X.dot(theta.reshape(-1, 1)))
 
    grad = (1.0 / m) * X.T.dot(h - y) + (l / m) * np.r_[[[0]], theta[1:].reshape(-1, 1)]
 
    return (grad.flatten())
 
#梯度下降并优化
def gradAscent(XX, y, l):
    initial_theta = np.zeros(XX.shape[1])
    cost = costFunctionReg(initial_theta, XX, y, l)
    print('Cost: \n', cost)
    # 最优化 costFunctionReg
    res2 = minimize(costFunctionReg, initial_theta, args=(XX, y, l), jac=compute_grad, options={'maxiter': 3000})
    return res2
 
 
def plotBestFit(data,res2,X,accuracy,l,axes):  #画出最终分类的图
    # 对X,y的散列绘图
    plotData(data, 'Microchip Test 1', 'Microchip Test 2', 'y = 1', 'y = 0', axes=None)
    # 画出决策边界
    x1_min, x1_max = X[:, 0].min(), X[:, 0].max(),
    x2_min, x2_max = X[:, 1].min(), X[:, 1].max(),
    xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))
    h = sigmoid(map_feature(xx1.ravel(), xx2.ravel()).dot(res2.x))
    h = h.reshape(xx1.shape)
    if axes == None:
        axes = plt.gca()
    axes.contour(xx1, xx2, h, [0.5], linewidths=1, colors='g');
    axes.set_title('Train accuracy {}% with Lambda = {}'.format(np.round(accuracy, decimals=2), l))
    plt.show()
 
def plotData(data, label_x, label_y, label_pos, label_neg, axes):
    # 获得正负样本的下标(即哪些是正样本，哪些是负样本)
    neg = data[:, 2] == 0
    pos = data[:, 2] == 1
    if axes == None:
        axes = plt.gca()
    axes.scatter(data[pos][:, 0], data[pos][:, 1], marker='+', c='k', s=60, linewidth=2, label=label_pos)
    axes.scatter(data[neg][:, 0], data[neg][:, 1], c='y', s=60, label=label_neg)
    axes.set_xlabel(label_x)
    axes.set_ylabel(label_y)
    axes.legend(frameon=True, fancybox=True)
 
def predict(theta, X):
    '''''Predict whether the label
    is 0 or 1 using learned logistic
    regression parameters '''
    m, n = X.shape
    p = zeros(shape=(m,1))
    h = sigmoid(X.dot(theta.T))
    for it in range(0, h.shape[0]):
        if h[it]>0.5:
            p[it,0]=1
        else:
            p[it,0]=0
    return p
 
def main():
    data, X, y = loadDataSet()
    #对给定的两个feature做一个多项式特征的映射
    mapped_fea = map_feature(X[:, 0], X[:, 1])
 
    # 决策边界，咱们分别来看看正则化系数lambda太大太小分别会出现什么情况
    # Lambda = 0 : 就是没有正则化，这样的话，就过拟合咯
    # Lambda = 1 : 这才是正确的打开方式
    # Lambda = 100 :正则化项太激进，导致基本就没拟合出决策边界
    l = 1
 
    res = gradAscent(mapped_fea, y, l)
    print(res)
 
    # 准确率
    accuracy = y[where(predict(res.x, mapped_fea) == y)].size / float(y.size)*100.0
    #画决策边界
    plotBestFit(data, res, X,  accuracy, l,axes=None)
 
if __name__ == '__main__':
    main()

代码运行结果如下: