標本に負の数が含まれる場合のジニ指数をPythonで計算【標本合計が負でもOK!】

2023年10月16日

佐々木朗希復帰戦 on 2023/09/10

標本(サンプル)の格差を表すのにジニ指数は用いられますが、全標本が正の数であることを前提としています。

例えば、標本が各社の売上高である場合、すべての標本は正の数なのでこの前提を満たします。

ところが標本が各社の営業利益である場合、負の数が含まれる可能性があります。

このような時には、まともに計算するとジニ指数は変な値になってしまいます。

0~1に収まるはずのジニ指数が1を超えたり、マイナスになったりするのです。

そのような例は【物流業界内の格差は大きいのか?】ジニ指数で他業界と比較してみたの最後の章で紹介していますので、興味のある方はご覧ください。

 

でも世の中には同じような問題意識を持った人がいて、このようなジニ指数の欠点を解消してくれています。

管理人が参考にしたのはこちらの論文です。

標本合計が負の場合へ拡張されたジニ係数の評価|新潟大学

 

本記事ではこの論文の要点を紹介し、Pythonに実装して負の数が含まれる営業利益率のジニ指数を計算してみます。

 

◆仕事や勉強の息抜きに。。。

ジニ指数を求める3つの計算式

まずは先の論文の要点だけを簡単に整理しておきます。

下記の3つの計算方法があります。

  • すべての標本が正の数の場合(=普通のジニ指数)
  • 標本に1つでも負の数が含まれる場合
  • 標本合計が負の数の場合

 

下に行くに従い汎用性が増し、1番目と2番目の式は3番目の式の特殊ケースということになります。

 

すべての標本が正の数の場合

これは普通のジニ指数の計算方法です。

【ジニ指数をExcelとPythonで計算する】テスト結果を例にわかりやすく解説に載せている方法と同じで次のように計算できます。

標本合計が負の場合へ拡張されたジニ係数の評価|新潟大学より抜粋

 

この図においてABの面積がわかれば、

ジニ指数=A(A+B)

で計算できます。

 

標本に1つでも負の数が含まれる場合

標本の中に少しだけ負の数が含まれている場合、もう少し正確にいうと、標本の中に1つ以上の負の数が含まれているけれども標本合計が正の数の場合は、下図におけるACの面積がわかれば下式でジニ指数が計算できます。

標本合計が負の場合へ拡張されたジニ係数の評価|新潟大学より抜粋

 

ジニ指数=(A+C)(A+B+C)

 

標本合計が負の数の場合

標本の中に負の数が多く含まれ標本合計が負の数になってしまう場合は、下図におけるAEの面積がわかれば下式でジニ指数が計算できます。

標本合計が負の場合へ拡張されたジニ係数の評価|新潟大学より抜粋

 

ジニ指数=(A+C)(A+B+C+D+E)

 

この3番目の式を拡張ジニ指数の計算式と呼ぶことにします。

 

拡張ジニ指数の計算式をPythonに実装する

分析データをImportする

Investing.comで上場企業の過去4年分の売上高/粗利益/営業利益/当期利益が取得できます。

また業界も52に分類されているので、業界ごとの売上高の分布もわかります。

これらのデータはInvestpyというライブラリーを使えばPythonで自動取得できるので、取得後、1年分だけのデータを残して下記のようなデータフレームに整理しました。

 

最終行に営業利益÷売上高で計算した営業利益率を入れています。

このデータを使って各業界ごとの営業利益率のジニ指数を求めていきます。

営業利益率は赤字の会社の場合はマイナスになるので、負の数の標本が含まれる可能性があり、また赤字の会社が多ければ標本合計も負の数になる可能性があります。

従って、拡張ジニ指数の計算式(3番目の式)をPythonに実装します。

外食&娯楽サービス”業界はコロナの影響で赤字の会社が多かったため、まずはこの業界を例に計算してみましょう。

 

ローレンツ曲線を描く

ローレンツ曲線を描いてみると、次のようにマイナス領域に大きく食い込んだ曲線になります。

