《人工智能》机器学习 - 第5章 逻辑回归（三 非线性逻辑回归实战）

5.3非线性逻辑回归实战

前文讲解了线性的逻辑回归，那么对于非线性的是否使用呢？答案是肯定的。只是在开始训练之前稍微处理一下数据。代码如下。

# 生成阶数为6的多项式
poly = PolynomialFeatures(6)
XX = poly.fit_transform(X[:,1:3])

另外在画边界图和线性的还有些不一样，具体请看代码。

# -*- coding:utf-8 -*-
import numpy as np
import time
import matplotlib.pyplot as plt
import matplotlib
from sklearn.preprocessing import PolynomialFeatures

matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False

"""
函数说明:读取数据

Parameters:
    filename - 文件名
Returns:
    xArr - x数据集
    yArr - y数据集
"""
def loadDataSet(filename):
    X = []
    Y = []
    
    with open(filename, 'rb') as f:
        for idx, line in enumerate(f):
            line = line.decode('utf-8').strip()
            if not line:
                continue
            eles = line.split()
            eles = list(map(float, eles))
            
            if idx == 0:
                numFea = len(eles)
            #去掉每行的最后一个数据再放入X中
            X.append(eles[:-1])
            #获取每行的租后一个数据
            Y.append([eles[-1]])
     
    # 将X,Y列表转化成矩阵
    xArr = np.array(X)
    yArr = np.array(Y)
    
    return xArr,yArr

"""
函数说明:定义sigmoid函数

Parameters:
    z
Returns:
    返回sigmoid函数
"""
def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

"""
函数说明:定义代价函数

Parameters:
    theta, X, Y, theLambda
Returns:
    返回
"""
def J(theta, X, Y, theLambda=0):
    m, n = X.shape
    h = sigmoid(np.dot(X, theta))
    J = (-1.0/m)*(np.dot(np.log(h).T,Y)+np.dot(np.log(1-h).T,1-Y)) + (theLambda/(2.0*m))*np.sum(np.square(theta[1:]))
    if np.isnan(J[0]):
        return np.inf
    # 其实J里面只有一个数值，需要取出该数值
    return J.flatten()[0]

"""
函数说明:梯度下降求解参数

Parameters:
    X, Y, options
Returns:
    返回参数
"""
def gradient(X, Y, options):
    '''
    options.alpha 学习率
    options.theLambda 正则参数λ
    options.maxLoop 最大迭代轮次
    options.epsilon 判断收敛的阈值
    options.method
        - 'sgd' 随机梯度下降
        - 'bgd' 批量梯度下降
    '''
    
    m, n = X.shape
    
    # 初始化模型参数，n个特征对应n个参数
    theta = np.zeros((n,1))
    
    cost = J(theta, X, Y)  # 当前误差
    costs = [cost]
    thetas = [theta]
    
    # Python 字典dict.get(key, default=None)函数返回指定键的值，如果值不在字典中返回默认值
    alpha = options.get('alpha', 0.01)
    epsilon = options.get('epsilon', 0.00001)
    maxloop = options.get('maxloop', 1000)
    theLambda = float(options.get('theLambda', 0)) # 后面有 theLambda/m 的计算，如果这里不转成float，后面这个就全是0
    method = options.get('method', 'bgd')
    
    # 定义随机梯度下降
    def _sgd(theta):
        count = 0
        converged = False
        while count < maxloop:
            if converged :
                break
            # 随机梯度下降，每一个样本都更新
            for i in range(m):
                h =sigmoid(np.dot(X[i].reshape((1,n)), theta))
                
                theta = theta - alpha*((1.0/m)*X[i].reshape((n,1))*(h-Y[i]) + (theLambda/m)*np.r_[[[0]], theta[1:]])
                thetas.append(theta)
                cost = J(theta, X, Y, theLambda)
                costs.append(cost)
                if abs(costs[-1] - costs[-2]) < epsilon:
                    converged = True
                    break
            count += 1
        return thetas, costs, count
    
    # 定义批量梯度下降
    def _bgd(theta):
        count = 0
        converged = False
        while count < maxloop:
            if converged :
                break
            
            h = sigmoid(np.dot(X, theta))
            theta = theta - alpha*((1.0/m)*np.dot(X.T, (h-Y)) + (theLambda/m)*np.r_[[[0]],theta[1:]])
            
            thetas.append(theta)
            cost = J(theta, X, Y, theLambda)
            costs.append(cost)
            if abs(costs[-1] - costs[-2]) < epsilon:
                converged = True
                break
            count += 1
        return thetas, costs, count
    
    methods = {'sgd': _sgd, 'bgd': _bgd}
    return methods[method](theta)  

