Python数据分析与机器学习-Python实现逻辑回归与梯度下降策略

源码下载：

http://download.csdn.net/download/adam_zs/10181481

'''
Logistic Regression
The data
我们将建立一个逻辑回归模型来预测一个学生是否被大学录取。
假设你是一个大学系的管理员，你想根据两次考试的结果来决定每个申请人的录取机会。
你有以前的申请人的历史数据，你可以用它作为逻辑回归的训练集。
对于每一个培训例子，你有两个考试的申请人的分数和录取决定。
为了做到这一点，我们将建立一个分类模型，根据考试成绩估计入学概率。
'''
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os

path = 'data' + os.sep + 'LogiReg_data.txt'
pdData = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])
# print(pdData.head())
'''
      Exam 1     Exam 2  Admitted
0  34.623660  78.024693         0
1  30.286711  43.894998         0
2  35.847409  72.902198         0
3  60.182599  86.308552         1
4  79.032736  75.344376         1
'''
# print(pdData.shape)
'''
(100, 3)
'''
positive = pdData[pdData["Admitted"] == 1]  # 正例
negative = pdData[pdData["Admitted"] == 0]  # 负例
fig, ax = plt.subplots(figsize=(10, 5))
# s点的大小
ax.scatter(positive['Exam 1'], positive['Exam 2'], s=30, c="b", marker="o", label="Admitted")
ax.scatter(negative['Exam 1'], negative['Exam 2'], s=30, c="r", marker="x", label="Not Admitted")
ax.legend()
ax.set_xlabel("Exam 1 Score")
ax.set_ylabel("Exam 2 Score")
# plt.show()

'''
The logistic regression 逻辑回归
目标：建立分类器（求解出三个参数 θ0θ1θ2θ0θ1θ2）
设定阈值，根据阈值判断录取结果
要完成的模块
sigmoid : 映射到概率的函数
model : 返回预测结果值
cost : 根据参数计算损失
gradient : 计算每个参数的梯度方向
descent : 进行参数更新
accuracy: 计算精度
sigmoid 函数¶
    g(z)=1/(1+e−z)
'''


# 定义sigmoid函数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))


nums = np.arange(-10, 10, step=1)
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(nums, sigmoid(nums), "r")


# plt.show()

# 定义模型
def model(X, theta):
    return sigmoid(np.dot(X, theta.T))  # np.do矩阵乘法


pdData.insert(0, 'Ones', 1)  # loc:插入列的位置 column:列名 value:列值
orig_data = pdData.as_matrix()
cols = orig_data.shape[1]  # 列数
X = orig_data[:, 0:cols - 1]
y = orig_data[:, cols - 1:cols]
theta = np.zeros([1, 3])  # 构造3个参数
# print(X[0:5])
'''
[[  1.          34.62365962  78.02469282]
 [  1.          30.28671077  43.89499752]
 [  1.          35.84740877  72.90219803]
 [  1.          60.18259939  86.3085521 ]
 [  1.          79.03273605  75.34437644]]
'''
# print(y[0:5])
'''
[[ 0.]
 [ 0.]
 [ 0.]
 [ 1.]
 [ 1.]]
'''
# print(theta)
# print(X.shape, y.shape, theta.shape)
'''(100, 3) (100, 1) (1, 3)'''


def cost(X, y, theta):
    '''
    损失函数
    将对数似然函数去负号
    '''
    left = np.multiply(-y, np.log(model(X, theta)))
    right = np.multiply(1 - y, np.log(1 - model(X, theta)))
    return np.sum(left - right) / (len(X))


# print(cost(X, y, theta))
'''0.69314718056'''


def gradient(X, y, theta):
    '''计算梯度'''
    grad = np.zeros(theta.shape)
    error = (model(X, theta) - y).ravel()
    for j in range(len(theta.ravel())):  # for each parmeter
        term = np.multiply(error, X[:, j])
        grad[0, j] = np.sum(term) / len(X)
    return grad


'''
Gradient descent
比较3中不同梯度下降方法
1.批量梯度下降
2.随机梯度下降
3.小批量梯度下降法
'''

# 3种停止策略
STOP_ITER = 0  # 按照迭代次数停止
STOP_COST = 1  # 两次迭代之间的损失值差异
STOP_GRAD = 2  # 梯度


def stopCriterion(type, value, threshold):
    # 设定三种不同的停止策略
    if type == STOP_ITER:
        return value > threshold
    elif type == STOP_COST:
        return abs(value[-1] - value[-2]) < threshold
    elif type == STOP_GRAD:
        return np.linalg.norm(value) < threshold


import numpy.random


# 迭代之前洗牌
def shuffleData(data):
    np.random.shuffle(data)
    cols = data.shape[1]
    X = data[:, 0:cols - 1]
    y = data[:, cols - 1:]
    return X, y


import time


