【数据挖掘案例】财政收入影响因素分析及预测模型

案例来源自《Python 数据分析与挖掘实战》第 13 章：财政收入影响因素分析及预测模型

案例目的：预测财政收入

1 数据的基本描述性分析

数据

字段含义
社会从业人数(x1 )、 在岗职工工资总额（x2）、社会消费品零售总额(x3)、城镇居民人均可支配收人(x4)、 城镇居民人均消费性支出(x5)、年末总人口(x6)、全社会固定资产投资额(x7)、地区生产总值（x8）、第一产业产值（x9）、税收（x10）、居民消费价格指数（x11）、第三产业与第二产业产值比（x12）地区生产总值(x8)和居民消费水平(x13)。

导包与读取数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pylab import mpl
# 正常显示中文标签
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 正常显示负号
mpl.rcParams['axes.unicode_minus'] = False
# 禁用科学计数法
pd.set_option('display.float_format', lambda x: '%.2f' % x) 

# 读入数据
data = pd.read_excel('data1.xlsx',index_col=0)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

数据的基本情况

data.shape # (20, 14)
data.info()
'''
<class 'pandas.core.frame.DataFrame'>
Int64Index: 20 entries, 1994 to 2013
Data columns (total 14 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x1      20 non-null     int64  
 1   x2      20 non-null     float64
 2   x3      20 non-null     float64
 3   x4      20 non-null     float64
 4   x5      20 non-null     float64
 5   x6      20 non-null     int64  
 6   x7      20 non-null     float64
 7   x8      20 non-null     float64
 8   x9      20 non-null     float64
 9   x10     20 non-null     float64
 10  x11     20 non-null     float64
 11  x12     20 non-null     float64
 12  x13     20 non-null     int64  
 13  y       20 non-null     float64
dtypes: float64(11), int64(3)
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# 描述性分析
data.describe().T
1
2

或更简洁一点

r = [data.min(), data.max(), data.mean(), data.std()]
r = pd.DataFrame(r, index=['Min', 'Max', 'Mean', 'STD']).T
r = np.round(r, 2)
r
1
2
3
4

结果

            Min         Max        Mean         STD
x1   3831732.00  7599295.00  5579519.95  1262194.72
x2       181.54     2110.78      765.04      595.70
x3       448.19     6882.85     2370.83     1919.17
x4      7571.00    42049.14    19644.69    10203.02
x5      6212.70    33156.83    15870.95     8199.77
x6   6370241.00  8323096.00  7350513.60   621341.85
x7       525.71     4454.55     1712.24     1184.71
x8       985.31    15420.14     5705.80     4478.40
x9        60.62      228.46      129.49       50.51
x10       65.66      852.56      340.22      251.58
x11       97.50      120.00      103.31        5.51
x12        1.03        1.91        1.42        0.25
x13     5321.00    41972.00    17273.80    11109.19
y         64.87     2088.14      618.08      609.25
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

变量的分布情况

from sklearn.preprocessing import MinMaxScaler
#实现归一化
scaler = MinMaxScaler() #实例化
scaler = scaler.fit(data) #fit，在这里本质是生成min(x)和max(x)
data_scale = pd.DataFrame(scaler.transform(data)) #通过接口导出结果
data_scale.columns = data.columns

import joypy
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
import seaborn as sns
fig, axes = joypy.joyplot(data_scale, alpha=.5, color='#FFCC99')#连续值的列为一个"脊"
1
2
3
4
5
6
7
8
9
10
11
12
13
14

data_scale.plot()
1

---

相关性分析

pear = np.round(data.corr(method = 'pearson'), 2)
pear
1
2

plt.figure(figsize=(12,12))
sns.heatmap(data.corr(), center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5},annot=True, fmt='.1f')
#设置x轴
plt.xticks(fontsize=15)
#设置y轴
plt.yticks(fontsize=15)
plt.tight_layout()
plt.savefig('a.png')
1
2
3
4
5
6
7
8
9

由图可知，居民消费价格指数（x11）与财政收入的线性关系不显著，而且呈现负相关。其余变量均与财政收入呈现高度的正相关关系。

2 数据的预处理

变量的筛选

分析方法的选择：以往对财政收入的分析会使用 多元线性回归模型和最小二乘估计方法来估计回归模型的系数，通过系数能否通过检验来检验它们之间的关系，但这样的结果对数据依赖程度很大，并且求得的往往只是局部最优解，后续的检验可能会失去应有的意义。因此本案例运用Adaptive-Lasso变量选择方法来研究

对于Lasso，这里参考书中的理论知识

Lasso全称最小绝对收缩和选择算子，和岭回归一样，Lasso是被创造来作用于多重共线性问题的算法，不过Lasso使用的是系数 $\beta$ 的L1范式（L1范式则是系数 $\beta$ 的绝对值）乘以正则化系数

不太严谨说法是，Lasso与岭回归非常相似，都是利用正则项来对原本的损失函数形成一个惩罚，以此来防止多重共线性。

Lasso 虽然是为了限制多重共线性被创造出来的，但其实并不使用它来抑制多重共线性，反而接受了它在其他方面的优势，

由于，L2正则化只会将系数压缩到尽量接近0，但L1正则化主导稀疏性，因此会将系数压缩到0。

因此，Lasso成为了线性模型中的特征选择工具首选。

Lasso回归的特点是在拟合广义线性模型的同时进行变量筛选和复杂度调整。 因此，不论目标因变量是连续的，还是二元或者多元离散的，都可以用Lasso回归建模然后预测。

这里的变量筛选是指不把所有的变量都放入模型中进行拟合，而是有选择的把变量放入模型从而得到更好的性能参数。

—关于Adaptive-Lasso变量选择方法（参考书籍）

Adaptive_lasso算法是近些年来被广泛应用于参数估计于变量选择的方法之一。Adaptive_Lasso算法能够解决最小二乘法和逐步回归局部最优解的不足，这是它的优点之一。


Adaptive_lasso算法计算出某变量的特征值非零，则表示该变量对预测变量存在较大影响，而如果某变量的特征值为零，则表示该变量对预测变量影响很小。

---

Lasso变量选择模型

这里没有找到AdaptiveLasso这个函数，用Lasso代替

from sklearn.linear_model import Lasso
model = Lasso(alpha=0.1, max_iter=100000)
model.fit(data.iloc[:, 0:13], data['y'])
q=model.coef_#各特征的系数
q=pd.DataFrame(q,index=data.columns[:-1])
q
'''
        0