"""
函数说明:绘制决策边界

Parameters:
    X, thetas
Returns:
    无
"""
def plotBoundary(X,Y,thetas,poly):
        
    # 绘制决策边界
    plt.figure(figsize=(6,4))
    for i in range(len(X)):
        x = X[i]
        if Y[i] == 1:
            plt.scatter(x[1], x[2], marker='*', color='blue', s=50)
        else:
            plt.scatter(x[1], x[2], marker='o', color='green', s=50)
    
    # 绘制决策边界
    x1Min,x1Max,x2Min,x2Max = X[:, 1].min(), X[:, 1].max(), X[:, 2].min(), X[:, 2].max()
    xx1, xx2 = np.meshgrid(np.linspace(x1Min, x1Max), np.linspace(x2Min, x2Max))
    h = sigmoid(poly.fit_transform(np.c_[xx1.ravel(), xx2.ravel()]).dot(thetas[-1]))
    h = h.reshape(xx1.shape)
    plt.contour(xx1, xx2, h, [0.5], colors='red')
    plt.xlabel(r'$x_1$')
    plt.ylabel(r'$x_2$')

"""
函数说明:绘制代价函数图像

Parameters:
    costs
Returns:
    返回参数
"""
def plotCost(costs):
    # 绘制误差曲线
    plt.figure(figsize=(8,4))
    plt.plot(range(len(costs)), costs)
    plt.xlabel(u'迭代次数')
    plt.ylabel(u'代价J')    

"""
函数说明:绘制参数曲线

Parameters:
    thetass - 参数
Returns:
    返回参数
"""
def plotthetas(thetas):
    # 绘制参数theta变化
    thetasFig, ax = plt.subplots(len(thetas[0]))
    thetas = np.asarray(thetas)
    for idx, sp in enumerate(ax):
        thetaList = thetas[:, idx]
        sp.plot(range(len(thetaList)), thetaList)
        sp.set_xlabel('Number of iteration')
        sp.set_ylabel(r'$\theta_%d$'%idx)

if __name__ == '__main__':
    ## Step 1: load data...
    print("Step 1: load data...")
    ori_X, Y = loadDataSet('non_data.txt')
    m, n = ori_X.shape
    X = np.concatenate((np.ones((m,1)), ori_X), axis=1)
    
    # 生成阶数为6的多项式
    poly = PolynomialFeatures(6)
    XX = poly.fit_transform(X[:,1:3])
    m, n = XX.shape
        
    ## Step 2: training bgd...
    print("Step 2: training bgd...")
    
    # 训练数据
    options = {
        'alpha': 1.0,
        'epsilon': 0.000001,
        'theLambda': 0, # 100
        'maxloop': 10000,
        'method': 'bgd'
    }
    
    start = time.time()
    bgd_thetas, bgd_costs, bgd_iterationCount = gradient(XX, Y, options)
    print("bgd_thetas:")
    print(bgd_thetas[-1])
    end = time.time()
    bgd_time = end - start
    print("bgd_time:"+str(bgd_time))
    
    ## Step 3: show boundary bgd...
    print("Step 3: show boundary bgd...")
    plotBoundary(X,Y,bgd_thetas,poly)

    ## Step 4: show costs bgd...
    print("Step 4: show costs bgd...")
    plotCost(bgd_costs)

    ## Step 5: show thetas bgd...
    print("Step 5: show thetas bgd...")
    plotthetas(bgd_thetas)
   
    print("======================================================")
    ## Step 2: training sgd...
    print("Step 2: training sgd...")
    
    # sgd随机梯度上升
    options = {
        'alpha':0.5,
        'epsilon':0.00000001,
        'theLambda': 0, # 100
        'maxloop':100000,
        'method':'sgd' 
    }

    # 训练模型
    start = time.time()
    sgd_thetas, sgd_costs, sgd_iterationCount = gradient(XX, Y, options)
    print("sgd_thetas:")
    print(sgd_thetas[-1])
    end = time.time()
    sgd_time = end - start
    print("sgd_time:"+str(sgd_time))
    
    ## Step 3: show boundary sgd...
    print("Step 3: show boundary sgd...")
    plotBoundary(X,Y,sgd_thetas,poly)

    ## Step 4: show costs sgd...
    print("Step 4: show costs sgd...")
    plotCost(sgd_costs)

    ## Step 5: show thetas sgd...
    print("Step 5: show thetas sgd...")
    plotthetas(sgd_thetas)

结果太多，不想贴出来，有兴趣的自己运行吧。
【完整代码参考附件2.Non_Logistic_Regression_Classify\Non_Logistic_Regression_Classify_v1
Non_Logistic_Regression_Classify_v1.0.py】

 调用sklearn库实现

使用sklearn代码如下所示。

# -*- coding:utf-8 -*-
import numpy as np
import time
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import matplotlib
from sklearn.preprocessing import PolynomialFeatures

matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False

"""
函数说明:读取数据

Parameters:
    filename - 文件名
Returns:
    xArr - x数据集
    yArr - y数据集
"""
def loadDataSet(filename):
    X = []
    Y = []
    
    with open(filename, 'rb') as f:
        for idx, line in enumerate(f):
            line = line.decode('utf-8').strip()
            if not line:
                continue
            eles = line.split()
            eles = list(map(float, eles))
            
            if idx == 0:
                numFea = len(eles)
            #去掉每行的最后一个数据再放入X中
            X.append(eles[:-1])
            #获取每行的租后一个数据
            Y.append([eles[-1]])
     
    # 将X,Y列表转化成矩阵
    xArr = np.array(X)
    yArr = np.array(Y)
    
    return xArr,yArr

