项目2-年收入判断

友情提示

同学们可以前往课程作业区先行动手尝试！！！

项目描述

二元分类是机器学习中最基础的问题之一，在这份教学中，你将学会如何实作一个线性二元分类器，来根据人们的个人资料，判断其年收入是否高于 50,000 美元。我们将以两种方法: logistic regression 与 generative model，来达成以上目的，你可以尝试了解、分析两者的设计理念及差别。
实现二分类任务：

• 个人收入是否超过50000元？

数据集介绍

这个资料集是由UCI Machine Learning Repository 的Census-Income (KDD) Data Set 经过一些处理而得来。为了方便训练，我们移除了一些不必要的资讯，并且稍微平衡了正负两种标记的比例。事实上在训练过程中，只有 X_train、Y_train 和 X_test 这三个经过处理的档案会被使用到，train.csv 和 test.csv 这两个原始资料档则可以提供你一些额外的资讯。

• 已经去除不必要的属性。
• 已经平衡正标和负标数据之间的比例。

特征格式

1. train.csv，test_no_label.csv。

• 基于文本的原始数据
• 去掉不必要的属性，平衡正负比例。

2. X_train, Y_train, X_test(测试)

• train.csv中的离散特征=>在X_train中onehot编码(学历、武功状态…)
• train.csv中的连续特征 => 在X_train中保持不变(年龄、资本损失…)。
• X_train, X_test : 每一行包含一个510-dim的特征，代表一个样本。
• Y_train: label = 0 表示 “<=50K” 、 label = 1 表示 " >50K " 。

项目要求

1. 请动手编写 gradient descent 实现 logistic regression
2. 请动手实现概率生成模型。
3. 单个代码块运行时长应低于五分钟。
4. 禁止使用任何开源的代码(例如，你在GitHub上找到的决策树的实现)。

数据准备

项目数据保存在：work/data/ 目录下。

环境配置/安装

无

Logistic回归

首先我们会做 Logistic回归

数据准备

下载资料，并且对每个属性做正规化，处理过后再将其切分为训练集与验证集。

import numpy as np
import pandas as pd

X_train_fpath = 'work/data/X_train'
with open(X_train_fpath) as f:
    X_train = np.array([line.strip('\n').split(',')[1:] for line in f])
    print(X_train)
    X_train = pd.DataFrame(X_train[1:],index=None,columns = X_train[0])
    print(X_train.head())
print(X_train.shape)
1
2
3
4
5
6
7
8
9
10

[['age' ' Private' ' Self-employed-incorporated' ...
  'weeks worked in year' ' 94' ' 95']
 ['33' '1' '0' ... ' 52' '0' '1']
 ['63' '1' '0' ... ' 52' '0' '1']
 ...
 ['16' '0' '0' ... ' 8' '1' '0']
 ['48' '1' '0' ... ' 52' '0' '1']
 ['48' '0' '0' ... ' 0' '0' '1']]
 
 
  age  Private  Self-employed-incorporated  State government  ...
0  33        1                           0                 0   
1  63        1                           0                 0   
2  71        0                           0                 0   
3  43        0                           0                 0   
4  57        0                           0                 0   
[5 rows x 510 columns]


(54256, 510)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

import numpy as np

np.random.seed(0)
X_train_fpath = 'work/data/X_train'
Y_train_fpath = 'work/data/Y_train'
X_test_fpath = 'work/data/X_test'
output_fpath = 'work/output_{}.csv'

# 将CSV文件解析为NumPy数组
with open(X_train_fpath) as f:
    next(f)
    X_train = np.array([line.strip('\n').split(',')[1:] for line in f], dtype = float)
with open(Y_train_fpath) as f:
    next(f)
    Y_train = np.array([line.strip('\n').split(',')[1] for line in f], dtype = float)
with open(X_test_fpath) as f:
    next(f)
    X_test = np.array([line.strip('\n').split(',')[1:] for line in f], dtype = float)

def _normalize(X, train = True, specified_column = None, X_mean = None, X_std = None):
    # 此函数用于归一化X的特定列。
    # 在处理测试数据时，将重复使用训练数据的平均值和标准差。
    # 参数：
    #     X：待处理数据。
    #     Train：处理训练数据时为‘True’，处理测试数据时为‘False。
    #     SPECIAL_COLUMN：将被规格化的列的索引。如果为‘None’，则为所有列。
    # 归一化。
    #     X_Mean：训练数据的平均值，当Train=‘False’时使用。
    #     X_STD：训练数据的标准差，当Train=‘False’时使用。
    # 输出：
    #     X：归一化数据。
    #     X_Mean：训练数据的计算平均值。
    #     X_STD：训练数据的计算标准差

    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):
    # 该功能将数据分解为训练集和验证集,dev_ratio代表验证集占的比例
    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 validation set: {}'.format(dev_size))
