【数学建模】熵值法确定权重及Python实现

2018-09-22 本文已影响167人 crossous

id=56316785

关键词：熵值法、Python、pandas

一、前言

本文的目的是用Python和类对熵值法确定权重的一系列行为进行封装
注意，本文的事例是运行在jupyter-notebook中，也可以直接运行到其他编辑器中。如果你看到直接输出变量（没用print）就出结果，这是jupyter的功能，你可以用print包裹代码最后的变量，能做到类似的效果。
注意，本文会讲核心代码按类封装，中间还有调试过程代码，如果不熟悉python，请注意不要把构建类的代码和调试代码混在一起，如果实在弄不懂，建议先搞清楚python面向对象（只用了解类就好）。
但因为本文是代码的封装，只要搞清楚传递数据的类型，就可以轻松移植代码用在其他题目上，不需要改变类内的代码。

二、熵值法估计权重原理

部分内容转自熵值法估计权重的详解并对原有内容进行删改，信息论和熵的原理直接参考原文。

熵值法估计权重模型

（1）选取 $n$ 个单位（省、市、国家），每个单位选取 $m$ 个指标作为变量，符号 $x_{ij}$ 表示第 $i$ 个单位的第 $j$ 个指标的数值 $(i=1,2,...,n;\quad j=1,2,...,m)$ 。
（2）指标的归一化处理：异质指标同质化，由于各项指标的计量单位并不统一，因此在用它们计算综合指标前，先要对它们进行标准化处理，即把指标的绝对值转化为相对值，并令 $x_{ij}=|x_{ij}|$ ，从而解决各项不同质指标值的同质化问题。而且，由于正向指标和负向指标数值代表的含义不同（正向指标数值越高越好，负向指标数值越低越好），因此，对于高低指标我们用不同的算法进行数据标准化处理。其具体方法如下:
正项指标： $x_{ij}^{'} = \frac{x_{ij}-min\{x_{1j},...,x_{nj}\}}{max\{x_{1j},...,x_{nj}\}-min\{x_{1j},...,x_{nj}\}}$ 负项指标： $x_{ij}^{'} = \frac{max\{x_{1j},...,x_{nj}\}-x_{ij}}{max\{x_{1j},...,x_{nj}\}-min\{x_{1j},...,x_{nj}\}}$ 则 $x_{ij}^{'}$ 为第 $i$ 个国家的第 $j$ 个指标的数值 $(i=1,2,...,n;\quad j=1,2,...,m)$ 。为了方便起见，归一化后的数据仍记为 $x_{ij}$ 。
（3）计算第 $j$ 项指标下第 $i$ 个样本值占该指标的比重： $p_{ij} = \frac{x_{ij}}{\sum_{i=1}^nx_{ij}}, \quad i=1,...,n,\quad j=1,...,m$ （4）计算第 $j$ 项指标的熵值： $e_j = -k\sum_{i=1}^np_{ij}\ln{(p_{ij})},\quad j = 1,..., m$ 其中 $k=\frac{1}{\ln\left( n \right)} > 0$ ，满足 $e_j \geq 0$ 。
（5）计算信息熵冗余度（差异）： $d_j = 1-e_j,\quad j=1,...,m$ （6）计算各项指标的权值： $w_j=\frac{d_j}{\sum_{j=1}^{m}d_j},\quad j=1,2,...,m$ （7）计算各样本的综合得分： $s_i = \sum_{j=1}^mw_j · p_{ij},\quad i=1,...,n$

三、按理及Python代码实现

1.题目

近年东北三省的GDP增速在全国各省中垫底，引起了全国的关注。搜集相关数据资料，建立模型对黑龙江省的经济发展现状做出量化评价。

2.解题思路

做这道问题的时候我特地去网上搜索东北经济状况发展不好的原因，其中有地理位置、资源产出、行政人员办事效率等多种因素共同的作用；不过这里我们是要量化分析，只要拿出各省的数据拿来用熵值法评价一番吧。

3.数据采集

首先是数据的采集，原题并没有给数据，此时我们就上网去搜，此类资料在国家统计局有不少，还提供csv下载；不过我是在百度文库下载的，此时放上数据：

2015年全年全国各省(市区)主要经济指标排名

地区,GDP总量,GDP总量增速,人口总量,人均GDP,人均GDP增速,地方财政收入总额,地方财政收入总额增速,固定资产投资,固定资产投资增速,社会消费品零售总额增速,进出口总额,进出口总额增速,城镇居民人均可支配收入,城镇居民人均可支配收入增速,农村居民人均可支配收入,农村居民人均可支配收入增速
广东,72812.55,8,10430.3,6.98,7.2,9364.76,16.2,29950.48,15.9,10.1,10229.5,-5,34757,8.1,13360,9.1
江苏,70116.4,8.5,7920,8.85,8.3,8028.29,11,45905.17,10.5,10.3,5456.1,-3.2,37173,8.2,16257,8.7
山东,63002.3,8,9579.3,6.58,7.5,5529.2,10,47381.46,13.9,10.6,2417.4,-12.7,31545,8,12930,8.8
浙江,42886.5,8,5613.7,7.64,4.4,4809.53,7.8,26664.72,13.2,10.9,3474.1,-2.2,43714,8.2,21125,9
河南,37010.25,8.3,9388,3.94,8.6,3009.6,9.9,34951.28,16.5,12.4,738.4,13.6,25576,8,10853,8.9
四川,30103.1,7.9,8050,3.74,7.9,3329.1,8.8,24965.56,10.2,12,515.9,-26.5,26205,8.1,10247,9.6
河北,29806.1,6.8,7185.4,4.15,4.4,2648.5,8.3,28905.74,10.6,9.4,514.8,-14,26152,8.3,11051,8.5
湖北,29550.19,8.9,5723.77,5.16,8.6,3005.39,17.1,26086.42,16.2,12.3,456.1,6,27051,8.8,11844,9.2
湖南,29047.2,8.6,7119.34,4.08,7.8,2513.1,11.1,24324.17,18.2,12.1,293.7,-4.8,28838,8.5,7675,9.8
辽宁,28743,3,4374.6,6.57,,2125.6,-33.4,17640.37,-27.8,7.7,960.9,-15.7,31126,,12057,
福建,25979.82,9,3689.4,7.04,8.7,2544.08,7.7,20973.98,17.4,12.4,1693.8,-4.5,33275,8.3,13793,9
上海,24964.99,6.9,2415.2,10.34,0.2,5519.5,13.3,6349.39,5.6,8.1,4492.4,-3.7,52962,8.4,23205,9.5
北京,22968.6,6.9,2170.5,10.58,1,4723.86,12.3,7446.02,8.3,7.3,3196.2,-23.1,52859,8.9,20569,9
安徽,22005.6,8.7,5950.05,3.7,8.6,2454,10.6,23803.93,12,12,488.1,-0.8,26936,8.4,10821,9.1
陕西,18171.86,8,3732.74,4.87,7,2059.87,12.1,18231.03,8.3,11.1,305,11.5,26420,8.4,8689,9.5
内蒙古,18032.79,7.7,2470.63,7.3,8.6,1963.5,6.5,13529.15,0.1,8,127.5,-12.4,30594,7.9,10776,8
广西,16803.12,8.1,4602.66,3.65,7.9,1515.08,6.5,15654.95,17.8,10,511.6,26.2,26416,7.1,9467,9
江西,16723.8,9.1,4503.93,3.7,9.2,2165.5,15.1,16993.9,16,11.4,426.1,-0.3,26500,9,11139,10.1
天津,16538.19,9.3,1516.81,10.9,9.3,2666.99,11.6,11814.57,12.6,10.7,1143.5,-14.6,34101,8.2,18482,8.6
重庆,15719.72,11,3016.55,5.21,10.8,2155.1,12.1,14208.15,17,12.5,749.4,-21.5,27239,8.3,10505,10.7
黑龙江,15083.7,5.7,3831.22,3.94,0.4,1165.2,-10.4,9884.28,3.6,8.9,209.8,-46.1,24203,7,11095,6.1
吉林,14274.11,6.5,2746,5.2,5.3,1229.29,2.2,12508.59,12.6,9.3,189.3,-28.3,24901,7.2,11326,5.1
云南,13717.88,8.7,4596.6,2.98,6.7,1808.14,6.5,13069.39,18,10.2,244,-17.6,26373,8.5,8242,10.5
山西,12802.58,3.1,3593.3,3.56,-2.8,1642.2,-9.8,13744.59,14.8,5.5,147.2,-9.3,25828,7.3,9454,7.3
贵州,10502.56,10.7,3475,3.02,9.9,1503.35,10,10676.7,21.6,11.8,123,14.2,24580,9,7387,10.7
新疆,9324.8,8.8,2181.33,4.27,5.2,1331,3.8,10525.42,10.1,7,196.8,-28.9,26275,13.2,9425,8
甘肃,6790.32,8.1,2557.53,2.66,6.8,743.9,10.6,8626.6,11.2,9,80.1,-7.3,23767,9,6936,10.5
海南,3702.8,7.8,867,4.27,5.1,627.7,8.7,3355.4,10.4,8.2,139.6,-12,26356,7.6,10858,9.5
宁夏,2911.77,8,630.14,4.62,,373.74,10,3426.42,10.7,7.1,37.9,-30.3,25186,8.2,9119,8.4
青海,2417.05,8.2,562.67,4.3,7.6,267.12,6.1,3144.17,12.7,11.3,19.3,12.6,24542,10,7933,8.9
西藏,1026.39,11,281,3.65,14.6,176,31,1295.68,21.2,12,9.1,-59.5,25457,15.6,8244,12

这里面有不少缺失数据，例如辽宁：

缺少人均GDP增速、城镇居民人均可支配收入增速、农村居民人均可支配收入增速。对于缺失数据，最简单的办法是直接去掉，但也可以去国家统计局去补全。

搜索辽宁并选择相关报表

找到需要的选项，并选择需要的时间点、时间段我们需要的是增速，这里好像没有，那么计算就好，用2015年的数据除以2014年的数据再减一。就这样补全了数据。

4.程序实现

（1）读入数据处理

首先，导入需要的包并读入数据

import pandas as pd
f = open('2015年全年全国各省(市区)主要经济指标排名.csv', encoding='utf8')
df = pd.read_csv(f)
df = df.dropna().reset_index(drop=True)
#如果编辑解释器不是jupyter-notebook或命令行，用print函数包裹下面的代码
df.head()

我们用pandas处理二维数据很方便，为何不直接用pd.read_csv读取文件？因为文件名是中文时read_csv读取会出错，而用_io.TextIOWrapper对象（既open读取返回的句柄）中转一下就不会出错。dropna意为删除NAN，就是缺失的数据所在的数据行，如果刚才没有补全数据，这条命令可以帮你删除，而reset_index是因为删除行后，pandas的行号不会改变，会对后面的调用产生一些影响，因此要重置行号。最后打印前五行。

此时数据还不能直接使用，首先地区的名称总不可作为判断的指标，因此要分离开来；其二，从数据本身考虑，GDP总量和GDP增速在作用上用重复，后面的数据也都是这样，因此要将后面的数据都二选一再进行判断。

indexs = ["GDP总量增速", "人口总量", "人均GDP增速", "地方财政收入总额", "固定资产投资", "社会消费品零售总额增速", "进出口总额","城镇居民人均可支配收入", "农村居民人均可支配收入"]

Positive = indexs
Negative = []

province = df['地区']
index = df[indexs]

显而易见，indexs为我们选择的指标，Positive、Negative分别为正、负向指标，而这里的指标统统为正向指标，因此Negative为空列表，province为Series对象，放着各个地区的名称。
现在数据方面处理完毕，我们分别保留着以下有用变量：

#变量名                 类型                    意义
province      #         series                  地区名称
Positive      #         list                    正向指标
Negative      #         list                    负向指标
index         #         dataframe               指标数据

（2）编写熵值法的类

首先，我们根据步骤，打好大致的框架：

class EmtropyMethod:
    def __init__(self, index, positive, negative, province):
        pass

    def uniform(self):
        pass

    def calc_probability(self):
        pass

    def calc_emtropy(self):
        pass

    def calc_emtropy_redundancy(self):
        pass

    def calc_Weight(self):
        pass

    def calc_score(self):
        pass

然后分析每个方法的功能来实现它：
__init__
__init__为初始化，我们只是为了将数据传入并检验格式是否正确。

def __init__(self, index, positive, negative, row_name):
    if len(index) != len(row_name):
        raise Exception('数据指标行数与行名称数不符')
    if sorted(index.columns) != sorted(positive+negative):
        raise Exception('正项指标加负向指标不等于数据指标的条目数')

    self.index = index.copy().astype('float64')
    self.positive = positive
    self.negative = negative
    self.row_name = row_name.copy()

传入的参数中，index为数据，dataframe类型；positive为正项指标，list类型；negative为负向指标，list类型；row_name为每行的含义，在这里指的就是地区。
如果地区数和指标的行数不同，或者如果正项指标+负向指标不等于全部指标，那显然就是出错了。（我这里比较指标的数量是将他们都排序，然后看他们是否相等）
如果都正确，那么将传递的值交给变量。
这里为了安全性考虑，DataFrame和Series类型都是引用赋值，为了安全性考虑，对象存的是传入参数的拷贝。
同时为了防止DataFrame类型的某些列是int类型，在后续计算出问题，将所有类型都转换为float64
uniform
归一化方法，将传递的数据归一化

def uniform(self):
    uniform_mat = self.index.copy()
    min_index = {column:min(uniform_mat[column]) for column in uniform_mat.columns}
    max_index = {column:max(uniform_mat[column]) for column in uniform_mat.columns}
    for i in range(len(uniform_mat)):
        for column in uniform_mat.columns:
            if column in self.negative:
                uniform_mat[column][i] = (uniform_mat[column][i] - min_index[column]) / (max_index[column] - min_index[column])
            else:
                uniform_mat[column][i] = (max_index[column] - uniform_mat[column][i]) / (max_index[column] - min_index[column])

    self.uniform_mat = uniform_mat
    return self.uniform_mat

对DataFrame中所有元素遍历，每次都判断是正向还是负向指标。按照公式计算。
可以暂时在中途调试（注意：这是调试代码）：

indexs = ["GDP总量增速", "人口总量", "人均GDP增速", "地方财政收入总额", "固定资产投资", "社会消费品零售总额增速", "进出口总额","城镇居民人均可支配收入", "农村居民人均可支配收入"]

Positive = indexs
Negative = []

province = df['地区']
index = df[indexs]

em = EmtropyMethod(index, Negative, Positive, province)
em.uniform()

归一化结果

calc_probability
对应原来计算指标比重的地方，既上面的 $p_{ij}$ ，注意这个公式中所用的 $x_{ij}$ 不再是原来数据中的，而是上一步归一化后的。

def calc_probability(self):
    try:
        p_mat = self.uniform_mat.copy()
    except AttributeError:
        raise Exception('你还没进行归一化处理，请先调用uniform方法')
    for column in p_mat.columns:
        sigma_x_1_n_j = sum(p_mat[column])
        p_mat[column] = p_mat[column].apply(lambda x_i_j: x_i_j / sigma_x_1_n_j if x_i_j / sigma_x_1_n_j != 0 else 1e-6)

    self.p_mat = p_mat
    return p_mat

这里会现将上一步的结果拷贝过来，如果uniform_mat不存在，显然是没进行归一化处理，就会向外抛出错误。
变量sigma_x_1_n_j是当前循环所在的 $j$ 列的所有指标加和，然后用apply将公式映射到当前列的每一个元素上。
可以看到匿名函数lambda表达式后面跟了一大串，意思是：传入 $x_{ij}$ ，将当前元素变为 $\frac{x_{ij}}{\sum_{i=1}{n}}$ ，如果 $x_ij$ 是0，那么就变成 $1 \times 10^{-6}$ 。
这是为什么？
首先看0是怎么来的，看看上一步归一化，如果（当前值是正向指标且是最小值）或（当前值是负向指标且是最大值），那么根据归一化公式就会等于 $0$ 。而到了这一步，分子为0，则算出来的结果比为0，也没什么问题，但下一步计算熵值的时候，会对 $p_{ij}$ 取对数 $ln$ ，如果这一步计算的哪个 $p_{ij}$ 是0，那么得出的值就会是无限大，在程序中体现就是NaN，因此要将值稍稍比1大一点。
同样，中途调试一下：

em = EmtropyMethod(index, Negative, Positive, province)
em.uniform()
em.calc_probability()

计算比重结果
calc_emtropy
计算熵值

def calc_emtropy(self):
    try:
        self.p_mat.head(0)
    except AttributeError:
        raise Exception('你还没计算比重，请先调用calc_probability方法')

    import numpy as np
    e_j  = -(1 / np.log(len(self.p_mat)+1)) * np.array([sum([pij*np.log(pij) for pij in self.p_mat[column]]) for column in self.p_mat.columns])
    ejs = pd.Series(e_j, index=self.p_mat.columns, name='指标的熵值')

    self.emtropy_series = ejs
    return self.emtropy_series

依旧是根据公式写代码，调试一下:

em = EmtropyMethod(index, Negative, Positive, province)
em.uniform()
em.calc_probability()
em.calc_emtropy()

计算信息熵结果
calc_emtropy_redundancy
计算信息熵冗余度

def calc_emtropy_redundancy(self):
    try:
        self.d_series = 1 - self.emtropy_series
        self.d_series.name = '信息熵冗余度'
    except AttributeError:
        raise Exception('你还没计算信息熵，请先调用calc_emtropy方法')

    return self.d_series

调试：

em = EmtropyMethod(index, Negative, Positive, province)
em.uniform()
em.calc_probability()
em.calc_emtropy()
em.calc_emtropy_redundancy()

信息熵冗余度

calc_Weight
计算权值
为了以后调度方便，可以一次调用获取，此方法内自动包含对之前方法的调用。

def calc_Weight(self):
    self.uniform()
    self.calc_probability()
    self.calc_emtropy()
    self.calc_emtropy_redundancy()
    self.Weight = self.d_series / sum(self.d_series)
    self.Weight.name = '权值'
    return self.Weight

调试：

em = EmtropyMethod(index, Negative, Positive, province)
em.calc_Weight()

权值

所有权值加和等于1。

calc_score
计算得分，同样，支持一次调用输出结果

def calc_score(self):
    self.calc_Weight()

    import numpy as np
    self.score = pd.Series(
        [np.dot(np.array(self.index[row:row+1])[0], np.array(self.Weight)) for row in range(len(self.index))],
        index=self.row_name, name='得分'
    ).sort_values(ascending=False)
    return self.score

将每一行的所有元素和权值逐个相乘再加和，是当前行的得分。
细节：dataframe取行我用的是切片，例如第0行就是self.index[0:1]，得出结果还是dataframe类型，我转化为numpy的array，因为dataframe是二维对象，因此得出结果是(1*n)的二维列表，因为我知道它只有一行，所以后面[0]转为一维列表，再与权值相乘求和，索引为地区，再进行降序排列。

（3）验证程序

数据上面已经处理好了，就是那四个重要数据，将它们传入类，然后调用计算得分方法就能出结果。

indexs = ["GDP总量增速", "人口总量", "人均GDP增速", "地方财政收入总额", "固定资产投资", "社会消费品零售总额增速", "进出口总额","城镇居民人均可支配收入", "农村居民人均可支配收入"]

Positive = indexs
Negative = []

province = df['地区']
index = df[indexs]

em = EmtropyMethod(index, Negative, Positive, province)
em.calc_score()

输出结果：

同时，因为封装成类，中间计算结果都得以保留，我们都可以用对象调用出来。

最终把整合好的类放出来：

class EmtropyMethod:
    def __init__(self, index, positive, negative, row_name):
        if len(index) != len(row_name):
            raise Exception('数据指标行数与行名称数不符')
        if sorted(index.columns) != sorted(positive+negative):
            raise Exception('正项指标加负向指标不等于数据指标的条目数')

        self.index = index.copy().astype('float64')
        self.positive = positive
        self.negative = negative
        self.row_name = row_name
        
    def uniform(self):
        uniform_mat = self.index.copy()
        min_index = {column:min(uniform_mat[column]) for column in uniform_mat.columns}
        max_index = {column:max(uniform_mat[column]) for column in uniform_mat.columns}
        for i in range(len(uniform_mat)):
            for column in uniform_mat.columns:
                if column in self.negative:
                    uniform_mat[column][i] = (uniform_mat[column][i] - min_index[column]) / (max_index[column] - min_index[column])
                else:
                    uniform_mat[column][i] = (max_index[column] - uniform_mat[column][i]) / (max_index[column] - min_index[column])

        self.uniform_mat = uniform_mat
        return self.uniform_mat
        
    def calc_probability(self):
        try:
            p_mat = self.uniform_mat.copy()
        except AttributeError:
            raise Exception('你还没进行归一化处理，请先调用uniform方法')
        for column in p_mat.columns:
            sigma_x_1_n_j = sum(p_mat[column])
            p_mat[column] = p_mat[column].apply(lambda x_i_j: x_i_j / sigma_x_1_n_j if x_i_j / sigma_x_1_n_j != 0 else 1e-6)

        self.p_mat = p_mat
        return p_mat
                 
    def calc_emtropy(self):
        try:
            self.p_mat.head(0)
        except AttributeError:
            raise Exception('你还没计算比重，请先调用calc_probability方法')

        import numpy as np
        e_j  = -(1 / np.log(len(self.p_mat)+1)) * np.array([sum([pij*np.log(pij) for pij in self.p_mat[column]]) for column in self.p_mat.columns])
        ejs = pd.Series(e_j, index=self.p_mat.columns, name='指标的熵值')

        self.emtropy_series = ejs
        return self.emtropy_series
        
    def calc_emtropy_redundancy(self):
        try:
            self.d_series = 1 - self.emtropy_series
            self.d_series.name = '信息熵冗余度'
        except AttributeError:
            raise Exception('你还没计算信息熵，请先调用calc_emtropy方法')

        return self.d_series
            
    def calc_Weight(self):
        self.uniform()
        self.calc_probability()
        self.calc_emtropy()
        self.calc_emtropy_redundancy()
        self.Weight = self.d_series / sum(self.d_series)
        self.Weight.name = '权值'
        return self.Weight
            
    def calc_score(self):
        self.calc_Weight()

        import numpy as np
        self.score = pd.Series(
            [np.dot(np.array(self.index[row:row+1])[0], np.array(self.Weight)) for row in range(len(self.index))],
            index=self.row_name, name='得分'
        ).sort_values(ascending=False)
        return self.score