【机器学习】李宏毅-判断年收入

李宏毅-判断年收入

1 实验目的

本次作业的数据是加州大学尔湾分校机器学习作业中下载得到，使用Classification中的生成模型generative model以及logistic regression解决二分类问题。

根据已有数据，判断该人年收入是否大于5万美元，最终得到预测结果输出到结果文件中比较准确率。

2 实验要求

• 不可以使用 tensorflow 或者 pytorch 库

• 请用概率生成模型generative model(分类第一节课)解决本问题

• 请用logistic regression(分类第二节课)解决本问题

• 上传格式为csv， 第一行必须为id，label，第二行开始为预测结果，每行分别为id以及预测的label，请以逗号分隔

3 实验环境

3.1 硬件环境

笔记本电脑、Intel Core i5

3.2 软件环境

windows10操作系统、Python 3.8 64-bit、Visual Studio Code

4 数据处理

4.1 数据文件描述

关于train.csv, test.csv：

age, workclass, fnlwgt (总人数), education, education num, marital-status, occupation, relationship, race, sex, capital-gain, capital-loss, hours-per-week, native-country, make over 50K a year or not

关于X_train, Y_train, X_test ：

1. discrete features in train.csv => one-hot encoding in X_train (work_class,education…)
2. continuous features in train.csv => remain the same in X_train (age,capital_gain…)
3. X_train, X_test : each row contains one 106-dim feature represents a sample
4. Y_train: label = 0 means “<= 50K”label = 1 means “ >50K ”

4.2 读取数据

在data_manager()类中定义了读取数据的函数read()：

#读取数据
def read(self,name,path):
    with open(path,newline = '') as csvfile:
        #注意去掉首行读取
        rows = np.array(list(csv.reader(csvfile))[1:] ,dtype = float) 

通过调用该函数读取X_train, Y_train, X_test 三个数据文件：

#读取文件
dm.read('X_train','data/X_train')
dm.read('Y_train','data/Y_train')
dm.read('X_test','data/X_test')

4.3 z-score特征标准化

z-score通过（x-μ）/σ将两组或多组数据转化为无单位的Z-Score分值，使得数据标准统一化，提高了数据可比性，削弱了数据解释性：

        if name == 'X_train':
                #求均值并转为1行,div=1*106
                self.mean = np.mean(rows,axis = 0).reshape(1,-1)
                #求方差并转为1行,div=1*106
                self.std = np.std(rows,axis = 0).reshape(1,-1)
                self.theta = np.ones((rows.shape[1] + 1,1),dtype = float)
                #对每一行数据进行标准化
                for i in range(rows.shape[0]):
                    rows[i,:] = (rows[i,:]-self.mean)/self.std
            elif name == 'X_test': 
                for i in range(rows.shape[0]):
                    rows[i,:] = (rows[i,:]-self.mean)/self.std
            #保存数据处理结果到类内字典索引里
            self.data[name] = rows

4.4 划分训练集和验证集

data_manager()类内函数partition()负责训练集和验证集的划分，其中num代表设定验证集里有多少笔验证数据，在此暂时设为3000：

#分割数据集为训练集和验证集
def partition(self):
    num=3000
    self.data['X_verify']=self.data['X_train'][:num]
    self.data['X_train']=self.data['X_train'][num:]
    self.data['Y_verify']=self.data['Y_train'][:num]
    self.data['Y_train']=self.data['Y_train'][num:]

由此可见训练集在self.data字典里的索引为X_train、Y_train，验证集在self.data字典里的索引为X_verify、Y_verify，测试集在self.data字典里的索引为X_test。

5 概率生成模型generative model

概率生成模型，简称生成模型（Generative Model），是概率统计和机器学习中的一类重要模型，指一系列用于随机生成可观测数据的模型。生成模型可以和贝叶斯概率公式进行结合，用于分类问题。

5.1 对训练数据二分类

由于概率生成模型需要计算样本选择的概率，所以需要先对训练数据进行二分类：

#记录分类的元素索引
class_0_id = []
class_1_id = []
#通过Y_train分类索引
for i in range(self.data['Y_train'].shape[0]):
    if self.data['Y_train'][i][0] == 0:
        class_0_id.append(i)
    else:
        class_1_id.append(i)
