Python でパレート分布、ジップの法則を計算する方法を調べてみました。 結果を記載します。

パレート分布とは

統計学入門 - 東京大学出版会 の 6.12 パレート分布 には以下の通り、記載されています。

対数正規分布が全集団の所得分布であるのにたいし、パレート分布は高額所得者の所得分布である。$x_{0}$ 以上の所得の確率は、対数正規分布よりもこの分布がよくあてはまる。パレート分布Pareto distribution はある定数$x_{0}$ 以上で存在し、そこから急に減少する図6.16 のような密度関数をもつ。

図は引用できませんが、パレート分布 - Wikipedia に似たような図の記載があります。
また、Wikipedia にはパレート分布について以下記載があります。

  • イタリアの経済学者ヴィルフレド・パレートが所得の分布をモデリングする分布として提唱した連続型の確率分布
  • 離散型はゼータ分布（ジップ分布）である。

ジップの法則について

Wikipedia に ジップの法則へのリンクがあり、気になったので調べてみました。
以下、目を通した記事になります。

Python でのパレート分布、ジップの法則の計算方法

調べた限り、以下で計算ができるようです。

各サンプルプログラムを実行してみます。

numpy.random.pareto

np.random.pareto でパレート分布に従う配列を取得できます。 

import numpy as np
import matplotlib.pyplot as plt
a, m = 6., 2.  # 分布の幅と、モードを指定   
s = (np.random.pareto(a, 1000) + 1) * m
#  ヒストグラムを描画
count, bins, _ = plt.hist(s, 100, density=True)

# 線グラフの描画   
fit = a * m ** a / bins ** (a + 1)
plt.plot(bins, max(count) * fit / max(fit), linewidth=2, color='r')
plt.show()

20190327_output_1_0.png - Google ドライブ

分布の幅(shape 、scale と言われている値)を調整すると、グラフの最大値に影響し、大きくすると最大値が大きくなります。 

import numpy as np
import matplotlib.pyplot as plt
a, m = 12., 2.  # 分布の幅と、モードを指定   
s = (np.random.pareto(a, 1000) + 1) * m
#  ヒストグラムを描画
count, bins, _ = plt.hist(s, 100, density=True)
# 線グラフの描画   
fit = a * m ** a / bins ** (a + 1)
plt.plot(bins, max(count) * fit / max(fit), linewidth=2, color='r')
plt.show()

20190327_output_3_0.png - Google ドライブ

モードを調整すると、x軸の開始位置が移動します。 

import numpy as np
import matplotlib.pyplot as plt
a, m = 6., 4.  # 分布の幅と、モードを指定   
s = (np.random.pareto(a, 1000) + 1) * m
#  ヒストグラムを描画
count, bins, _ = plt.hist(s, 100, density=True)

# 線グラフの描画   
fit = a * m ** a / bins ** (a + 1)
plt.plot(bins, max(count) * fit / max(fit), linewidth=2, color='r')
plt.show()

20190327_output_5_0.png - Google ドライブ

scipy.stats.pareto

scipy.stats.pareto を使います。 

from scipy.stats import pareto
b = 2.62
# モーメントの計算、出力    
mean, var, skew, kurt = pareto.stats(b, moments='mvsk')
print(mean, var, skew, kurt)    

1.6172839506172838 1.6101990746886534 nan nan

パレート分布 - Wikipedia に記載がありますが、b <= 3 のため、歪度、尖度 は 末定義になります。 

from scipy.stats import pareto
b = 4.5
# モーメントの計算、出力    
mean, var, skew, kurt = pareto.stats(b, moments='mvsk')
print(mean, var, skew, kurt)    

1.2857142857142858 0.1469387755102041 5.465943944999486 146.44444444444446

b > 4 を設定すると、歪度、尖度 の計算結果が返ります。 

import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)
b = 2.62
# 確率密度関数（PDF）、パーセント点関数 (PPF) 
# パレート分布を描画する     
x = np.linspace(pareto.ppf(0.01, b),pareto.ppf(0.99, b), 100)    
ax.plot(x, pareto.pdf(x, b),'r-', lw=5, alpha=0.6, label='pareto pdf')     
rv = pareto(b)    
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

[<matplotlib.lines.Line2D at 0x12dd5d8d0>]

20190327_output_11_1.png - Google ドライブ

ppf、pdf に loc と scale を指定すると、グラフ切片、グラフのスケールが変わります。 

import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)
b = 2.62
# グラフのx軸がずれる
loc = 30
# scale で結果が割られるので、0.5を指定すると、パレート分布y軸の値が2倍になる
scale = 0.5
# 確率密度関数（PDF）、パーセント点関数 (PPF) 
# パレート分布を描画する     
x = np.linspace(pareto.ppf(0.01, b, loc, scale),pareto.ppf(0.99, b, loc, scale), 100)    
ax.plot(x, pareto.pdf(x, b, loc, scale),'r-', lw=5, alpha=0.6, label='pareto pdf')     

# parete関数を使用する場合は、後続の pdf ppf 関数には、loc、scale は指定不要
rv = pareto(b, loc, scale)    
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