"""
函数说明:定义sigmoid函数

Parameters:
    z
Returns:
    返回sigmoid函数
"""
def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))


"""
函数说明:绘制决策边界

Parameters:
    X, thetas
Returns:
    无
"""
def plotBoundary(X,Y,thetas,poly):
        
    # 绘制决策边界
    plt.figure(figsize=(6,4))
    for i in range(len(X)):
        x = X[i]
        if Y[i] == 1:
            plt.scatter(x[1], x[2], marker='*', color='blue', s=50)
        else:
            plt.scatter(x[1], x[2], marker='o', color='green', s=50)
    
    # 绘制决策边界
    x1Min,x1Max,x2Min,x2Max = X[:, 1].min(), X[:, 1].max(), X[:, 2].min(), X[:, 2].max()
    xx1, xx2 = np.meshgrid(np.linspace(x1Min, x1Max), np.linspace(x2Min, x2Max))
    h = sigmoid(poly.fit_transform(np.c_[xx1.ravel(), xx2.ravel()]).dot(thetas))
    h = h.reshape(xx1.shape)

    plt.contour(xx1, xx2, h, [0.5], colors='red')
    plt.xlabel(r'$x_1$')
    plt.ylabel(r'$x_2$')



"""
函数说明:绘制预测值和实际值比较图

Parameters:
    X_test,Y_test,y_predict
Returns:
    无
"""
def plotfit(X_test,Y_test,y_predict):
    
    p_x = range(len(X_test))
    
    plt.figure(figsize=(6,4),facecolor="w")
    plt.ylim(0,6)
    plt.plot(p_x,Y_test,"ro",markersize=8,zorder=3,label=u"真实值")
    plt.plot(p_x,y_predict,"go",markersize=14,zorder=2,label=u"预测值,$R^2$=%.3f" %lr.score(X_test,Y_test))
    plt.legend(loc="upper left")
    plt.xlabel(u"编号",fontsize=12)
    plt.ylabel(u"类型",fontsize=12)
    plt.title(u"Logistic算法对数据进行分类",fontsize=12)
    #plt.savefig("Logistic算法对数据进行分类.png")
    plt.show()
    
#测试
if __name__ == '__main__':
    ## Step 1: load data...
    print("Step 1: load data...")

    ori_X, Y = loadDataSet('non_data.txt')
    m, n = ori_X.shape
    X = np.concatenate((np.ones((m,1)), ori_X), axis=1)
    
    # 生成阶数为6的多项式
    poly = PolynomialFeatures(6)
    XX = poly.fit_transform(X[:,1:3])
    m, n = XX.shape
    #print(XX)
    #print(Y)
        
    ## Step 2: split data...
    print("Step 2: split data...")
    X_train,X_test,Y_train,Y_test=train_test_split(XX,Y,test_size=0.2,random_state=0)

    start = time.time()
    ## Step 3: init LogisticRegression...
    print("Step 3: init LogisticRegression...")
    lr=LogisticRegression(C=1000.0,random_state=0)
 
    ## Step 4: training...
    print("Step 4: training ...")
    #输出y做出ravel()转换。不然会有警告
    lr.fit(X_train,Y_train.ravel())
    end = time.time()
    sk_time = end - start
    print('time:',sk_time)
    
    ## Step 5: show thetas...
    print("Step 5: show thetas...")
    #模型效果获取
    r = lr.score(X_train,Y_train)
    print("R值(准确率):",r)
       
    print("参数:",lr.coef_)
    print("截距:",lr.intercept_)
 
    w0 = np.array(lr.intercept_)
    w1 = np.array(lr.coef_.T[1:28])
    
    #合并为一个矩阵
    thetas = np.vstack((w0,w1))
    print('thetas:')
    print(thetas)
        
    ## Step 6: show sigmoid p...
    print("Step 6: show sigmoid p...")
    #sigmoid函数转化的值，即：概率p
    predict_p = lr.predict_proba(X_test)
    print(predict_p)     #sigmoid函数转化的值，即：概率pprint(y_predict)     #sigmoid函数转化的值，即：概率p
    
    ## Step 7: show boundary...
    print("Step 7: show boundary...")
    plotBoundary(X,Y,thetas,poly)
    
    ## Step 8: predict ...
    print("Step 8: predict...")
    y_predict = lr.predict(X_test) 
    print(y_predict)
    
    ## Step 9: show plotfit ...
    print("Step 9: show plotfit...")
    #画图对预测值和实际值进行比较
    plotfit(X_test,Y_test,y_predict)

结果如下所示。

【完整代码参考附件2.Non_Logistic_Regression_Classify\Non_Logistic_Regression_Classify-sklearn_v2\Non_Logistic_Regression_Classify-sklearn_v2.0.py】

本章参考代码
点击进入下载