Dagum基尼系数分解的MATLAB程序代码(更新)

很多同学在初次使用Dagum基尼系数法做研究时苦于找不到计算软件，手动计算又耗时费力。去年写论文时，参考Dagum(1997)论文中提供的算例，自己编写了一个MATLAB版的计算程序，在这里分享给大家。转载和使用请注明出处！

一、Dagum基尼系数分解简介

1、总体基尼系数

$G=\frac{1}{2 \overline{y} n^{2}}\left(\sum_{i=1}^{n} \sum_{r=1}^{n}\left|y_{i}-y_{r}\right|\right)=\sum_{j=1}^{k} \sum_{h=1}^{k} \sum_{i=1}^{n_{j}} \sum_{r=1}^{n_{h}}\left|y_{j i}-y_{h r}\right| / 2 n^{2} \overline{y}$

2、子群内部基尼系数

$\mathrm{G}_{\mathrm{jj}}=\frac{1}{2 \overline{y}_{j}}\left(\sum_{i=1}^{n_{j}} \sum_{r=1}^{n_{j}}\left|y_{j i}-y_{j r}\right|\right) / n_{j}^{2}$

3、子群内差异对总体基尼系数的贡献

$G_{w}=\sum_{i=1}^{n_{j}} G_{j j} P_{j} S_{j}$

4、子群h和子群j之间的基尼系数

$\mathrm{G}_{\mathrm{jh}}=\left(\sum_{i=1}^{n_{j}} \sum_{r=1}^{n_{h}}\left|y_{j i}-y_{h r}\right|\right) / n_{j} n_{h}\left(\overline{y}_{j}+\overline{y}_{h}\right)$

5、子群间净值差距对总体基尼系数的贡献度

$G_{n b}=\sum_{j=2}^{k} \sum_{h=1}^{j-1} \mathrm{G}_{\mathrm{jh}}\left(P_{j} S_{h}+P_{h} S_{j}\right) D_{j h}$

6、超变密度的贡献

$\mathrm{G}_{\mathrm{t}}=\sum_{j=2}^{k} \sum_{h=1}^{j-1} \mathrm{G}_{\mathrm{jh}}\left(P_{j} S_{h}+P_{h} S_{j}\right)\left(1-D_{j h}\right)$
其中：

$D_{j h}=\frac{d_{j h}-p_{j h}}{d_{j h}+p_{j h}}$

$d_{j h}=\int_{0}^{\infty} d F_{j}(y) \int_{0}^{y}(y-x) d F_{h}(x)$
$p_{j h}=\int_{0}^{\infty} d F_{h}(y) \int_{0}^{y}(y-x) d F_{j}(x)$

7、Dagum基尼系数分解

$\mathrm{G}=\mathrm{G}_{\mathrm{w}}+G_{n b}+G_{t}$

二、Dagum基尼系数代码

1、主函数文件（Dagum_Gini.m)

%Dagum基尼系数分解主函数文件
function [G_Total,G_W,G_sub,G_nb,G_jh,G_t,G_test]=Dagum_Gini(varargin)
%% 输入变量和输出变量说明
%
% 输入变量：假定有子群：A，B,C，...，首先，按各子群均值由大到小排序，
%         如:B的均值>A的均值>C的均值>...,则输入命令为：Dagum_Gini(B,A,C,...)
%         注：所有子群变量均以列向量形式输入。
%
% 输出变量：
%         G_Total：   总体基尼系数
%         G_W：       子群内差异贡献
%         G_sub：     子群内基尼系数
%         G_nb：      子群间差异贡献
%         G_jh：      子群间基尼系数
%         G_t：       超变密度贡献
%         G_test：    检验G_W+G_nb+G_t是否等于G_Total

%% 检查变量输入格式
n_sub_group=nargin;                     %输入变量数（子群数）
all_sub_group=varargin;                 %输入的所有子群（以元胞数组形式存储）

for i_group=1:n_sub_group
    [~,check_size]=size(all_sub_group{i_group});
    if check_size~=1
        error('输入格式错误：输入的各变量必须为列向量！')
    end
end

for i_group=1:n_sub_group-1
    if mean(all_sub_group{i_group})>=mean(all_sub_group{i_group+1})
        %disp('变量输入格式正确！')
    else
        error('输入格式错误：请按子群均值从大到小依次输入变量！')
    end
end