x1  -0.00
x2  -0.59
x3   0.44
x4  -0.13
x5   0.17
x6   0.00
x7   0.27
x8   0.03
x9  -7.56
x10 -0.09
x11  3.38
x12  0.00
x13 -0.01
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

计算出某变量的特征值非零，则表示该变量对预测变量存在较大影响，而如果某变量的特征值为零，则表示该变量对预测变量影响很小。

调整参数值（参照：https://blog.csdn.net/weixin_43746433/article/details/100047231）

from sklearn.linear_model import Lasso
lasso = Lasso(1000)  #调用Lasso()函数，设置λ的值为1000
lasso.fit(data.iloc[:,0:13],data['y'])
print('相关系数为：',np.round(lasso.coef_,5))  #输出结果，保留五位小数
## 计算相关系数非零的个数
print('相关系数非零个数为：',np.sum(lasso.coef_ != 0))
mask = lasso.coef_ != 0  #返回一个相关系数是否为零的布尔数组
print('相关系数是否为零：',mask)

new_reg_data = data.iloc[:,:13].iloc[:,mask]  #返回相关系数非零的数据
new_reg_data = pd.concat([new_reg_data,data.y],axis=1)
new_reg_data.to_excel('new_reg_data.xlsx')
1
2
3
4
5
6
7
8
9
10
11
12