print('Size of testing set: {}'.format(test_size))
print('Dimension of data: {}'.format(data_dim))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

Size of training set: 48830
Size of validation set: 5426
Size of testing set: 27622
Dimension of data: 510
1
2
3
4

一些有用的函数

这几个函数可能会在训练回圈中被重复使用到。

def _shuffle(X, Y):
    # 这个函数将两个相等长度的列表/数组X和Y一起打乱
    randomize = np.arange(len(X))
    np.random.shuffle(randomize)
    return (X[randomize], Y[randomize])

def _sigmoid(z):
    # Sigmoid函数可以用来计算概率。
    # 为了避免溢出，设置了最小输出值1e-8，最大输出值1 - (1e-8)。
    return np.clip(1 / (1.0 + np.exp(-z)), 1e-8, 1 - (1e-8))

def _f(X, w, b):
    # 这是逻辑回归函数，参数为w和b
    # 参数:
    #   X: input data, shape = [batch_size, data_dimension]
    #   w: 权重向量，形状= [data_dimension，]
    #   b: 偏置值，标量
    # 输出:
    #   np.matmul返回两个数组的矩阵乘积：z=x*w+b
    #   X：每一行被正标记的预测概率，shape = [batch_size，]
    return _sigmoid(np.matmul(X, w) + b)

def _predict(X, w, b):
    # 这个函数返回X每一行的真值预测
    # 通过对逻辑回归函数的结果进行四舍五入。
    return np.round(_f(X, w, b)).astype(np.int)
    
def _accuracy(Y_pred, Y_label):
    # 这个函数计算预测精度
    # np.abs求绝对值
    acc = 1 - np.mean(np.abs(Y_pred - Y_label))
    return acc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

梯度与损失

def _cross_entropy_loss(y_pred, Y_label):
    # 这个函数计算交叉熵。
    # 输入:
    #     y_pred: 概率预测，浮点向量
    #     Y_label: 真实标签，bool向量
    # 输出:
    #     np.dot作用是向量点积和矩阵乘法
    #     log下什么都不写默认是求自然对数
    #     cross_entropy：交叉熵，标量，cross_entropy = y真*lny预 - (1-y真)*ln(1-y预)
    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和偏置值b的交叉熵损失梯度
    # np.sum中，当axis为0时,是压缩行,即将每一列的元素相加,将矩阵压缩为一行
    #           当axis为1时,是压缩列,即将每一行的元素相加,将矩阵压缩为一列
    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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

模型训练

一切准备就绪，开始训练吧!

我们使用小批次梯度下降法来训练。训练资料被分为许多小批次，针对每一个小批次，我们分别计算其梯度以及损失，并根据该批次来更新模型的参数。当一次回圈完成，也就是整个训练集的所有小批次都被使用过一次以后，我们将所有训练资料打散并且重新分成新的小批次，进行下一个回圈，直到事先设定的回圈数量达成为止。

#  权值和偏差的零初始化
# np.zeros((data_dim,))：得到一行全为0的列表，个数为：data_dim
w = np.zeros((data_dim,)) 
b = np.zeros((1,))

# 训练的一些参数
max_iter = 10
batch_size = 8
learning_rate = 0.2

# 保持每次迭代的损耗和精度，进行绘图
train_loss = []
dev_loss = []
train_acc = []
dev_acc = []

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

# 迭代训练
for epoch in range(max_iter):
    # 在每轮的迭代中随机打乱训练集X_train和验证集Y_train
    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 - learning_rate/np.sqrt(step) * w_grad
        b = b - learning_rate/np.sqrt(step) * b_grad

        step = step + 1
            
    # 计算训练集和验证集的损耗和精度
    y_train_pred = _f(X_train, w, b)
    Y_train_pred = np.round(y_train_pred)
    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('validation loss: {}'.format(dev_loss[-1]))
print('Training accuracy: {}'.format(train_acc[-1]))
print('validation accuracy: {}'.format(dev_acc[-1]))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

Training loss: 0.271355435246406
validation loss: 0.28963596750262866
Training accuracy: 0.8836166291214418
validation accuracy: 0.8733873940287504
1
2
3
4

绘制损失和精度曲线

%matplotlib inline
import matplotlib.pyplot as plt

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

# 精度曲线
plt.plot(train_acc)
plt.plot(dev_acc)
plt.title('Accuracy')
plt.legend(['train', 'dev'])
plt.savefig('acc.png')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

预测测试标签

预测测试集的资料标签并且存在 output_logistic.csv 中。

# 预测测试集标签
predictions = _predict(X_test, w, b)
with open(output_fpath.format('logistic'), 'w') as f:
    f.write('id,label\n')
    for i, label in enumerate(predictions):
        f.write('{},{}\n'.format(i, label))

