様々な分布の定義を確認して、実際のデータが分布に従っているか確認する方法を見ていく。連続分布に従って生成される実際の分布は離散的なので\(\rm PDF\)は使い勝手が悪い。\(\rm CDF\)について考える。
指数分布
\[ \begin{align} \rm (pd)f &= &\lambda e^{-\lambda x}\\ \rm (CD)F &= &1 - e^{-\lambda x}\\ mean(F) &= &1/\lambda\\ median(F) &= &\log(2)/\lambda \end{align} \]
である。
このCDFを変形してやると \[ log(CCDF) = log(1-CDF) = -\lambda x \] なのでCCDFをy軸をlogスケールにしてプロットすれば直線になるはずである。演習問題4-1に合わせて、平均32.6、44個の標本も含めてプロットすると
from __future__ import division import numpy as np import matplotlib.pyplot as plt fig = plt.figure() ax = fig.add_subplot(111) ax.set_ylabel('ccdf') ax.set_yscale('log') for i, size in enumerate((44, 100, 1000, 100000)): pmf_5, bins = np.histogram(np.random.exponential(scale=32.6, size=size), bins=50, range=(0,250), density=True) ccdf = 1 - np.cumsum(pmf_5*5) bins = (bins + 2.5)[:-1] ax.plot(bins, ccdf, color=('k','r','g','b')[i], label=[str(j) for j in range(4)]) ax.legend([str(j)+'sample' for j in (44, 100, 1000, 100000)], loc='best') plt.show()
となる。
パレート分布
wikipediaによる式はかなりややこしい。簡略化されたものがテキストにはあって、 \[ \begin{align} \rm pdf &= &\dfrac{\alpha}{x_{m}}\biggl(\dfrac{x}{x_{m}}\biggr)^{-\alpha-1}\\ \rm cdf &= &1 - \biggl(\dfrac{x}{x_{m}}\biggr)^{-\alpha}\\ \end{align} \] ただし、\(x\)の最小値は\(x_{m}\)である。cdfを変形して、 \[ \log(\rm ccdf) = -\alpha(\log{x} - log{x_{m}}) \]
パレート分布は年収の分布を表すのに使えるらしい。京大卒の平均年収は677万円、先進国のパレート係数は1.5くらい(出典が分からない)というのが見えた。実際にシミュレーションしてみる。手元で計算してみたところ、\(x_{m}=400\)くらいになり、あまりにも現実にそぐわない気がするので撤退。
現実は以下
参考: 全300校 出身大学別年収データ | 20代の”はたらき”データベース『キャリアコンパス』- powered by DODA -
パレート分布はスケールフリー性があるらしいのでプロットしてみる。
from __future__ import division import numpy as np from matplotlib import pylab L = np.sort(np.random.pareto(1.5, size=100000)) pylab.scatter(L, L/np.sum(L)) pylab.show()
以下上から下へ順に拡大していっている。
正規分布
\(\rm erf(x) = 2/\sqrt{\pi} \displaystyle\int_{0}^{x} e^{-t^{2}}\)
対数正規分布
\[ \rm cdf_{\rm lognormal} \equiv \rm cdf(log(x)) \]
対数正規分布の仕組みによると、むしろ日本の年収は対数正規分布に従うらしい。
from __future__ import division import numpy as np from matplotlib import pylab pylab.hist(np.random.lognormal(mean=2,sigma=0.1,size=10000),bins=100, range=(0,20)) pylab.show() salary = 300 + 10 * np.random.random_sample(1000) pylab.hist(salary, range=(250,500), bins=500) pylab.show() for i in range(100): if i % 20 == 0: pylab.hist(salary, range=(250,1000), bins=500) pylab.show() up = [] for j in 1000*np.random.random_sample(20): j = int(j) if j not in up: salary[j] = salary[j]*1.1 up.append(j)
例に従ってプログラムを組んでみた結果は
すごくギザギザな分布になった。どこかおかしいらしいがまたあす以降に検証したい。