#按照索引去X_train找到对应元素
class_0 = self.data['X_train'][class_0_id]
class_1 = self.data['X_train'][class_1_id] 

5.2 计算均值和协方差

参考均值的计算公式：

在代码中实现：

#分别求出两类均值
mean_0 = np.mean(class_0,axis = 0)
mean_1 = np.mean(class_1,axis = 0) 

不同的类别共享协方差矩阵，因为协方差的规模是和特征数量的平方成正比的，如果特征数量很大，协方差的尺寸也大，由于不同的高斯分布有不同的协方差矩阵，这样模型的参数就会很多，就会过拟合。计算协方差的公式为：

在代码中实现：

n = class_0.shape[1]
#两类的协方差,div=106*106
cov_0 = np.zeros((n,n))
cov_1 = np.zeros((n,n))
#计算类1的协方差
for i in range(class_0.shape[0]):
    #需要进行矩阵转置,div从1*106转置为106*1
    cov_0 += np.dot(np.transpose([class_0[i]-mean_0]),[class_0[i]-mean_0])
#计算类2的协方差
    for i in range(class_1.shape[0]):
    #需要进行矩阵转置,div从1*106转置为106*1
    cov_1 += np.dot((class_1[i]-mean_1).transpose(),class_1[i]-mean_1)
#类中的行数
n_0 = class_0.shape[0]
n_1 = class_1.shape[0]
#两个类共享一个协方差
cov = (cov_0+cov_1)/(n_0+n_1)

5.3 计算参数w和b

参考计算参数w和b的公式：

从公式中可以看出，我们不光需要刚才计算出的均值和协方差，还需要协方差的逆，代码实现如下：

#协方差的逆
cov_inv = inv(cov)
#根据公式求w和b
self.w = np.dot((mean_0-mean_1).transpose(),cov_inv).transpose()
self.b = (-0.5)*np.dot(np.dot(mean_0.transpose(),cov_inv),mean_0)+0.5*np.dot(np.dot(mean_1.transpose(),cov_inv),mean_1)+np.log(float(n_0)/n_1)

5.4 后验概率函数sigmoid

后验概率参考公式：

代码实现如下，注意这里的np.exp()会存在溢出风险，所以使用了拆分运算（在逻辑回归模型中改进了更好的办法）：

#后验概率函数sigmoid
def func(self,x):
    arr = np.empty([x.shape[0],1],dtype=float)
    ##拆分运算，防止exp溢出
    for i in range(x.shape[0]):
        z = x[i,:].dot(self.w) + self.b
        z *= (-1)
        arr[i][0] = 1 / (1 + np.exp(z))
    #避免由于数值太小，通过numpy的clip函数，将其转换为了界于1e-8和1-(1e-8)的数字
    return np.clip(arr, 1e-8, 1-(1e-8))

5.5 分类模型输出

通过5.4中后验概率的计算结果，我们需要对模型进行输出，后验概率值＞0.5则输出1，后验概率值＜0.5则输出0，定义类内函数predict()来完成这一过程：

#函数输出
def predict(self,x):
    ans = np.ones([x.shape[0],1],dtype=int)
    for i in range(x.shape[0]):
        if x[i] > 0.5:
            ans[i] = 0
    return ans

5.6 训练集和验证集上的准确率

计算出w和b之后就可以直接计算出训练的生成模型在训练集上的准确率：

#用求得的w和b计算sigmoid()函数结果
result = self.func(self.data['X_train'])
#经过0-1函数处理得到预测值
answer = self.predict(result)
#计算训练集上的准确率
accuracy=1-np.mean(np.abs(self.data['Y_train']-answer))
print("accuracy on train set = ",accuracy)

另外我们用类内函数verify()计算在验证集上的准确率，验证集来自于4.4环节划分得到：

#验证集验证
def verify(self):
    #用求得的w和b计算sigmoid()函数结果
    result = self.func(self.data['X_verify'])
    #经过0-1函数处理得到预测值
    answer = self.predict(result)
    #计算验证集上的准确率
    accuracy=1-np.mean(np.abs(self.data['Y_verify']-answer))
    print("accuracy on verify set = ",accuracy)

运行程序得到训练的生成模型在训练集和验证集上各自的准确率：

5.7 模型预测结果