# 打印出最重要的10个权重w
# np.argsort返回的是元素值从小到大排序后的索引值的数组
# [::-1]将数组倒序
ind = np.argsort(np.abs(w))[::-1]
with open(X_test_fpath) as f:
    content = f.readline().strip('\n').split(',')
features = np.array(content)
for i in ind[0:10]:
    print(features[i], w[i])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

 Not in universe -4.031960278019252
 Spouse of householder -1.6254039587051405
 Other Rel <18 never married RP of subfamily -1.4195759775765409
 Child 18+ ever marr Not in a subfamily -1.2958572076664745
 Unemployed full-time 1.1712558285885908
 Other Rel <18 ever marr RP of subfamily -1.1677918072962366
 Italy -1.0934581438006177
 Vietnam -1.0630365633146412
num persons worked for employer 0.9389922773566517
 1 0.822661492211719
1
2
3
4
5
6
7
8
9
10

多变量生成模型

接着我们将实作基于 generative model 的二元分类器。

数据准备

训练集与测试集的处理方法跟 logistic regression 一模一样，然而因为 generative model 有可解析的最佳解，因此不必使用到验证集。

# 将CSV文件解析为NumPy数组
with open(X_train_fpath) as f:
    next(f)
    X_train = np.array([line.strip('\n').split(',')[1:] for line in f], dtype = float)
with open(Y_train_fpath) as f:
    next(f)
    Y_train = np.array([line.strip('\n').split(',')[1] for line in f], dtype = float)
with open(X_test_fpath) as f:
    next(f)
    X_test = np.array([line.strip('\n').split(',')[1:] for line in f], dtype = float)

# 归一化
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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

平均值和协方差

在 generative model 中，我们需要分别计算两个类别内的资料平均与共变异。

# 计算类内平均值
X_train_0 = np.array([x for x, y in zip(X_train, Y_train) if y == 0])
X_train_1 = np.array([x for x, y in zip(X_train, Y_train) if y == 1])

mean_0 = np.mean(X_train_0, axis = 0)
mean_1 = np.mean(X_train_1, axis = 0)  

# 计算类内协方差
cov_0 = np.zeros((data_dim, data_dim))
cov_1 = np.zeros((data_dim, data_dim))

for x in X_train_0:
    cov_0 += np.dot(np.transpose([x - mean_0]), [x - mean_0]) / X_train_0.shape[0]
for x in X_train_1:
    cov_1 += np.dot(np.transpose([x - mean_1]), [x - mean_1]) / X_train_1.shape[0]

# 共享协方差是单个类内协方差的加权平均值。
cov = (cov_0 * X_train_0.shape[0] + cov_1 * X_train_1.shape[0]) / (X_train_0.shape[0] + X_train_1.shape[0])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

计算权重和偏差

权重矩阵与偏差向量可以直接被计算出来。

# 计算协方差矩阵的逆矩阵。
# 由于协方差矩阵可能近乎奇异，np.linalg.inv()可能会给出较大的数值误差。
# 通过SVD分解，可以高效准确地得到矩阵的逆。
u, s, v = np.linalg.svd(cov, full_matrices=False)
inv = np.matmul(v.T * 1 / s, u.T)

# 直接计算权重w和偏差b
w = np.dot(inv, mean_0 - mean_1)
b =  (-0.5) * np.dot(mean_0, np.dot(inv, mean_0)) + 0.5 * np.dot(mean_1, np.dot(inv, mean_1))\
    + np.log(float(X_train_0.shape[0]) / X_train_1.shape[0]) 

# 计算训练集的准确度
Y_train_pred = 1 - _predict(X_train, w, b)
print('Training accuracy: {}'.format(_accuracy(Y_train_pred, Y_train)))
1
2
3
4
5
6
7
8
9
10
11
12
13
14

Training accuracy: 0.8671114715423179
1

预测测试标签

预测测试集的资料标签并且存在 output_generative.csv 中。

# 预测测试集标签
predictions = 1 - _predict(X_test, w, b)
with open(output_fpath.format('generative'), 'w') as f:
    f.write('id,label\n')
    for i, label in  enumerate(predictions):
        f.write('{},{}\n'.format(i, label))

# 打印出最重要的10个权重w
ind = np.argsort(np.abs(w))[::-1]
with open(X_test_fpath) as f:
    content = f.readline().strip('\n').split(',')
features = np.array(content)
for i in ind[0:10]:
    print(features[i], w[i])
1
2
3
4
5
6
7
8
9
10
11
12
13
14

 Retail trade 7.67333984375
 Midwest -6.3125
 34 -5.835205078125
 37 -5.489013671875
 Child <18 ever marr not in subfamily -5.4759521484375
 Other service -5.00390625
 Different county same state 4.66796875
 33 -3.91015625
 Private household services 3.8623046875
 32 -3.51953125
1
2
3
4
5
6
7
8
9
10