相关系数为：
 [-1.8000e-04 -0.0000e+00  1.2414e-01 -1.0310e-02  6.5400e-02  1.2000e-04
  3.1741e-01  3.4900e-02 -0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00
 -4.0300e-02]
相关系数非零个数为： 8
相关系数是否为零： [ True False  True  True  True  True  True  True False False False False  True]
1
2
3
4
5
6

根据非零的系数，最终通过Lasso筛选出的变量如下

变量筛选好之后，接下来的开始建模

3 建立财政收入预测模型

这里使用的预测方法是，灰色预测 + 神经网络预测 的组合模型

基本思路是

（1）首先是在前面的变量选择结果基础上，对单个选定的影响因素建立灰色预测模型，得到它们在2014年及2015年的预测值。

（2）然后，通过神经网络对历史数据（$X_{1994-2003}$，$Y_{1994-2003}$）建立训练模型。

（3）把$X_{1994-2003}$以及灰色预测的数据结果$X_{2004-2005}$一并代人训练好的模型中，得到预测值。

---

灰色预测模型

灰色模型学习资料（自己推荐自己*-*）

灰色预测函数

def GM11(x0): #自定义灰色预测函数
    import numpy as np
    x1 = x0.cumsum() # 生成累加序列
    z1 = (x1[:len(x1)-1] + x1[1:])/2.0 # 生成紧邻均值（MEAN）序列，比直接使用累加序列好，共 n-1 个值
    z1 = z1.reshape((len(z1),1))
    B = np.append(-z1, np.ones_like(z1), axis = 1)    # 生成 B 矩阵
    Y = x0[1:].reshape((len(x0)-1, 1))    # Y 矩阵
    [[a],[u]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)    #计算参数
    f = lambda k: (x0[0]-u/a)*np.exp(-a*(k-1))-(x0[0]-u/a)*np.exp(-a*(k-2))    #还原值
    delta = np.abs(x0 - np.array([f(i) for i in range(1,len(x0)+1)]))    # 计算残差
    C = delta.std()/x0.std()
    P = 1.0*(np.abs(delta - delta.mean()) < 0.6745*x0.std()).sum()/len(x0)
    return f, a, u, x0[0], C, P #返回灰色预测函数、a、b、首项、方差比、小残差概率
1
2
3
4
5
6
7
8
9
10
11
12
13

---

data.index = range(1994, 2014)
data.loc[2014] = None
data.loc[2015] = None
1
2
3

# 模型精度评价
# 被lasso筛选出来的6个变量
l = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']
for i in l:
    GM = GM11(data[i][list(range(1994, 2014))].values)
    f = GM[0]
    c = GM[-2]
    p = GM[-1]
    data[i][2014] = f(len(data)-1)
    data[i][2015] = f(len(data))
    data[i] = data[i].round(2)
    if (c < 0.35) & (p > 0.95):
        print('对于模型{}，该模型精度为---好'.format(i))
    elif (c < 0.5) & (p > 0.8):
        print('对于模型{}，该模型精度为---合格'.format(i))
    elif (c < 0.65) & (p > 0.7):
        print('对于模型{}，该模型精度为---勉强合格'.format(i))
    else:
        print('对于模型{}，该模型精度为---不合格'.format(i))

data[l+['y']].to_excel('data2.xlsx')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

对于模型x1，该模型精度为---好
对于模型x3，该模型精度为---好
对于模型x4，该模型精度为---好
对于模型x5，该模型精度为---好
对于模型x6，该模型精度为---好
对于模型x7，该模型精度为---好
对于模型x8，该模型精度为---好
对于模型x13，该模型精度为---好
1
2
3
4
5
6
7
8

预测值如下：

---

神经网络预测模型

下面用历史数据建立神经网络模型

