数据挖掘专栏二-Matlab-消费者人群画像信用智能评分_用户画像 matlab代码

文章对应的代码和实验报告在百度云和github上都放了一份。百度云地址上个专栏有提。github仓库地址:https://github.com/gamblerInCoding/Portrait

项目名称

    消费者人群画像-信用智能评分

数据集说明

数据集名称

    消费者人群画像-信用智能评分数据集

数据集来源

    https://www.datafountain.cn/competitions/337/datasets

数据集介绍

    数据的原始来源是中国移动福建公司提供2018年某月份的样本数据（脱敏），包括客户的各类通信支出、欠费情况、出行情况、消费场所、社交、个人兴趣等丰富的多维度数据。该数据集包含了2个csv文件，train_dataset.csv文件和test_dataset.csv文件。其中train_dataset.csv文件是用来给我们提供训练机器学习算法模型的原始数据文件，包含了某月中各个用户的上网行为信息及其信用评价的分数。数据集大小为5.43MB。整个数据集有50000条记录，用户编码，用户实名制是否通过核实，用户年龄，是否大学生落户，是否黑名单落户，是否4G不健康客户，用户网龄，用户最近一次缴费距今时长，缴费用户最近一次缴费金额，用户近6个月平均消费话费，用户账单当月总费用，用户当月账户余额，缴费用户当前是否欠费缴费，用户话费敏感度等30个特征，其中信用分作为反应变量，其他29个特征作为解释变量，在解释特征中有一个用来标识不同用户上网行为记录的编码。每条记录都由这30个字段组成，代表着一名用户在该月的行为特征及其获得的信用分数。

数据集字段说明

   数据集的字段说明如下表所示。

实验环境

    MATLAB 版本：9.10.0.1602886 (R2021a)
    MATLAB 许可证编号: 968398
    操作系统: Microsoft Windows 10 家庭中文版 Version 10.0 (Build 19043)
    Java 版本: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot™ 64-Bit Server VM mixed mode

实验步骤

实验过程

数据集探索

数据集获取和读入

    我选择了importdata这个函数来读入原始csv文件中的数据。

datas = importdata("train_dataset.csv");
1

数据集基本信息

    鉴于Matlab里面并没有直接提供针对数据集描述性的函数，所以我编写了2个函数来针对数据集进行描述性统计。一个函数是dts，统计了对应列的所有数据的均值，方差，标准差，极差，变异系数，偏度，峰度。另外一个函数是fws，统计了对应列的所有数据的中位数，上下四分位数，四分位极差。这是针对列进行数据统计的，所以传进来的x是datas中的一列。

function dts(x)
a = x(:);
nans = isnan(a);
ind = find (nans); %nan是0/0.
a(ind)=[];
xbar= mean(a);
disp(['均值是：',num2str(xbar)]);
s2 = var(a);
disp(['方差是：',num2str(s2)]);
s = std(a);
disp(['标准差是：',num2str(s)]);%数据里必须是元素的类型一样，所以要有num2str()函数转一下。
R = range(a);
disp(['极差是：',num2str(R)]);
cv = 100*s./xbar;%它是一个相对的数且没有量纲，所以更具有说明性。
disp(['变异系数是：',num2str(cv)]);
g1 = skewness(a,0);
disp(['偏度：',num2str(g1)]);
g2=kurtosis(a,0);
disp(['峰度',num2str(g2)]);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

function fws(x)
a = x(:);
a(isnan(a))=[];
ss5 = prctile(a,50);
disp(['中位数是：',num2str(ss5)]);
ss25 = prctile(a,25);
disp(['下四分位数是：',num2str(ss25)]);
ss75 = prctile(a,75);
disp(['上四分位数是：',num2str(ss75)]);
RS = ss75-ss25;
disp(['四分位极差：',num2str(RS)]);
end
1
2
3
4
5
6
7
8
9
10
11
12

数据预处理

重复值处理

    由于我把由于标识用户的唯一标记用户编码字段放在了name里，所以在data里面的只是纯数据，我们要去查看data中是否存在重复值，如果数据集中存在了大量的重复值会影响机器学习算法模型的计算和求解过程。在Matlab中我们可以使用unique这个函数，这个函数能直接得到原始数据去掉完全相同的行后剩下的数据。

