Python で
パレート分布とは
統計学入門 - 東京大学出版会 の6.12 パレート分布 には
対数正規分布が
全集団の 所得分布であるのに たいし、 パレート分布は 高額所得者の 所得分布である。 $x_{0}$ 以上の 所得の 確率は、 対数正規分布よりも この 分布が よく あてはまる。 パレート分布 Pareto distributionはある 定数$x_{0}$ 以上で 存在し、 そこから 急に 減少する 図6.16 のような 密度関数を もつ。
図は引用できませんが、
また、
- イタリアの
経済学者ヴィルフレド・パレートが 所得の 分布を モデリングする 分布と して 提唱した 連続型の 確率分布 - 離散型は
ゼータ分布（ジップ分布）である。
ジップの
Wikipedia に
以下、
Python での
調べた
パレート分布
一般化パレート分布
generalized Pareto distributions, GPD
scipy に一般化パレート分布 generalized Pareto distributions, GPDの計算が できる モジュールが ありました。 ジップの
法則
scipy、numpy で 計算モジュールが ありました。
各サンプルプログラムを
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()
分布の
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()
モードを
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()
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>]
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>]
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>]
c はc = 0 の
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>]
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>]
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>]
numpy.random.zipf
続いて、numpy.random.zipf で
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()
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()
参考
以下、
scipy.stats.genpareto、scipy.stats.zipf Web解析用途で
以上です。
コメント