类内函数wirte_file()可以把模型的预测结果导出到格式为csv的文件里， 第一行为id，label，第二行开始为预测结果，每行分别为id以及预测的label，以逗号分隔：

#预测结果并输出到目标文件中
def write_file(self,path):
    result = self.func(self.data['X_test'])
    answer = self.predict(result)
    with open(path, 'w', newline='') as csvfile:
        writer = csv.writer(csvfile) 
        writer.writerow(['id','label']) 
        for i in range(answer.shape[0]):
            writer.writerow([i+1,answer[i][0]])

运行程序得到预测结果，详见结果文件generative_model_output.csv。

6 逻辑回归模型logistic regression

逻辑回归假设数据服从伯努利分布（即数据的标签为0或者1），通过极大化似然函数的方法，运用梯度下降来求解参数，来达到将数据二分类的目的。

使用的基础模型为：

6.1 参数设定

在梯度下降的过程中，我们需要诸多参数，在此预先定义好：

#学习率
learning_rate = 1
#迭代次数
iter_time = 1000
#初始化权重w和偏移b
self.w = np.zeros([len(self.data['X_train'][0]),1])
self.b = 0
#初始化adagrad参数
w_adagrad = np.ones([len(self.data['X_train'][0]),1])
b_adagrad = 1

6.2 梯度下降迭代过程

迭代过程类似上一次线性回归的实验，但是使用的偏微分公式发生了变化：

因此参数更新的公式变为：

另外为了更好地训练模型，使用Adagrad方法更新学习率：

代码实现如下：

#进度条对象
p=progressbar.ProgressBar()
#梯度下降法迭代
for i in p(range(iter_time)):
     result = self.func(self.data['X_train'])
     #求w梯度
     w_gradient = -np.dot(self.data['X_train'].transpose(),self.data['Y_train']-result)
     w_adagrad += (w_gradient/self.data['X_train'].shape[0]) ** 2
     #求b梯度
     b_gradient = -np.sum(self.data['Y_train']-result)
     b_adagrad += (b_gradient/self.data['X_train'].shape[0]) ** 2
     #更新w
     self.w = self.w - learning_rate * w_gradient / np.sqrt(w_adagrad)
     #更新b
     self.b = self.b - learning_rate * b_gradient / np.sqrt(b_adagrad)

6.3 后验概率函数sigmoid

在概率生成模型中，由于运算规模较小，并且使用了拆分运算，所以没有出现exp的overflow问题，在逻辑回归模型的最开始的调试过程中我们发现了这个问题，引用scipy.special模块中的expit函数解决了溢出的问题：

from scipy.special import expit

#sigmoid函数
def func(self,x):
    arr = np.empty([x.shape[0],1],dtype=float)
    #拆分运算，防止exp溢出
    for i in range(x.shape[0]):
        z = x[i,:].dot(self.w) + self.b
        z *= (-1)
        arr[i][0] = 1 / (1 + expit(z))
    #避免由于数值太小，通过numpy的clip函数，将其转换为了界于1e-8和1-(1e-8)的数字
    return np.clip(arr, 1e-8, 1-(1e-8))

6.4 训练集和验证集上的准确率

计算出w和b之后就可以直接计算出训练的逻辑回归模型在训练集上的准确率：

#用求得的w和b计算sigmoid()函数结果
result = self.func(self.data['X_train'])
#经过0-1函数处理得到预测值
answer = self.predict(result)
#计算训练集上的准确率
accuracy=1-np.mean(np.abs(self.data['Y_train']-answer))
print("accuracy on train set = ",accuracy)

另外我们用类内函数verify()计算在验证集上的准确率，验证集来自于4.4环节划分得到：

#验证集验证
def verify(self):
    #用求得的w和b计算sigmoid()函数结果
    result = self.func(self.data['X_verify'])
    #经过0-1函数处理得到预测值
    answer = self.predict(result)
    #计算验证集上的准确率
    accuracy=1-np.mean(np.abs(self.data['Y_verify']-answer))
    print("accuracy on verify set = ",accuracy)

运行程序得到训练的逻辑回归模型在训练集和验证集上各自的准确率：

6.5 正则化regularization

正则化的思想是在损失函数中加入限制项，公式如下：