data=unique(data,'rows');
1

缺失值处理

    Matlab里面提供了isnan这个函数来判断数据的缺失情况。对于缺失值，我们可以使用删除法、替换法、填补法。我们可以用如下命令判断缺失值数量

sum(isnan(data))
1

异常值处理

    异常值指的是数据极大或者极小，和正常情况下的数据分布情况偏离了太多，针对这种数据我们可以先使用箱型图对数据分布进行观察。Matlab中提供了boxplot命令让我们绘制箱型图。
    由于原始数据有29个不同的特征，绘制在一张图上所以我绘制了3个不同的figure，每个figure上分别绘制10个特征。绘制箱型图的代码如下所示：

disp("异常值")
figure(1);boxplot(data(:,1:10))
title('变量1-10箱形图','fontsize',12)
figure(2);boxplot(data(:,11:20))
title('变量11-20箱形图','fontsize',12)
figure(3);boxplot(data(:,21:29))
title('变量21-29箱形图','fontsize',12)
1
2
3
4
5
6
7

    绘制出的箱型图图像如下所示。
                                                                                                                            常见的针对异常值进行处理的方法都需要数据满足一定的分布，即正态分布或者类似正态分布。而有些方法需要使用作图工具箱而，所以我在接下来处理异常值的时候使用了肖维勒方法来剔除异常值。其基本思想是规定一个置信水平,确定一个置信限度,凡是超过该限度的误差, 就认为它是异常值,从而予以剔除。

% 使用肖维勒方法（等置信概率）剔除异常值
[m n] = size(data);
Y = [];            
w = 1 + 0.4*log(m);    % 肖维勒系数（近似计算公式）
for i = 1:n
   x = data(:,i);    
   YiChang = abs(x-mean(x)) > w*std(x);
   Y(:,i) = YiChang;
end
[u v] = find(Y() == 1);   % 找出异常值所在的行与列
ls = size(u,1);
uu = unique(u);    % 剔除重复的行数
now = size(uu,1);
disp("剔除异常值的数量:");
size(uu,1)
data(uu,:) = [ ];   %令异常值所在行为空,即剔除异常值
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

    查看统计结果，发现删除了3355个异常值。

数据可视化探索

分类变量分析-用户话费敏感度

    用户话费敏感度是一个比较重要的特征，下面对它进行探索性分析。我们可以先使用Matlab提供的tabulate函数来统计离散变量的频数分布，这个函数能返回每个函数值和函数值对应的频数值。

target=data(:,29);
phoneprice=data(:,13);
result = tabulate(phoneprice(:));
lab = num2str(result(:,1));
figure(4)
subplot(1,2,1)
bar(result(:,1),result(:,2))
title('用户话费敏感度分布条形图','fontsize',12)
subplot(1,2,2)
pie(result(:,2))
title('用户话费敏感度分布饼状图','fontsize',12)
legend(lab);
1
2
3
4
5
6
7
8
9
10
11
12

分类变量分析-用户最近一次缴费距今时长

    用户最近一次缴费距今时长越短，说明用户缴费的意愿越高，该用户的收入情况肯定也会越好，其信用分的水平也会越高。我们在这个变量的可视化分析中，为了更进一步看出用户最近一次缴费距今时长这个变量和信用分之间的关系，我统计了用户最近一次缴费距今时长相同的用户的信用分的平均值。

disp("用户最近一次缴费距今时长（月）")
ms = [0,0];
timetonow=data(:,7);
for i=1:46645
    if timetonow(i,1)==0
        ms(1,1)=ms(1,1)+target(i,1);
    end
    if timetonow(i,1)==1
        ms(1,2)=ms(1,2)+target(i,1);
    end
end
result = tabulate(timetonow(:));
for i=1:2
    ms(1,i)=ms(1,i)/result(i,2);