(外食&娯楽サービスの業界ID23

import matplotlib.pyplot as plt
import numpy as np

rev_list = df[df['ID_Industry'] == '23']['OP Rate1'].values
rev_list = np.sort(rev_list)

x = np.linspace(0, 1, len(rev_list)+1)
x = x.tolist()

rel = rev_list / np.sum(rev_list)
rel_cum = rel.cumsum()
rel_cum = rel_cum.tolist()
rel_cum = [0] + rel_cum

plt.plot(x, rel_cum)
plt.xlabel("Company")
plt.ylabel("Cumulative Revenue")

 

AEの面積を求める

これからAEの面積を求めるわけですが、まずはACから求めます。

一見難しそうですが、scipyライブラリのintegrate.trapz関数を使えばローレンツ曲線の積分(ローレンツ曲線とx軸で囲まれた部分の面積)を計算できるので、次のように求めることができます。

ローレンツ曲線とx軸との交点を(k,0)とすると、

C=ローレンツ曲線を0からkまで積分した値×(-1)

B=ローレンツ曲線をkから1まで積分した値

A0.5B

 

 

このようにローレンツ曲線のy座標のリストを見ていって、値が負でなくなるのは何番目かを調べます。

そしてx座標のリストについても同じ数だけ要素を取り出します。

このようにしてx座標のリスト、y座標のリストともに2つに分けると、前半のx座標-y座標のリストの組からCの面積が、後半のそれからはBの面積が求まります。

以上をPythonでは次のようにコーディングできます。

rel_cum1, rel_cum2 = np.hsplit(rel_cum, [np.count_nonzero(rel_cum < 0) + 1])
x1, x2 = np.hsplit(x, [np.count_nonzero(rel_cum <  0) + 1])
rel_cum1 = rel_cum1.tolist()
area1 = integrate.trapz(rel_cum1, x1)
rel_cum2 = rel_cum2.tolist()
area2 = integrate.trapz(rel_cum2, x2)
C = - area1
B = area2
A = 0.5 - area2

 

次にDEの面積を求めます。

今度はローレンツ曲線のy座標のリストを見ながら極小点が何番目にあるかを調べ、同じ数だけの要素をx座標のリストからも取り出します。

ここで求めたいのはローレンツ曲線とその極小点で接する直線x=ρで囲まれる面積なので、ローレンツ曲線をρだけ上に平行移動します。

こうすることによって、求めたい面積DEがローレンス曲線の積分値になるからです。

 

以上をPythonではつぎのようにコーディングできます。

rel_cum3, rel_cum4 = np.hsplit(rel_cum, [np.where(rel_cum == rel_cum.min())[0][0] + 1])
rel_cum3 = rel_cum3 - rel_cum.min()
rel_cum4 = rel_cum4 - rel_cum.min()
x3, x4 = np.hsplit(x, [np.where(rel_cum == rel_cum.min())[0][0] + 1])
rel_cum3 = rel_cum3.tolist()
area3 = integrate.trapz(rel_cum3, x3)
rel_cum4 = rel_cum4.tolist()
area4 = integrate.trapz(rel_cum4, x4)
D = area3
E = area4 – B

 

拡張ジニ指数を求める

これでAEの面積が求まったので、あとはこれらを先程の拡張ジニ指数の公式に入れるだけです。

 

拡張ジニ指数

(A+C)(A+B+C+D+E)

.79

 

外食&娯楽サービス業界の営業利益率のジニ指数は.79になりました。

 

各業界のジニ指数を求める

それでは他の全業界のジニ指数を求めてみましょう。

全業界の業界IDを入れたファイルを作り、for文で繰り返せば良いだけです。

# 業界IDのリストの各要素を文字列にする
ids0 = df_industry['ID_Industry'].tolist()
ids = []
for id in ids0:
    a = str(id)
    ids.append(a)

# 業界ごとのジニ指数と会社数を入れる空データフレームを作る
df_gini = pd.DataFrame(columns=['gini', 'cnt'])
# 各業界の会社数iをゼロに初期化する
i = 0
for id in ids:
# 指定した業界の会社ごとの営業利益率のリストを作り、ソートする
    rev_list = df[df['ID_Industry'] == id]['OP Rate1'].values
    rev_list = np.sort(rev_list)

# 上記リストに含まれる要素数+1の数だけxとして0~1の要素を等間隔に作る
    x = np.linspace(0, 1, len(rev_list)+1)

# 営業利益率のリストを累積相対数に変換する
    rel = rev_list / np.sum(rev_list)
    rel.sort()
    rel_cum = rel.cumsum()
    rel_cum = rel_cum.tolist()
    rel_cum = [0] + rel_cum
    rel_cum = np.array(rel_cum)

# A~Eの面積を求める
    rel_cum1, rel_cum2 = np.hsplit(rel_cum, [np.count_nonzero(rel_cum < 0) + 1])
    rel_cum3, rel_cum4 = np.hsplit(rel_cum, [np.where(rel_cum == rel_cum.min())[0][0] + 1])
    rel_cum3 = rel_cum3 - rel_cum.min()
    rel_cum4 = rel_cum4 - rel_cum.min()
    x1, x2 = np.hsplit(x, [np.count_nonzero(rel_cum <  0) + 1])
    x3, x4 = np.hsplit(x, [np.where(rel_cum == rel_cum.min())[0][0] + 1])

    rel_cum1 = rel_cum1.tolist()
    area1 = integrate.trapz(rel_cum1, x1)
    rel_cum2 = rel_cum2.tolist()
    area2 = integrate.trapz(rel_cum2, x2)
    rel_cum3 = rel_cum3.tolist()
    area3 = integrate.trapz(rel_cum3, x3)
    rel_cum4 = rel_cum4.tolist()
    area4 = integrate.trapz(rel_cum4, x4)

    C = - area1
    B = area2
    A = 0.5 - area2
    D = area3
    E = area4 - B

# 拡張ジニ指数を求める
    gini_ex = (A + C) / (A + B + C + D + E)

# データフレームに結果を書き込む
    df_gini.loc[i, 'gini'] = gini_ex
    df_gini.loc[i, 'cnt'] = len(rev_list)

    i = i + 1

# データフレームに業界IDを加える
df_gini = df_gini.reset_index()
df_gini = df_gini.rename({'index': 'ID_Industry'}, axis='columns')
df_gini['ID_Industry'] = df_gini['ID_Industry'].astype(str)
df_gini.columns = ['ID_Industry', 'gini', 'cnt']

# データフレームに業界名と順位を加える
df_gini = pd.merge(df_gini, df, on = 'ID_Industry', how = 'left')
df_gini = df_gini[['ID_Industry', 'Industry', 'gini', 'cnt']].drop_duplicates().reset_index().drop('index', axis=1)
df_gini['gini_rank'] = df_gini['gini'].rank(ascending = False).astype(int)
df_gini = df_gini.sort_values('gini_rank', ascending = True)
df_gini

 

物流業界52業界の中で45番目に低い.36でした

外食&娯楽サービス業界に比べて、営業利益率の格差が遥かに小さい業界ということができます。

ちなみにローレンツ曲線を描いてみると、こんな感じでy=xに近い曲線になっています。

ちなみに縦軸を累積相対営業利益率ではなく、営業利益率でプロットしてみると次のようなグラフになります。

このように見ると、上位10%くらいの会社は飛び抜けて高い営業利益率を叩き出しているように見えますが、これらの会社はすべてコロナ特需のあった海運会社です。

これらの会社を除くと、営業利益率010%の間にほぼ均等に収まっていることがわかります。

 

高ジニ指数の業界を見てみる

ジニ指数上位の業界のうち、会社数が多い“外食&娯楽サービス”を見てみましょう。

ジニ指数は先ほど計算したように.79で全体の6位です。

ローレンツ曲線は次のようになります。

 

かなり特殊なローレンス曲線になっています。

大きくマイナスにせり出しているのは、それだけ赤字の会社の多いということです。

先ほどの物流業界と同じように、縦軸を累積相対営業利益率ではなく、営業利益率でプロットしてみましょう。

 

一番業績の悪かった会社は営業利益率が約-80%で、一番業績の良かった会社は約40%と大きな開きがあり、ジニ指数が大きいのも頷けます。

また赤字の会社グループ、普通の会社のグループ、高収益の会社のグループの3つにはっきり分かれていることもわかります。

業績の悪かった10社は次のようになっています。

 

逆に業績の良かった10社は次の通りです。(下から順に良かった会社)

 

ほとんどの会社の営業利益率が0~10%に収まっている物流業界は、コロナ禍でも安定した業界だったと言えるでしょう。