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

2023年9月11日2023年10月16日

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

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

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

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

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

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

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

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

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

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

1. ジニ指数を求める３つの計算式
2. 拡張ジニ指数の計算式をPythonに実装する
3. 高ジニ指数の業界を見てみる

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

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

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

すべての標本が正の数の場合（＝普通のジニ指数）
標本に１つでも負の数が含まれる場合
標本合計が負の数の場合

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

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

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

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

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

この図においてAとBの面積がわかれば、

ジニ指数＝A／(A+B)

で計算できます。

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

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

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

ジニ指数＝(A+C)／(A+B+C)

標本合計が負の数の場合

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

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

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

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

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

分析データをImportする

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

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

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

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

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

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

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

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

ローレンツ曲線を描く

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

（外食＆娯楽サービスの業界IDは23）

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")

A～Eの面積を求める

これからA～Eの面積を求めるわけですが、まずはA～Cから求めます。

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

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

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

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

A＝0.5－B

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

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

このようにしてx座標のリスト、y座標のリストともに２つに分けると、前半の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

次にDとEの面積を求めます。

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

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

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

以上を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

拡張ジニ指数を求める

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

拡張ジニ指数

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

＝０.７９

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

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

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

全業界の業界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)

# 上記リストに含まれる要素数＋１の数だけxとして０～１の要素を等間隔に作る
    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