end
lab = num2str(result(:,1));
figure(5)
subplot(1,3,1)
bar(result(:,1),result(:,2))
title('用户最近一次缴费距今时长分布条形图','fontsize',12)
subplot(1,3,2)
pie(result(:,2))
title('用户最近一次缴费距今时长分布饼状图','fontsize',12)
legend(lab);
subplot(1,3,3)
plot(result(:,1),ms,'-*r');
title('用户最近一次缴费距今时长分布与信用分关系折线图','fontsize',12)
disp("用户实名制是否通过核实")
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

连续变量分析-用户近6个月平均消费值

    用户近6个月消费值也是一个重要的特征。不同的人具有不同的消费能力，往往消费能力高的用户的收入水平也较高，它具有比较高的可支配财富，所以他抵御风险的能力也较强，相比于用户当月账单总费用，近6个月消费值对应的时间周期更长，反映的情况也更稳定，信用分也应该越高。

figure(9)
shop6=data(:,9);
plot(shop6,target,'-og')
xlabel('用户近6个月平均消费值')
ylabel('信用分')
title('用户近6个月平均消费值与信用分之间关系的折线图','fontsize',12)
1
2
3
4
5
6

连续变量分析-用户年龄分布分析

    用户年龄分布也是一个十分影响结果的变量，对于年轻人来说，因为刚刚步入社会，容易遇到各种各样的生活压力，而这个时候他们自己还没有稳定的工作，所以十分容易发生债务违约的现象，但是对于中年人和老年人来说，他们已经成家立业，步入社会，遇到生活压力时候的抗风险能力也会越强。所以接下来我对用户年龄分布进行了分析。在Matlab中提供了ksdensity函数来让我们绘制核密度图，核密度图显示年龄的概率分布，能直观展现数据分布情况。

figure(11)
[f1, x1] = ksdensity(year);
t1 = area(x1, f1);t1.FaceColor='b'; t1.FaceAlpha=0.5;
title('用户年龄核密度分布图');
xlabel("年龄");
1
2
3
4
5

                                                                                                                            我们可以把年龄按照10岁为一个区间划分成十个区间，这样能更方便的观察到结果。划分成十个区间后，我们可以作出帕累托图，观察数据的累计分布的情况。在Matlab中我们可以使用pareto函数作出帕累托图。我们在按年龄段统计用户的数量的时候可以直接一次搞定，不需要出现10个if，因为用户的年龄数据除以10后取值得到的结果，就是对应以10为单位的年龄段的下标，直接把下标对应的数加1就行。

figure(10)
year=data(:,2);
result = zeros(1,10);
for i=1:46645
    result(1,floor(year(i,1)/10)+1)=result(1,floor(year(i,1)/10)+1)+1;
end
pareto(result);
xlabel('用户年龄')
ylabel('人数')
1
2
3
4
5
6
7
8
9

相关性分析

    为了直观地得到不同变量之间的关系，我们可以使用相关系数热力图来观察变量之间的相关系数。在Matlab中我们可以使用corr求得相关系数，然后使用colormap绘制出热力图。

figure(12)
rho = corr(data);
string_name={'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'};
xvalues = string_name;
yvalues = string_name;
h = heatmap(xvalues,yvalues,rho, 'FontSize',10, 'FontName','Times New Roman');
h.Title = 'Correlation Coefficient';
colormap(jet)
1
2
3
4
5
6
7
8

特征工程

特征创造

    在本问题中根据有些关于欠费和账单的特征能进行符合实际生活场景的特征创造，生成新的特征。对于缴费用户最近一次缴费金额和用户近6个月平均消费值这两个特征。我的理解是对于用户来说，如果用户缴纳的费用比用户账单金额高，说明该用户收入水平比较高足以付清该用户的账单。所以进行特征创造的时候可以使用缴费用户最近一次缴费金额减去用户近6个月平均消费值作为一个新的特征字段用来代表用户的收入水平支付当月消费的账单的能力，这个值是正数代表用户该月收入足以偿还当月的账单，这个值如果是负数代表用户该月的收入不足以偿还当月消费的账单。而且这个只越大代表用户偿还账单的能力越强，可以反映出该用户的信用情况是越好的。
    上面的特征创造我们是基于用户的收入和花销间的关系，除此之外我们可以考虑用户的花销本身。如果用户的花销越来越小，说明用户违约的可能性也越来越小。基于次我考虑使用用户账单当月总费用减去用户近6个月平均消费值作为一个新的特征字段用来代表用户花销的变化情况。这个值越小代表用户的花销的变化情况是越来越小的。进行特征创造的代码如下所示，在进行特征创造时我使用了向量化编程的思想方法，不使用for循环去遍历46635行的原始数据，从而提高程序运行效率：

