ML——手写logistic-regression 实现线性二元分类器

参考李宏毅老师ML作业二：https://colab.research.google.com/drive/1JaMKJU7hvnDoUfZjvUKzm9u-JLeX6B2C#scrollTo=RJuoQ_R2jUmX

　　　　　　　　　　　https://github.com/PoRuiHuang/ML2020SPRING

"""
2020/4/11
cfancy
logistic-regression
实做一个线性二元分类器，来根据人们的个人资料，判断其年收入
是否高于50000$
所用数据是已经做了稍微的处理
X_train ：训练集
Y_trian ：目标集(训练标签集)
X_trian ：测试集
手写  gradient descent
"""
# 导入库
import sys
import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
import csv

#常用的基础函数
 #平均值
def mean(x):
  return sum(x) / len(x)

 #值与平均值之差
def de_mean(x):
  x_bar = mean(x)
  return [x_i - x_bar for x_i in x]

 #矩阵点乘
def dot(v, w):
    return sum(v_i * w_i for v_i, w_i in zip(v, w))
 #python 中 zip()函数：打包为元组的列表

 #
def sum_of_squares(v):
    """
    :param v:
    :return:
    """
    return dot(v, v)

 # 方差
def variance(x):
    """
    方差
    :param x:
    :return:
    """
    n = len(x)
    deviations = de_mean(x)
    return sum_of_squares(deviations) / (n - 1)
 #标准差
def standard_deviation(x):
    """
    标准差
    :param x:
    :return:
    """
    return math.sqrt(variance(x))

 #协方差
def covariance(x,y):
    """
    协方差
    :param x:
    :param y:
    :return:
    """
    n = len(x)
    return dot(de_mean(x), de_mean(y)) / (n - 1)

 #相关系数
def correlation(x, y):
    """
    相关系数
    :return:
    """
    stdev_x = standard_deviation(x)
    stdev_y = standard_deviation(y)
    if stdev_x > 0 and stdev_y > 0:
        return covariance(x, y) / stdev_x / stdev_y
    else:
        return 0

# 一些有用的函数
 # 打乱索引
def _shuffle(X, Y):
    """
    将两个长度相等的列表或者数组打乱其索引
    打乱训练集（即打乱索引）
    :param X:
    :param Y:
    :return:
    """
    randomize = np.arange(len(X))
    np.random.shuffle(randomize)
    return (X[randomize], Y[randomize])
# A = np.arange(10)
# B = np.arange(10)
#
# A ,B = _shuffle(A,B)
# print(A)  #[3 5 2 9 6 7 0 8 4 1]
#
# print(B)  #[3 5 2 9 6 7 0 8 4 1]

 #
def _sigmoid(z):
    """
    Sigmoid 函数用来计算概率
    :param z:
    :return:
    """
    return np.clip(1 / (1.0 + np.exp(-z)), 1e-8, 1 - (1e-8))  #clip 限制数组中的值

# activate function
def _f(X, w, b):
    """
    logistic_regression 函数
    :param X: 输入数据  shape = [batch_size, data_dimension
    :param w: 参数 权重 w 向量  shape = [data_dimension, ]
    :param b: 参数 偏置 b 数值
    :return: predicted probability of each row of X being positively labeled, shape = [batch_size, ]
    """
    return _sigmoid(np.matmul(X, w) + b)

def _predict(X, w, b):
    """
    结果通过四舍五入的
    :param X:
    :param w:
    :param b:
    :return:返回预测的值
    """
    return np.round(_f(X, w,b)).astype(np.int)

def _accuracy(Y_pred, Y_label):
    # This function calculates prediction accuracy
    acc = 1 - np.mean(np.abs(Y_pred - Y_label))
    return acc


"""--------------------准备数据--------------------"""
np.random.seed(0)   ##随机数生成器
# 读取数据 X_train Y_train X_test  并对数据进行正规化（归一化）normalization
X_train_fpath = './data/X_train'  #训练集
Y_train_fpath = './data/Y_train'  #目标集（标签集）
X_test_fpath = './data/X_test'  # 测试集
output_fpath = './output_{}.csv'  # 预测保存

#解析csv文件转为numpy数组
drop_list = []
with open(Y_train_fpath) as f:
    # next(f)
    # Y_train = np.array([line.strip('