# data数据 theta参数 batchSize每次取样本数(1:随机梯度下降,最大值:批量梯度下降,1-最大值:小批量梯度下降法)
# stopType停止策略 thresh策略对应的域值 alpha学习率
def descent(data, theta, batchSize, stopType, thresh, alpha):
    # 梯度下降求解

    init_time = time.time()
    i = 0  # 迭代次数
    k = 0  # batch
    X, y = shuffleData(data)
    grad = np.zeros(theta.shape)  # 计算的梯度
    costs = [cost(X, y, theta)]  # 损失值

    while True:
        grad = gradient(X[k:k + batchSize], y[k:k + batchSize], theta)
        k += batchSize  # 取batch数量个数据
        if k >= n:
            k = 0
            X, y = shuffleData(data)  # 重新洗牌
        theta = theta - alpha * grad  # 参数更新
        costs.append(cost(X, y, theta))  # 计算新的损失
        i += 1

        if stopType == STOP_ITER:
            value = i
        elif stopType == STOP_COST:
            value = costs
        elif stopType == STOP_GRAD:
            value = grad
        if stopCriterion(stopType, value, thresh): break

    return theta, i - 1, costs, grad, time.time() - init_time


def runExpe(data, theta, batchSize, stopType, thresh, alpha):
    # import pdb; pdb.set_trace();
    theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, alpha)
    # 画图展示
    name = "Original" if (data[:, 1] > 2).sum() > 1 else "Scaled"
    name += " data - learning rate: {} - ".format(alpha)
    if batchSize == n:
        strDescType = "Gradient"
    elif batchSize == 1:
        strDescType = "Stochastic"
    else:
        strDescType = "Mini-batch ({})".format(batchSize)
    name += strDescType + " descent - Stop: "
    if stopType == STOP_ITER:
        strStop = "{} iterations".format(thresh)
    elif stopType == STOP_COST:
        strStop = "costs change < {}".format(thresh)
    else:
        strStop = "gradient norm < {}".format(thresh)
    name += strStop
    print("***{}\nTheta: {} - Iter: {} - Last cost: {:03.2f} - Duration: {:03.2f}s".format(
        name, theta, iter, costs[-1], dur))
    fig, ax = plt.subplots(figsize=(12, 4))
    ax.plot(np.arange(len(costs)), costs, 'r')
    ax.set_xlabel('Iterations')
    ax.set_ylabel('Cost')
    ax.set_title(name.upper() + ' - Error vs. Iteration')
    return theta


'''
不同的停止策略
设定迭代次数
'''
# 选择的梯度下降方法是基于所有样本的
n = 100  # 批量梯度下降,STOP_ITER按照迭代次数停止,迭代5000次
runExpe(orig_data, theta, n, STOP_ITER, thresh=5000, alpha=0.000001)
# plt.show()
'''
根据损失值停止
设定阈值 1E-6, 差不多需要110 000次迭代
'''
# 两次迭代之间的损失值差异
runExpe(orig_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001)
# plt.show()

'''
根据梯度变化停止
设定阈值 0.05,差不多需要40 000次迭代
'''
runExpe(orig_data, theta, n, STOP_GRAD, thresh=0.05, alpha=0.001)

'''
对比不同的梯度下降方法
'''
# 随机梯度下降
runExpe(orig_data, theta, 1, STOP_ITER, thresh=5000, alpha=0.001)

# 有点爆炸。。。很不稳定,再来试试把学习率调小一些
# 随机梯度下降 调小学习率
runExpe(orig_data, theta, 1, STOP_ITER, thresh=15000, alpha=0.000002)
# 速度快，但稳定性差，需要很小的学习率

# 小批量梯度下降法
runExpe(orig_data, theta, 16, STOP_ITER, thresh=15000, alpha=0.001)

'''
浮动仍然比较大，我们来尝试下对数据进行标准化 将数据按其属性(按列进行)减去其均值，然后除以其方差。最后得到的结果是，对每个属性/每列来说所有数据都聚集在0附近，方差值为1
'''
from sklearn import preprocessing as pp

# 数据标准化
scaled_data = orig_data.copy()
scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3])

runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.001)

# 它好多了！原始数据，只能达到达到0.61，而我们得到了0.38个在这里！ 所以对数据做预处理是非常重要的

runExpe(scaled_data, theta, n, STOP_GRAD, thresh=0.02, alpha=0.001)
# 更多的迭代次数会使得损失下降的更多！
theta = runExpe(scaled_data, theta, 1, STOP_GRAD, thresh=0.002 / 5, alpha=0.001)

# 随机梯度下降更快，但是我们需要迭代的次数也需要更多，所以还是用batch的比较合适！！！
runExpe(scaled_data, theta, 16, STOP_GRAD, thresh=0.002 * 2, alpha=0.001)

'''精度'''


# 设定阈值
def predict(X, theta):
    return [1 if x >= 0.5 else 0 for x in model(X, theta)]


scaled_X = scaled_data[:, :3]
y = scaled_data[:, 3]
predictions = predict(scaled_X, theta)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
accuracy = (sum(map(int, correct)) % len(correct))
print('accuracy = {0}%'.format(accuracy))