数学建模：使用灰色预测和 SVR 构建财政收入模型

阳光穿透心脏的1/2处 2023-02-28 05:28 3阅读 0赞

灰色预测算法：[https://blog.csdn.net/weixin\_41987016/article/details/107448710][https_blog.csdn.net_weixin_41987016_article_details_107448710]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 1][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 2][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 3][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 4][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 5][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 6][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 7][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 8][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 9][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 10][]

完整代码：

#GM.py
    def GM(x0): #自定义灰色预测函数
      import numpy as np
      x1 = x0.cumsum() #1-AGO序列
      z1 = (x1[:len(x1)-1] + x1[1:])/2.0 #紧邻均值（MEAN）生成序列
      z1 = z1.reshape((len(z1),1))
      B = np.append(-z1, np.ones_like(z1), axis = 1)
      Yn = x0[1:].reshape((len(x0)-1, 1))
      [[a],[b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Yn) #计算参数
      f = lambda k: (x0[0]-b/a)*np.exp(-a*(k-1))-(x0[0]-b/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.6*745*x0.std()).sum()/len(x0)
      return f, a, b, x0[0], C, P #返回灰色预测函数、a、b、首项、方差比、小残差概率

#pearson相关系数
    import numpy as np
    import pandas as pd
    inputfile = 'data.csv' ## 输入的数据文件
    data = pd.read_csv(inputfile) ## 读取数据
    ## 保留两位小数
    # print('相关系数矩阵为：',np.round(data.corr(method = 'pearson'), 2))
    data = np.round(data.corr(method = 'pearson'), 2)
    outputfile = '相关系数矩阵.csv'
    data.to_csv(outputfile)  #存储数据
    #查看正相关和负相关
    data.sort_values(by=['y'],ascending=False,inplace=True)
    data.reset_index(drop = True)
    print(data['y'])

#使用Lasso回归选取财政收入预测的关键特征
    # 代码 8-2
    import numpy as np
    import pandas as pd
    from sklearn.linear_model import Lasso
    inputfile = 'data.csv' #输入的数据文件
    data = pd.read_csv(inputfile) #读取数据
    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)
    
    outputfile = 'new_reg_data.csv'  #输出的数据文件
    new_reg_data = data.iloc[:, mask]  #返回相关系数非零的数据
    new_reg_data.to_csv(outputfile)  #存储数据
    print('输出数据的维度为：',new_reg_data.shape)  #查看输出数据的维度

#构建灰色预测模型，并预测2014年和2015年的政财收入
    import numpy as np
    import pandas as pd
    from GM import GM #引入自编的灰色预测函数
    inputfile = 'new_reg_data.csv' #输入的数据文件
    inputfile1 = 'data.csv' #输入的数据文件
    new_reg_data = pd.read_csv(inputfile) #读取经过特征选择后的数据
    data = pd.read_csv(inputfile1) #读取经过特征选择后的数据
    new_reg_data.index = range(1994,2014)
    new_reg_data.loc[2014] = None
    new_reg_data.loc[2015] = None
    l = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']
    for i in l:
        #new_reg_data.loc[range(1994,2014),i]获取l(i)的列数据
        f = GM(new_reg_data.loc[range(1994,2014),i].as_matrix())[0]
        ##将2014和2015的预测结果放入new_reg_data的l对应列表里
        new_reg_data.loc[2014,i] = f(len(new_reg_data)-1)#2014年预测结果
        new_reg_data.loc[2015,i] = f(len(new_reg_data)) ##2015年预测结果
        new_reg_data[i] = new_reg_data[i].round(2) ## 保留两位小数
    outputfile = 'new_reg_data_GM.xls' ## 灰色预测后保存的路径
    y = list(data['y'].values) ## 提取财政收入列，合并至新数据框中
    y.extend([np.nan,np.nan])
    new_reg_data['y'] = y
    new_reg_data.to_excel(outputfile) ## 结果输出
    print('预测结果为：',new_reg_data.loc[2014:2015,:]) ##预测结果展示
    for i in l:
        C = GM(new_reg_data.loc[range(1994,2014),i].as_matrix())[4]
        p = GM(new_reg_data.loc[range(1994,2014),i].as_matrix())[5]
        print('   P   C')
        print(i,p,C)

import numpy as np
    from sklearn.svm import LinearSVR
    import matplotlib.pyplot as plt
    inputfile = 'new_reg_data_GM11.xls' #灰色预测后保存的路径
    data = pd.read_excel(inputfile) #读取数据
    data.index = range(1994,2016)
    feature = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13'] 
    data_train = data.loc[range(1994,2014)].copy()#取2014年前的数据建模
    data_mean = data_train.mean()
    data_std = data_train.std()
    data_train = (data_train - data_mean)/data_std #数据标准化
    
    x_train = data_train[feature].as_matrix() #特征数据
    y_train = data_train['y'].as_matrix() #标签数据
    linearsvr = LinearSVR()   #调用LinearSVR()函数
    linearsvr.fit(x_train,y_train)
    x = ((data[feature] - data_mean[feature])/ \
    data_std[feature]).as_matrix()  #预测，并还原结果。
    data[u'y_pred'] = linearsvr.predict(x) * \
    data_std['y'] + data_mean['y']
    ## SVR预测后保存的结果
    outputfile = 'new_reg_data_GM11_revenue.xls'
    data.to_excel(outputfile)
    print('真实值与预测值分别为：',data[['y','y_pred']])
    
    print('预测图为：',data[['y','y_pred']].plot(subplots = True,
          style=['b-o','r-*']))

#!usr/bin/env python
    #encoding:utf-8
    from __future__ import division
     
     
    '''
    功能：计算回归分析模型中常用的四大评价指标
    '''
     
    from sklearn.metrics import explained_variance_score, mean_absolute_error, median_absolute_error, r2_score
    def calPerformance(y_true,y_pred):
        '''
        模型效果指标评估
        y_true：真实的数据值
        y_pred：回归模型预测的数据值
        explained_variance_score：解释回归模型的方差得分，其值取值范围是[0,1]，越接近于1说明自变量越能解释因变量
        的方差变化，值越小则说明效果越差。
        mean_absolute_error：平均绝对误差（Mean Absolute Error，MAE），用于评估预测结果和真实数据集的接近程度的程度
        ，其其值越小说明拟合效果越好。
        r2_score：判定系数，其含义是也是解释回归模型的方差得分，其值取值范围是[0,1]，越接近于1说明自变量越能解释因
        变量的方差变化，值越小则说明效果越差。
        '''
        model_metrics_name=[mean_absolute_error, median_absolute_error,explained_variance_score, r2_score]  
        tmp_list=[]  
        for one in model_metrics_name:  
            tmp_score=one(y_true,y_pred)  
            tmp_list.append(tmp_score)  
        print(['mean_absolute_error','median_absolute_error','explained_variance_score','r2_score'])
        print(tmp_list)
        return tmp_list
     
    if __name__=='__main__':
        inputfile = 'predict_data.csv' #输入的数据文件
        data = pd.read_csv(inputfile)
        data.drop(data[np.isnan(data['y'])].index, inplace=True)
        y_pred = data['y_pred']
        y_true = data['y']
        calPerformance(y_true,y_pred)

[https_blog.csdn.net_weixin_41987016_article_details_107448710]: https://blog.csdn.net/weixin_41987016/article/details/107448710
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70]: /images/20230209/1c087c65bd414a739f4ff20161ab9822.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 1]: /images/20230209/a89b7ddbcd1b487abd7e28542da269aa.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 2]: /images/20230209/e799f3fbbc824be1ad639f637d8eccd0.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 3]: /images/20230209/87ab5dcea2eb42e782f2f43b9af8a585.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 4]: /images/20230209/5508522fabbb4b2aa4b3113c34271222.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 5]: /images/20230209/0c6d1197a9c047de9efd652cf3d7f2a7.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 6]: /images/20230209/69d8c6d28b7b44daae4d6773c7d9b942.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 7]: /images/20230209/55b01315c6f445768f7d2c2bc080dc9d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 8]: /images/20230209/1c3245ef1bb94d44ba4893686a9a082f.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 9]: /images/20230209/aa8542e0ef7f470aba3e49b070c724e5.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTk4NzAxNg_size_16_color_FFFFFF_t_70 10]: /images/20230209/c034fdb3e8f54c21a275326c3f8a010c.png