disp('特征创造')
newFea = zeros(46645,1);
newFea(:,1)=data(:,8)-data(:,9);
data = [data newFea];
newFea(:,1)=data(:,10)-data(:,9);
data = [data newFea];
1
2
3
4
5
6

特征筛选

    有些字段的数据是完全一样的这对于建立机器学习模型没有任何帮助，反而会加重机器学习算法模型的负担。比如说我们之前讨论过的用户实名制是否通过核实字段，这个字段的数据都是1，并不会对区分不同用户的信用分带来任何帮助，所以我们去除该字段的数据。

disp('特征筛选')
data(:,1) = [];
1
2

特征标准化-zscore标准化

    许多机器学习的算法中目标函数的基础都是假设所有的特征都是零均值并且具有同一阶数上的方差。如果某个特征的方差比其他特征大几个数量级，那么它就会在学习算法中占据主导位置，会影响学习器的效率，导致学习器并不能像我们说期望的那样，从其他特征中学习。特征标准化是让不同维度之间的特征在数值上有一定比较性，可以大大提高分类器的准确性。但是对于回归算法来说，本身就需要同时考虑各个不同的特征，如果不同特征之间数量级相差太大，将会影响学习器学习的效率，所以我们应该对原始数据进行特征标准化。在Matlab中提供给我们一种函数zscore，这种函数支持了0标准化。

disp('特征标准化')
traindata = data;
traindata(:,28) = [];
newData = zscore(traindata);
thetarget = zscore(target);
1
2
3
4
5

特征降维-pca降维

    由于经过处理后的数据有多达30个特征，有很多特征只是噪声数据，所以我们应该要把它进行pca降维，把降维后的数据作为机器学习算法建模的原始数据。在Matlab中提供函数pca来让我们实现pca降维，我们在进行pca降维后要绘制出降维后的各个主成分的解释贡献率帕累托图来决定我们需要选取那几个主成分进行回归算法的建模。

figure(13)
disp('特征降维')
[COEFF,SCORE,latent,tsquared,explained,mu]=pca(newData);
data_PCA=newData*COEFF(:,1:10);
latent1=100*latent/sum(latent);%将latent总和统一为100，便于观察贡献率
pareto(latent1);%调用matla画图 pareto仅绘制累积分布的前95%，因此y中的部分元素并未显示
xlabel('Principal Component');
ylabel('Variance Explained (%)');
iris_pac=data_PCA(:,1:5) ;
1
2
3
4
5
6
7
8
9

数据集划分

    对于机器学习建模来说，通常情况下应该要划分成是三个数据集，训练集用来训练数据，测试集用来调节参数，验证集在未参与模型训练的情况下用来评估泛化能力。所以我先按3:1的大小关系划分出验证集和其他数据集。

assess_x = iris_pac1(size(iris_pac1,1)*0.75:end,:);
assess_y = target(size(target,1)*0.75:end,:);
iris_pac1 = iris_pac1(1:size(iris_pac1,1)*0.75,:);
target = target(1:size(target,1)*0.75,:);
1
2
3
4

决策树回归

    在Matlab中提供了fitrtree函数用来构建决策树的机器学习模型，我设置了参数CrossVal代表交叉验证，通过循环调整minleaf这个参数的取值来调整参数，然后绘制出了每个参数下的决策树模型对应的RMSE的值。

leafs = logspace(1,2,10);
rng('default')
N = numel(leafs);
err = zeros(N,1);
for n=1:N
        t =  fitrtree(iris_pac1,target,'CrossVal','On', 'MinLeaf',leafs(n));%交叉验证法估计算法精度
        err(n) = kfoldLoss(t);%计算误差