其参数设置为误差精度107，学习次数10000次，神经元个数为Lasso变量选择方法选择的变量个数8。

'''神经网络'''
data2 = pd.read_excel('data2.xlsx', index_col=0)
# 提取数据
feature = list(data2.columns[:len(data2.columns)-1]) # ['x1', 'x2', 'x3', 'x4', 'x5', 'x7']
train = data2.loc[list(range(1994, 2014))].copy()
mean = train.mean()
std = train.std() 
train = (train - mean) / std    # 数据标准化，这里使用标准差标准化
x_train = train[feature].values
y_train = train['y'].values

# 建立神经网络模型 
from keras.models import Sequential
from keras.layers.core import Dense, Activation

model = Sequential()
model.add(Dense(input_dim=8, units=12))
model.add(Activation('relu'))
model.add(Dense(input_dim=12, units=1))
model.compile(loss='mean_squared_error', optimizer='adam')
model.fit(x_train, y_train, epochs=10000, batch_size=16)
model.save_weights('net.model')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

训练好模型后，将$X_{1994-2005}$全部带入模型中作预测，结果如下：

# 将整个变量矩阵标准化
x = ((data2[feature] - mean[feature]) / std[feature]).values
# 预测，并还原结果
data2['y_pred'] = model.predict(x) * std['y'] + mean['y']
data2.to_excel('data3.xlsx')
1
2
3
4
5

预测结果

绘制真实值与预测值之间的折线图

import matplotlib.pyplot as plt
p = data2[['y', 'y_pred']].plot(style=['b-o', 'r-*'])
p.set_ylim(0, 2500)
p.set_xlim(1993, 2016)
plt.show()
1
2
3
4
5

从结果中，比较预测值与真实值基本高度吻合

为了与神经网络预测结果有一个对比，下面使用其他预测模型查看其结果如何

from sklearn.linear_model import LinearRegression # 线性回归
from sklearn.neighbors import KNeighborsRegressor # K近邻回归
from sklearn.neural_network import MLPRegressor # 神经网络回归
from sklearn.tree import DecisionTreeRegressor # 决策树回归
from sklearn.tree import ExtraTreeRegressor # 极端随机森林回归
from xgboost import XGBRegressor # XGBoot
from sklearn.ensemble import RandomForestRegressor # 随机森林回归
from sklearn.ensemble import AdaBoostRegressor  # Adaboost 集成学习
from sklearn.ensemble import GradientBoostingRegressor # 集成学习梯度提升决策树
from sklearn.ensemble import BaggingRegressor # bagging回归
from sklearn.linear_model import ElasticNet

from sklearn.metrics import explained_variance_score,\
mean_absolute_error,mean_squared_error,\
median_absolute_error,r2_score

models=[LinearRegression(),KNeighborsRegressor(),MLPRegressor(alpha=20),DecisionTreeRegressor(),ExtraTreeRegressor(),XGBRegressor(),RandomForestRegressor(),AdaBoostRegressor(),GradientBoostingRegressor(),BaggingRegressor(),ElasticNet()]
models_str=['LinearRegression','KNNRegressor','MLPRegressor','DecisionTree','ExtraTree','XGBoost','RandomForest','AdaBoost','GradientBoost','Bagging','ElasticNet']


data2 = pd.read_excel('data2.xlsx', index_col=0)
# 提取数据
feature = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']
train = data2.loc[list(range(1994, 2014))].copy()
mean = train.mean()
std = train.std() 
train = (train - mean) / std    # 数据标准化，这里使用标准差标准化
x_train = train[feature].values
y_train = train['y'].values
# 将整个变量矩阵标准化
x = ((data2[feature] - mean[feature]) / std[feature]).values