%% 计算整体和子群的样本数、均值，子群样本数占整体比重等
Total_group=cell2mat(all_sub_group');    %将子群合并成总体
n_T=length(Total_group);                 %计算总体样本数
m_T=mean(Total_group);                   %计算总体均值
        
% 预分配内存
n_sample_sub_group=zeros(n_sub_group,1); %存储各子群样本数量
m_sub_group=n_sample_sub_group;          %存储各子群均值
P=n_sample_sub_group;                    %存储各子群样本数占总体样本数比重
S=n_sample_sub_group;                    %存储各子群元素之和占总体元素之和的比重

for i_group=1:n_sub_group
    n_sample_sub_group(i_group)=length(all_sub_group{i_group});%b{i}是矩阵（列向量）
    m_sub_group(i_group)=mean(all_sub_group{i_group});
    P(i_group)=length(all_sub_group{i_group})/n_T;
    S(i_group)=length(all_sub_group{i_group})*mean(all_sub_group{i_group})/(n_T*m_T);
end

%% 载入基尼系数求解子函数
my_GINI=GINI_COMPUTE;

%% 计算各子群基尼系数
G_sub=zeros(n_sub_group,1);%预分配内存
for i_group=1:n_sub_group
    %循环计算各子群的基尼系数，并将结果存储在G_sub矩阵中
    G_sub(i_group)=my_GINI.GINI_1(all_sub_group{i_group},all_sub_group{i_group});
end

%% 使用for循环+矩阵上三角化计算G_jh、D_jh和P_S
%预分配内存
G_jh_T=zeros(n_sub_group,n_sub_group);
D_jh_T=zeros(n_sub_group,n_sub_group);
P_S_T=zeros(n_sub_group,n_sub_group);

for i_group=1:n_sub_group
    for j_group=1:n_sub_group
        %循环计算子群间基尼系数矩阵
        G_jh_T(i_group,j_group)=my_GINI.GINI_1(all_sub_group{i_group},all_sub_group{j_group});
        %循环计算子群间的相对影响矩阵
        D_jh_T(i_group,j_group)=my_GINI.D_jh(all_sub_group{i_group},all_sub_group{j_group});
        %循环计算权重矩阵
        P_S_T(i_group,j_group)=P(i_group)*S(j_group);
    end
end

% G_jh、D_jh和P_S的矩阵上三角化
G_jh=triu(G_jh_T,1);
D_jh=triu(D_jh_T,1); 
D_jh_1=1-D_jh;
D_jh_1=triu(D_jh_1,1);
P_S=triu(P_S_T,1)+triu(P_S_T',1); %（重要的一步）原矩阵与转置矩阵相加

%% 计算基尼系数分解结果
G_Total=my_GINI.GINI_1(Total_group,Total_group);   %计算总体基尼系数
G_nb=sum(sum(P_S.*D_jh.*G_jh));        %计算子群间差异贡献额
G_t=sum(sum(P_S.*D_jh_1.*G_jh));       %计算超变密度  
G_W=sum(P.*S.*G_sub);                  %计算组内基尼系数
G_test=G_W+G_nb+G_t;                   %检验总体基尼系数（G_Total）是否等于G_W+G_nb+G_t

%% 显示计算结果
disp('%%%%%%%%%%%%%%%%Dagum基尼系数分解结果%%%%%%%%%%%%%%%')
disp(' ')
disp(['总体基尼系数（G_T）: ',num2str(G_Total)])
disp(' ')
disp(['子群内差异贡献（G_W）: ',num2str(G_W)])
disp('------------各子群基尼系数------------')
for i_group=1:n_sub_group
    disp(['子群内基尼系数','（G_sub(',num2str(i_group),')):', num2str(G_sub(i_group))])
end
disp('-----------------------------------')
disp(' ')
disp([char('子群间差异贡献（G_nb）: ') num2str(G_nb)])

disp('------------子群间基尼系数------------')
for i_group=1:n_sub_group
    for j_group=1:n_sub_group
        if j_group>i_group
        disp(['子群间基尼系数','（G_jh_',num2str(i_group),'_',num2str(j_group),'):', num2str(G_jh(i_group,j_group))])
        end
    end
end 

disp('-----------------------------------')
disp(' ')
disp([char('超变密度贡献（G_t）: ') num2str(G_t)])
disp(' ')

disp('===========各部分贡献率(%)===========')
disp(['子群内差异贡献率:',num2str(G_W/G_Total*100)])
disp(['子群间差异贡献率:',num2str(G_nb/G_Total*100)])
disp(['超变密度贡献率:',num2str(G_t/G_Total*100)])
disp('===================================')
disp(' ')

disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')

end

2、子函数文件(GINI_COMPUTE.m)

%% 计算整体、子群内、子群间基尼系数和子群间相对影响

function GINI_MAIN=GINI_COMPUTE

%将一组功能独立的函数写在同一个函数文件中
%主函数GINI_MAIN(输出)=GINI_COMPUTE(输入)
%注1：运行该函数时可使用：my_GINI=GINI_COMPUTE
%注2：对应的子函数输入形式为：my_GINI.GINI_1

GINI_MAIN.GINI_1=@GINI_1;%计算整体、子群内和子群间间基尼系数
GINI_MAIN.D_jh=@D_jh;    %计算子群间相对影响（D_ij）
end

%% 子函数1：计算整体、子群内和子群间基尼系数

function GINI1=GINI_1(xx,yy)
n1=length(xx);
n2=length(yy);
m1=mean(xx);
m2=mean(yy);
D_ij_2=zeros(n1,n2);

for ii=1:n1
    for jj=1:n2
        D_ij_2(ii,jj)=xx(ii)-yy(jj);%列向量元素两两相减的差值矩阵
    end
end

sum_d2=sum(sum(abs(D_ij_2))); %对矩阵所有元素先取绝对值，然后求和 
GINI1=sum_d2/(n1*n2*(m1+m2)); %计算总体、组内和组间基尼系数
end

%% 子函数2:计算子群j和子群h变量的相对影响D_jh

function D_jh1=D_jh(xx,yy)
n1=length(xx);
n2=length(yy);
D_ij_3=zeros(n1,n2);

for ii=1:n1
    for jj=1:n2
        D_ij_3(ii,jj)=xx(ii)-yy(jj);
    end
end
%从矩阵中找出数值大于0的元素 
E11=sort(D_ij_3(find(D_ij_3>0)));
%子群j和子群h间差值(算术平均)
d_jh=sum(E11)/(n1*n2);   

%从矩阵中找出数值小于0的元素并对他们取绝对值
E12=abs(sort(D_ij_3(find(D_ij_3<0))));
%超变一阶矩（算术平均）
p_jh=sum(E12)/(n1*n2);

%j区域和h区域变量的相对影响
D_jh1=(d_jh-p_jh)/(d_jh+p_jh);
%注1：若子群j的样本均值和子群h的样本均值相同，则D_jh1=0
%注2：若子群j的所有元素均大于子群子群h的元素，则D_jh1=1
%注3：其他情况下，0<D_jh1<1
end

三、示例：

将一个整体区域划分为三个子区域A、B、C，计算该整体区域人均GDP的Dagum基尼系数分解。

第一步：保存函数文件

点击MATLAB窗口左上角的新建脚本按钮，复制上面的主函数文件（Dagum_Gini.m）和子函数文件（GINI_COMPUTE.m）代码，分别粘贴到两个MATLAB新建脚本中，然后点击保存。（注：保存在同一个文件夹中。直接点保存，不用重新命名）。

第二步：运行

运行方式1：在MATLAB的命令行窗口中输入下列示例代码，按回车键运行。
运行方式2：将下列示例代码保存到一个新建脚本文件中，对该脚本自由命名，例如Dagum_GINI_example.m，然后在命令行窗口输入Dagum_GINI_example，按回车键运行。建议使用运行方式2。

示例代码：

%子群A人均GDP(列向量)
A=[62905,45056,64746,43965,66318,30137,19014,30432,41683,26687,52590,50474,38357,32047,43442,39605,29063,30105,27566,10814,42063,25933,13077,40116,10949,10589]';
%子群B人均GDP（列向量）
B=[35500,19409,20263,5891,8828,10099,20355,8001,14157,14478,33711,20387,19171,17799,10020,15901,12232,11493,30460,12590,16899,20156,25013,9478,17241,8228,9570,8559]';
%子群C人均GDP（列向量）
C=[16629,30006,394,403,648,673,304,374,452,508,343,529,338,176,373,510]';

%经计算A的均值>B的均值>C的均值，所以Dagum基尼系数分解的输入指令如下：
Dagum_Gini(A,B,C);

注：当数据量较大时，手动输入比较麻烦。可从excel中导入数据，并进行变量命名。或在MATLAB的主页栏中找到并点击新建变量按钮，然后将数据粘贴到数据框中并命名。具体操作读者可自行百度。

输出结果：

输出结果1.jpg

注1：G_sub(1)为子群A基尼系数，G_sub(2)为子群B基尼系数。按变量输入顺序逐一对应。
注2：G_jh_1_2为A和B的子群间基尼系数，G_jh_1_3为A和C的子群间基尼系数。同样按输入变量顺序逐一对应。
注3：上述计算结果仅为单一年份的结果，要计算不同年份的基尼系数分解结果，需根据每年的子群数据逐一使用该程序计算。
注4：该程序支持3个以上的子群划分。例如存在A、B、C、D四个子群，且有B的均值>A的均值>D的均值>C的均值，则输入Dagum_Gini(B,A,D,C)即可。

四、特殊情形

1、子群的分布完全不重合：

例如，有A、B、C三个子群，子群C的所有元素均大于子群B的元素，子群B的所有元素均大于子群A的元素。此时，子群间的相对影响 $D_{jh}=1$ ，超变密度 $G_{t}=0$ 。

测试代码：

A=[1:0.5:10]';
B=[11:0.8:20]';
C=[21:0.6:30]';
%C的均值>B的均值>A的均值
Dagum_Gini(C,B,A);

输出结果：

输出结果2.jpg

2、子群的均值相同：

例如，有A、B、C三个子群，存在 $\overline A=\overline B=\overline C$ ，则，子群间相对影响 $D_{jh}=0$ ，组间差异贡献 $G_{nb}=0$ 。

测试代码：

a=magic(10);%生成魔方矩阵（每一列均值都相等）
A=a(:,1);
B=a(:,2);
C=a(:,3);

%A的均值=B的均值=C的均值
Dagum_Gini(A,B,C);

输出结果

输出结果3.jpg