').split(',')[1] for line in f], dtype = float)
    # print(Y_train.shape)
    # print(type(Y_train))
    # print(Y_train)

    data2 = pd.read_csv(f)
    # print(data2.shape)
    data2 = data2.iloc[:, 1:]
    # data2.drop(columns=['migration code-change in msa','migration code-change in reg','migration code-move within reg','migration prev res in sunbelt','country of birth father','country of birth mother','country of birth self'],inplace=True)

    # next(f)
    Y_train = data2.to_numpy(dtype=float)
    Y_train = Y_train.reshape((len(Y_train),))
# Parse csv files to numpy array
with open(X_train_fpath) as f:
    # next(f)
    # X_train = np.array([line.strip('
').split(',')[1:] for line in f], dtype = float)
    # print(X_train.shape)
    # print(X_train)

    data1 = pd.read_csv(f)
    # print(data1.shape)

    data1 = data1.iloc[:, 1:]
    # data1.drop(columns=['migration code-change in msa','migration code-change in reg','migration code-move within reg','migration prev res in sunbelt','country of birth father','country of birth mother','country of birth self'],inplace=True)
    # data.drop(columns=['id', 'migration code-change in msa','migration code-change in reg','migration code-move within reg','migration prev res in sunbelt','country of birth father','country of birth mother','country of birth self'],inplace=True)
    X_train = data1.to_numpy(dtype=float)

    for i in range(len(X_train[0])):
        corr = correlation(X_train[:, i], Y_train)
        # print(i,corr,sep=": ")
        if (abs(corr) < 0.01):
            drop_list.append(i)

    print(drop_list)
    # print(len(drop_list))
    print(X_train.shape)
    X_train = np.delete(X_train, drop_list, 1)


with open(X_test_fpath) as f:
    # next(f)
    # X_test = np.array([line.strip('
').split(',')[1:] for line in f], dtype = float)
    # print(X_test.shape)
    # print(X_test)

    data3 = pd.read_csv(f)
    # print(data3.shape)
    data3 = data3.iloc[:, 1:]
    # data3.drop(columns=['migration code-change in msa','migration code-change in reg','migration code-move within reg','migration prev res in sunbelt','country of birth father','country of birth mother','country of birth self'],inplace=True)
    X_test = data3.to_numpy(dtype=float)
    X_test = np.delete(X_test, drop_list, 1)
print("------------------------")
print(X_train)
print("------------------------")
print(Y_train)
print("------------------------")
print(X_test)

#正规化（归一化）
def _normalize(X, train = True , specified_column = None, X_mean = None, X_std = None):
    """
    此函数归一化X的特定列
    :param X: 数据
    :param train:  train = true 表示 处理training data,False 表示 testing data
    :param specified_column: 将被标准化的列的索引  None 表示所有列都被归一化
    :param X_mean: 当train = False 时training data 的平均值
    :param X_std:当train = False 时training data 的标准差
    :return: 返回归一化后的数据，training data的平均值 ，training data 的标准差
    处理测试数据时，将重用训练数据的均值和标准方差
    """
    if specified_column == None:
        specified_column = np.arange(X.shape[1])
    if train:
        X_mean = np.mean(X[:, specified_column], 0).reshape(1, -1)
        X_std   = np.std(X[:, specified_column], 0).reshape(1, -1)

    X[:, specified_column] = (X[:, specified_column] - X_mean) / (X_std + 1e-8)

    return X, X_mean, X_std

#将处理后的数据切分为训练集和发展集
def _train_dev_split(X, Y, dev_ratio = 0.25):
    """
    归一化后的数据切分为训练集与发展集
    :param X:
    :param Y:
    :param dev_ration:
    :return:
    """
    train_size = int(len(X) * (1 - dev_ratio))
    return X[:train_size], Y[:train_size], X[train_size:], Y[train_size:]

# 归一化训练集和测试集
X_train, X_mean, X_std = _normalize(X_train, train=True)
X_test, _, _ = _normalize(X_test, train=False, specified_column=None, X_mean=X_mean, X_std=X_std)

#将数据分为训练集和发展集
dev_ratio = 0.1
X_train, Y_train, X_dev, Y_dev = _train_dev_split(X_train, Y_train, dev_ratio = dev_ratio)

train_size = X_train.shape[0]
dev_size = X_dev.shape[0]
test_size = X_test.shape[0]
data_dim = X_train.shape[1]
print('Size of training set: {}'.format(train_size))
print('Size of development set: {}'.format(dev_size))
print('Size of testing set: {}'.format(test_size))
print('Dimension of data: {}'.format(data_dim))