for name,model in zip(models_str,models):
    print('开始训练模型：'+name)
    model=model   #建立模型
    a = 'y_pred_'+ name
    data2[a] = model.fit(x_train,y_train).predict(x) * std['y'] + mean['y']
    df=data2[:-2]
    print('平均绝对误差为：',mean_absolute_error(df['y'].values,df[a].values))
    print('均方误差为：',mean_squared_error(df['y'],df[a]))
    print('中值绝对误差为：',median_absolute_error(df['y'],df[a]))
    print('可解释方差值为：',explained_variance_score(df['y'],df[a]))
    print('R方值为：',r2_score(df['y'],df[a]))
    print('*-*'*15)
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

结果

开始训练模型：LinearRegression
平均绝对误差为： 32.624948062446734
均方误差为： 1596.8011348674995
中值绝对误差为： 29.89818813145061
可解释方差值为： 0.9954717568606967
R方值为： 0.9954717568606967
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：KNNRegressor
平均绝对误差为： 63.916899999999984
均方误差为： 17521.95548179999
中值绝对误差为： 26.000000000000043
可解释方差值为： 0.9517094362717755
R方值为： 0.950310860278652
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：MLPRegressor
平均绝对误差为： 235.42036857561698
均方误差为： 76090.58720607273
中值绝对误差为： 196.3296364219621
可解释方差值为： 0.7845740946711923
R方值为： 0.7842206697141109
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：DecisionTree
平均绝对误差为： 2.7711166694643907e-14
均方误差为： 1.0551803468236254e-26
中值绝对误差为： 0.0
可解释方差值为： 1.0
R方值为： 1.0
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：ExtraTree
平均绝对误差为： 2.7711166694643907e-14
均方误差为： 1.0551803468236254e-26
中值绝对误差为： 0.0
可解释方差值为： 1.0
R方值为： 1.0
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：XGBoost
平均绝对误差为： 0.5735757446289022
均方误差为： 0.5008986120857999
中值绝对误差为： 0.6280151367187159
可解释方差值为： 0.9999985795409094
R方值为： 0.9999985795408995
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：RandomForest
平均绝对误差为： 24.30998000000009
均方误差为： 2541.844172424009
中值绝对误差为： 8.976250000000078
可解释方差值为： 0.9932099301805153
R方值为： 0.9927917834076989
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：AdaBoost
平均绝对误差为： 13.087633333333354
均方误差为： 386.4120509666666
中值绝对误差为： 7.9691666666666805
可解释方差值为： 0.9991155201544625
R方值为： 0.9989042043617541
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：GradientBoost
平均绝对误差为： 0.023211721877213876
均方误差为： 0.0014299207277676446
中值绝对误差为： 0.011657245535729999
可解释方差值为： 0.9999999959449999
R方值为： 0.9999999959449999
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：Bagging
平均绝对误差为： 34.284249999999986
均方误差为： 4470.224954049998
中值绝对误差为： 7.295999999999992
可解释方差值为： 0.9874093338106849
R方值为： 0.9873232395460447
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
开始训练模型：ElasticNet
平均绝对误差为： 294.6018692079905
均方误差为： 118127.42764924551
中值绝对误差为： 249.10992216728664
可解释方差值为： 0.665011689849139
R方值为： 0.6650116898491392
*-**-**-**-**-**-**-**-**-**-**-**-**-**-**-*
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
66
67
68
69
70
71
72
73
74
75
76
77

注：存在严重的过拟合问题

看一下，线性回归、XGBoost 预测值对比图

线性回归

XGBoost

这里做法存在一定的问题，没有划分训练集和测试集，存在严重过拟合，模型失效。

下面试一下

线支持向量机线性回归（LinearSVR）

from sklearn.svm import LinearSVR
linearsvr = LinearSVR().fit(x_train,y_train)  #调用LinearSVR()函数
# 预测，并还原结果
data2['y_pred_linearsvr'] = linearsvr.predict(x) * std['y'] + mean['y']
1
2
3
4

相比其他回归模型，神经网络预测模型在这里效果更好一点

至此。