end
plot(leafs,err);
grid on
xlabel('Min Leaf Size');
ylabel('cross-validated error');
title('决策树minleaf调参')
OptimalTree = fitrtree(iris_pac1,target,'minleaf',46);
resuberror = resubLoss(OptimalTree) %衡量误差，默认均方差算法，此处可以设置损失函数
lossOpt = kfoldLoss(crossval(OptimalTree))
Ynew = predict(OptimalTree,iris_pac1);
trainResult = fix(Ynew);
newYnew = predict(OptimalTree,assess_x);
testResult = fix(newYnew);
rmse=0;
trainResult1 = trainResult;
testResult1 = testResult;
for i=1:34983
    rmse = rmse + (trainResult(i,1)-target(i,1))^2;
end
rmse=rmse/34983;
rmse1= 0;
for i=1:11662
    rmse1 = rmse1 + (testResult(i,1)-assess_y(i,1))^2;
end
rmse1=rmse1/11662;
disp("决策树回归训练集上的rmse:")
disp(rmse)
disp("决策树回归测试集上的rmse:")
disp(rmse1)
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

高斯回归

    在Matlab中提供了fitrgp函数用来进行高斯回归，我在高斯回归中使用的都是默认参数。

%高斯回归模型

t =  fitrgp(iris_pac1,target);
Ynew = predict(t,iris_pac1);
trainResult = fix(Ynew);
newYnew = predict(t,assess_x);
testResult = fix(newYnew);
rmse=0;
trainResult2 = trainResult;
testResult2 = testResult;
for i=1:34983
    rmse = rmse + (trainResult(i,1)-target(i,1))^2;
end
rmse=rmse/34983;
rmse1= 0;
for i=1:11662
    rmse1 = rmse1 + (testResult(i,1)-assess_y(i,1))^2;
end
rmse1=rmse1/11662;
disp("高斯回归训练集上的rmse:")
disp(rmse)
disp("高斯回归测试集上的rmse:")
disp(rmse1)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

线性回归

    在Matlab中提供了fitrlinear函数用来进行线性回归，我在线性回归中使用的也都是默认参数。

%线性回归模型

t =  fitrlinear(iris_pac1,target);
Ynew = predict(t,iris_pac1);
trainResult = fix(Ynew);
newYnew = predict(t,assess_x);
testResult = fix(newYnew);
rmse=0;
trainResult3 = trainResult;
testResult3 = testResult;
for i=1:34983
    rmse = rmse + (trainResult(i,1)-target(i,1))^2;
end
rmse=rmse/34983;
rmse1= 0;
for i=1:11662
    rmse1 = rmse1 + (testResult(i,1)-assess_y(i,1))^2;
end
rmse1=rmse1/11662;
disp("线性回归训练集上的rmse:")
disp(rmse)
disp("线性回归测试集上的rmse:")
disp(rmse1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

模型融合-基于回归树集成算法的stacking堆叠法

    Stacking堆叠法的基本思想是基于几个基本学习器的结果，在它们这些结果的上面一层再构造一层学习器。我stacking堆叠法选择的是回归树集成算法。

%模型融合-stacking/回归树集成
trainDatas = [trainResult1 trainResult2 trainResult3];
testDatas = [testResult1 testResult2 testResult3];

t =  fitrensemble(trainDatas,target);
Ynew = predict(t,trainDatas);
trainResult = fix(Ynew);
newYnew = predict(t,testDatas);
testResult = fix(newYnew);
rmse=0;
for i=1:34983
    rmse = rmse + (trainResult(i,1)-target(i,1))^2;
end
rmse=rmse/34983;
rmse1= 0;
for i=1:11662
    rmse1 = rmse1 + (testResult(i,1)-assess_y(i,1))^2;
end
rmse1=rmse1/11662;
disp("模型融合训练集上的rmse:")
disp(rmse)
disp("模型融合测试集上的rmse:")
disp(rmse1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

总结反思

    在结果中我们可以看到我尝试了好几种回归算法，只有决策树回归算法的效果是最好的，传统的线性回归和基于概率统计的高斯回归算法的效果都是比较差的，而基于模型融合技术的stacking堆叠法对于提升模型准确率有极大帮助。表现比决策树算法的结果要优秀很多。