"""-------------------梯度函数和损失函数-----------------------"""
def _cross_entropy_loss(y_pred, Y_label,):
    """
    交叉熵损失函数  根据推导的公式
    :param y_pred:  概率预测  浮点型向量
    :param Y_label: 真实标签值  布尔值向量
    :return: 返回的是数值
    """
    cross_entropy = -np.dot(Y_label, np.log(y_pred)) - np.dot((1 - Y_label), np.log(1 - y_pred))
    return cross_entropy

def _gradient(X, Y_label, w, b):
    """
    计算交叉熵的梯度  根据笔记推到的公式 w_grad
    :param X:
    :param Y_label:
    :param w:
    :param b:
    :return:
    """
    y_pred = _f(X, w, b)
    pred_error = Y_label - y_pred
    w_grad = -np.sum(pred_error * X.T, 1)
    b_grad = -np.sum(pred_error)
    return w_grad, b_grad

"""-------------------Training--------------------------"""
"""
使用小批次梯度下降法来训练（MBGD——Mini-batch gradient descent）,训练集被分为许多小批次，
我们分别计算其梯度以及损失，并根据该批次来更新模型的参数。当一次回圈(epoch)完成，也就是整个训练集的所有
小批次都被使用过一次后，我们将所有训练资料打散并且重新分成新的小的批次，进行下一个回圈，知道事先设定好的回圈数量为止。
"""
#初始化参数 w,b
w = np.zeros((data_dim,))
b = np.zeros((1,))

#训练时需要的参数 最大迭代次数、每个batch的大小以及学习率
max_iter = 40
batch_size = 10
lr = 0.1

#保存每次迭代时的损失和精度，以进行绘制
train_loss = []
dev_loss = []
train_acc = []
dev_acc = []

#计算参数更新的次数
count = 1

# 迭代训练
for epoch in range(max_iter):
    #在每个回圈开始时的随机变动
    X_train, Y_train = _shuffle(X_train, Y_train)

    #小批次训练
    for idx in range(int(np.floor(train_size / batch_size))):
        X = X_train[idx * batch_size:(idx+1) * batch_size]
        Y = Y_train[idx * batch_size:(idx+1) * batch_size]

        #计算梯度
        w_grad, b_grad = _gradient(X, Y, w, b)

        #梯度下降参数更新
        #学习率随时间衰减
        w = w - lr/np.sqrt(count) * w_grad
        b = b - lr/np.sqrt(count) * b_grad

        count = count + 1

    #计算训练集和发展集的损失和精度
    y_train_pred = _f(X_train, w, b)
    Y_train_pred = np.round(y_train_pred) #round ()四舍五入
    train_acc.append(_accuracy(Y_train_pred, Y_train))
    train_loss.append(_cross_entropy_loss(y_train_pred, Y_train,) / train_size)

    y_dev_pred = _f(X_dev, w, b)
    Y_dev_pred = np.round(y_dev_pred)
    dev_acc.append(_accuracy(Y_dev_pred, Y_dev))
    dev_loss.append(_cross_entropy_loss(y_dev_pred, Y_dev) / dev_size)
    print('Training loss: {}'.format(train_loss[-1]))
    print('Development loss: {}'.format(dev_loss[-1]))
    print('Training accuracy: {}'.format(train_acc[-1]))
    print('Development accuracy: {}'.format(dev_acc[-1]))

"""------------------------画图----------------------------"""
# 损失曲线(Loss curve)
plt.plot(train_loss)
plt.plot(dev_loss)
plt.title('Loss')
plt.legend(['train', 'dev'])
plt.savefig('loss.png')
plt.show()

# 精度曲线(Accuracy curve)
plt.plot(train_acc)
plt.plot(dev_acc)
plt.title('Accuracy')
plt.legend(['train', 'dev'])
plt.savefig('acc.png')
plt.show()


"""------------------------------保存-------------------------------"""
# 保存预测测试集标签 (Predict testing labels)
predictions = _predict(X_test, w, b)
with open(output_fpath.format('logistic'), 'w') as f:
    f.write('id,label
')
    for i, label in  enumerate(predictions):
        f.write('{},{}
'.format(i, label))

# 对测试集归一化处理  打印出最重要的权重w
ind = np.argsort(np.abs(w))[::-1]  #从后往前
with open(X_test_fpath) as f:
    content = f.readline().strip('
').split(',')
features = np.array(content)
for i in ind[0:10]:
    print('*************')
    print(features[i], w[i])