[<matplotlib.lines.Line2D at 0x12dfccdd8>]

20190327_output_13_1.png - Google ドライブ

scipy.stats.genpareto

scipy.stats.genpareto では 一般化パレート分布 の計算が行えます。使い方はscipy.stats.pareto と同様です。 

from scipy.stats import genpareto
c = 0.1
# モーメントの計算、出力    
mean, var, skew, kurt = genpareto.stats(c, moments='mvsk')
print(mean, var, skew, kurt) 

1.1111111111111116 1.543209876543199 2.8110568859996445 14.828571428571088

import matplotlib.pyplot as plt
from scipy.stats import genpareto
c = 0.1
fig, ax = plt.subplots(1, 1)
x = np.linspace(genpareto.ppf(0.01, c),
                genpareto.ppf(0.99, c), 100)
ax.plot(x, genpareto.pdf(x, c),'r-', lw=5, alpha=0.6, label='genpareto pdf')
rv = genpareto(c)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

[<matplotlib.lines.Line2D at 0x12e0836d8>]

20190327_output_16_1.png - Google ドライブ

c はパレート指数を表します。c = 0 の場合は指数分布、c = -1 の場合は一様分布になり、 c < 0 の場合は タイプ2 のパレート分布になります。 

import matplotlib.pyplot as plt
from scipy.stats import genpareto
# c = 0 の場合は指数分布となる
c = 0
fig, ax = plt.subplots(1, 1)
x = np.linspace(genpareto.ppf(0.01, c),
                genpareto.ppf(0.99, c), 100)
ax.plot(x, genpareto.pdf(x, c),'r-', lw=5, alpha=0.6, label='genpareto pdf')
rv = genpareto(c)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

[<matplotlib.lines.Line2D at 0x12e5afc88>]

20190327_output_18_1.png - Google ドライブ 

import matplotlib.pyplot as plt
from scipy.stats import genpareto
# c = -1 の場合は一様分布になる
c = -1
fig, ax = plt.subplots(1, 1)
x = np.linspace(genpareto.ppf(0.01, c),
                genpareto.ppf(0.99, c), 100)
ax.plot(x, genpareto.pdf(x, c),'r-', lw=5, alpha=0.6, label='genpareto pdf')
rv = genpareto(c)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

[<matplotlib.lines.Line2D at 0x12e77c940>]

20190327_output_19_1.png - Google ドライブ 

import matplotlib.pyplot as plt
from scipy.stats import genpareto
# c = -0.80 の場合はタイプ2 のパレート分布になる。c > 0 の場合とは形が逆になる   
c = -0.80    
fig, ax = plt.subplots(1, 1)
x = np.linspace(genpareto.ppf(0.01, c),
                genpareto.ppf(0.99, c), 100)
ax.plot(x, genpareto.pdf(x, c),'r-', lw=5, alpha=0.6, label='genpareto pdf')
rv = genpareto(c)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

[<matplotlib.lines.Line2D at 0x12eedb898>]

20190327_output_20_1.png - Google ドライブ

numpy.random.zipf

続いて、numpy.random.zipf で zip の法則(分布)の計算をします。 

import numpy as np
import matplotlib.pyplot as plt
from scipy import special
a = 2. # parameter
# zip の法則に従う配列を生成
s = np.random.zipf(a, 1000)
# ヒストグラムを描画
count, bins, ignored = plt.hist(s[s<50], 50, density=True)
x = np.arange(1., 50.)
# special.zetac はゼータ分布の計算に用いる関数    
y = x**(-a) / special.zetac(a)
plt.plot(x, y/max(y), linewidth=2, color='r')
plt.show()

20190327_output_23_0.png - Google ドライブ

scipy.stats.zipf

scipy.stats.zipf の使い方は、scipy.stats.pareto 等と似ています。 

from scipy.stats import zipf
a = 6.5
# モーメントの計算、出力    
mean, var, skew, kurt = zipf.stats(a, moments='mvsk')
print(mean, var, skew, kurt)

1.0130420979439105 0.015940729460401704 12.561822558088188 279.4461024601007

from scipy.stats import zipf
import matplotlib.pyplot as plt
a = 3
fig, ax = plt.subplots(1, 1)
# zipfの 1%点、99%点を取得   
x = np.arange(zipf.ppf(0.01, a),
               zipf.ppf(0.99, a))
# 確率質量関数(Probability Mass Function, PMF)をplotする
ax.plot(x, zipf.pmf(x, a), 'bo', ms=8, label='zipf pmf')
# 縦棒を引く
ax.vlines(x, 0, zipf.pmf(x, a), colors='b', lw=5, alpha=0.5)
rv = zipf(a)
# 縦棒を引く
ax.vlines(x, 0, rv.pmf(x), colors='k', linestyles='-', lw=1, label='frozen pmf')
ax.legend(loc='best', frameon=False)
plt.show()

20190327_output_27_0.png - Google ドライブ

参考

以下、記事作成中に参考にしました。

scipy.stats.genpareto、scipy.stats.zipf Web解析用途でよく使いそうに思いました。
以上です。

コメント

カテゴリー

最近の投稿