定义正则化系数λ：

regularization_lambda=1

修改梯度下降迭代过程，考虑正则项的微分：

#梯度下降法迭代过程
for i in p(range(iter_time)):
     result = self.func(self.data['X_train'])
     #求w梯度
     w_gradient = -np.dot(self.data['X_train'].transpose(),self.data['Y_train']-result)
     w_gradient += regularization_lambda*w_gradient
     w_adagrad += (w_gradient/self.data['X_train'].shape[0]) ** 2
     #求b梯度
     b_gradient = -np.sum(self.data['Y_train']-result)
     b_adagrad += (b_gradient/self.data['X_train'].shape[0]) ** 2
     #更新w
     self.w = self.w - learning_rate * w_gradient / np.sqrt(w_adagrad)
     #更新b
     self.b = self.b - learning_rate * b_gradient / np.sqrt(b_adagrad)

运行程序得到新的结果（对于正则化的分析评估将在7.3展开）：

6.6 模型预测结果

类内函数wirte_file()可以把模型的预测结果导出到格式为csv的文件里， 第一行为id，label，第二行开始为预测结果，每行分别为id以及预测的label，以逗号分隔。

代码相较于生成模型作出了适当修改：

#预测结果并输出到目标文件中
def write_file(self,path):
    result = self.func(self.data['X_test'])
    answer = np.around(result)
    with open(path, 'w', newline='') as csvfile:
        writer = csv.writer(csvfile) 
        writer.writerow(['id','label']) 
        for i in range(answer.shape[0]):
            writer.writerow([i+1,int(answer[i][0])])

运行程序得到预测结果，详见结果文件logistic_regression_output.csv。

7 分析评估

7.1 不同分类模型的影响

用逻辑回归的方法生成的模型叫做判别模型，用高斯公式描述后验概率称为生成模型。逻辑回归和生成模型用的一样的模型。前者直接在训练数据的特征空间进行学习，直接找出w，b，后者会考虑数据的生成情况，即每个类别下的数据分布规律，用概率分布算出w和b。

虽然是同一个函数，但是用两种方法得到的w和b结果是不同的。因为在生成式模型中，我们对概率分布是有假设的，假设是高斯分布或伯努利分布。

下面是不同模型在不同数据集上的的准确率的汇总：

可以看到逻辑回归模型相较于概率生成模型的拟合效果更好，在训练集和验证集上的表现也更出色。

但是逻辑回归模型的表现依赖于其一系列参数的选取，例如学习率大小和正则化系数等等，其次样本量也对测试结果有一定影响，因此并不能保证说逻辑回归模型优于概率生成模型。

7.2 特征标准化feature scaling对模型准确率的影响

z-score通过（x-μ）/σ将两组或多组数据转化为无单位的Z-Score分值，使得数据标准统一化，提高了数据可比性，削弱了数据解释性。

7.3 正则化regularization对模型准确率的影响

加入正则化的目的就是让模型的参数变小，使模型对数据敏感度下降，所以正则项的加入就是降低模型变化率的过程，变化率降低直观表现为模型更加平滑，使模型更加集中且增大了与真实模型之间的距离，也就是增大了偏移量bias，从而增强了高阶过拟合模型的泛化能力。

下面是引入正则项前后逻辑回归模型在不同数据集上的准确率的汇总：

可以看到，正则化对模型在数据集上的准确率都有略微提升，但是提升的效果并不是非常明显。

究其原因依然是正则化对高次项的过拟合现象有较好的改良效果，在本次实验中的逻辑回归模型并没有使用高次项作为参数，因此使用正则化之后对结果的影响并不是特别明显。

7.4 猜测对结果影响最大的特征

在逻辑回归模型得到的结果中，找到权重w里绝对值最大的一项，该w对应的特征可能就是对结果影响最大的特征。定义类内函数find_max()来完成这一过程：

#猜测影响最大的特征
def find_max(self):
    max=np.max(dm.w)
    for i in range(len(self.w)):   
        if self.w[i]==max:
            print("The most important feature is No.",i)
            break

得到函数输出：

可以看到第三个特征对结果的影响最大，对应的特征是capital_gain。

8 References

[1] python.org

[2] kaggle.com

[3] 百度百科

[4] CSDN专业开发者技术社区